Silent sites in EMBOSS

EMBOSS includes a program called silent that determines positions within a coding sequence which can be mutated to give a convenient restriction site without changing the encoded protein.

$ silent hemA.txt Find restriction sites to insert (mutate) with no translation change Comma separated enzyme list [all]: Output report [hema.silent]: emboss.results.txt

The column headings and a sample line are as shown, when run on the sequence of the Salmonella typhimurium hemA gene:

Start End Strand EnzymeName RS-Pattern Base-Posn AAs Silent Mutation 327 332 + SpeI ACTAGT 330 L.L Yes G->A

In an early post I blogged about this and posted some code.

It seemed like a good idea to compare results from the two programs. One problem is the redundancy of restriction enzymes with respect to cleavage sites. Every position reported by silent in EMBOSS has a number of enzymes, some as many two dozen or so (each one listed as a separate hit).

If you've ever worked with these enzymes you know that they have personalities, some are highly active, some produce overhangs, some are stable, some not so good. Filtering out the redundancy in a smart way is a bit of a pain. That's why one of the optional inputs to the program is a list of enzymes that you like. I ran silent with the default setting (lots of hits), and then for this analysis I made a short list of "good_enzymes."

To repeat what we did previously I copied out Restriction_Dictinary.py from the biopython-1.55 source. Then I went back to the old post and got three files (REnzymes, GeneticCode2, extrasites).

I ran REnzymes.py and it looks like it works even though the format of Restriction_Dictinary.py has changed (and is really unrecognizable to me). I modified the old extrasites.py slightly to make the sequence uppercase. Rather than mess with that script any more, I write the results to disk as my.results.txt.

$ python extrasites.py > my.results.txt

So now the problem is to load to the two different textfiles with results and compare the data. A classic everyday bioinformatics problem. My output looks like this:

codon 123 GCG => GCT AAAAAA GCG TTTGCG AAAAAA GCT TTTGCG HindIII AAGCTT

The original sequence is on the second line and the mutated sequence on the third. The affected codon is set off from the surrounding sequence by a space on each side.

The code to analyze the differences is pretty ugly. I listed it below but I hope you don't look at it unless you're stuck. The output, shown next, reveals that for the most part the two sets of results are congruent. EMBOSS results are output on a single line starting with 'E', while my results are output on three lines, the first starting with 'M'. The results were sorted by order of position in the sequence. The EMBOSS Base-Posn has been converted to a codon for consistency.

A few significant differences can be seen. Mostly they involve a single enzyme, KasI. I haven't sorted that out yet. Overall, I think the agreement is very good.

output:

$ python analyze.py E (40, 13, 'GGCGCC', 'KasI') M --- 13 KasI AAAACG GCA CCTGTA AAAACG GCG CCTGTA E (145, 48, 'GTCGAC', 'SalI') M --- 48 SalI GTGCTG TCA ACCTGT GTGCTG TCG ACCTGT E (199, 66, 'CTGCAG', 'PstI') M --- 66 PstI AACCTG CAA GAAGCG AACCTG CAG GAAGCG E (322, 107, 'TCTAGA', 'XbaI') M --- 107 XbaI AGCGGT CTG GATTCA AGCGGT CTA GATTCA E (331, 110, 'ACTAGT', 'SpeI') M --- 110 SpeI GATTCA CTG GTGCTG GATTCA CTA GTGCTG E (373, 124, 'AAGCTT', 'HindIII') M --- 124 HindIII AAAAAA GCG TTTGCG AAAAAA GCT TTTGCG E (481, 160, 'GGCGCC', 'KasI') M --- 160 KasI ATCGGC GCT AGCGCC ATCGGC GCC AGCGCC E (526, 175, 'AGATCT', 'BglII') M --- 175 BglII GCCCGC CAA ATCTTT GCCCGC CAG ATCTTT E (541, 180, 'GTCGAC', 'SalI') M --- 180 SalI GAATCG CTC TCGACG GAATCG CTG TCGACG M --- 180 SalI GAATCG CTC TCGACG GAATCG TTG TCGACG E (571, 190, 'GGCGCC', 'KasI') E (589, 196, 'ACTAGT', 'SpeI') M --- 196 SpeI ATTGAA CTG GTGGCG ATTGAA CTA GTGGCG E (685, 228, 'GGCGCC', 'KasI') E (719, 239, 'CTGCAG', 'PstI') M --- 240 PstI GCCCGT TTG CAGGAT GCCCGT CTG CAGGAT M --- 248 SalI ATTATC AGT TCGACC ATTATC TCG TCGACC M --- 295 MluI GCGAAC GCT TATCTT GCGAAC GCG TATCTT E (904, 301, 'GTCGAC', 'SalI') M --- 301 SalI AGCGTC GAT GATTTA AGCGTC GAC GATTTA M --- 303 PstI GATGAT TTA CAGAGC GATGAT CTG CAGAGC E (952, 317, 'CTGCAG', 'PstI') E (952, 317, 'CTGCAG', 'PstI') M --- 317 PstI CAGGCT GCG GCAGTA CAGGCT GCA GCAGTA M --- 317 PstI CAGGCT GCG GCAGTA CAGGCT GCT GCAGTA E (1021, 340, 'GGCGCC', 'KasI') M --- 340 KasI GCCCAG GGG GCCAGC GCCCAG GGC GCCAGC E (1130, 376, 'CTGCAG', 'PstI') M --- 377 PstI GCCATC TTG CAGGAT GCCATC CTG CAGGAT M --- 377 PstI GCCATC TTG CAGGAT GCCATC CTG CAGGAT E (1135, 378, 'AGATCT') M --- 378 BglII ATCTTG CAG GATCTG ATCTTG CAA GATCTG

code listing:

import REnzymes RE = REnzymes.REnzymes() good_enzymes = ['PstI', 'XbaI', 'SpeI','HindIII', 'SalI', 'KasI', 'SpeI', 'PstI', 'BglII', 'MluI','Bcl'] def load_data(fn): FH = open(fn,'r') data = FH.read() FH.close() return data.strip() #=================================== # part 1: EMBOSS results data = load_data('emboss.results.txt') data = data.split('\n\n')[2] assert data[:7] == ' Start' eL = list() for e in data.strip().split('\n')[1:]: e = e.split() x, j, flag, enz, seq, i = e[:6] if not enz in good_enzymes: continue if not '+' in flag: continue codon = int(i)/3 i = int(i)+1 eL.append((i,codon, seq, enz)) #=================================== # part 2: my results data = load_data('my.results.txt') data = data.split('\n\n')[:-1] mL = list() for entry in data: lines = entry.split('\n') codon = lines[0].split()[1] # bug in original codon = int(codon) + 1 enz,seq = lines[3].split() if not enz in good_enzymes: continue L = ['---', str(codon), enz] L += [lines[1].strip()] L += [lines[2].strip()] mL.append(L) #=================================== # part 3: show e = eL.pop(0) m = mL.pop(0) while eL or mL: if int(m[1]) < int(e[1]): print 'M', ' '.join(m[:3]) print ' ', '\n '.join(m[3:]) m = mL.pop(0) else: print 'E', e e = eL.pop(0) if not mL: print 'E', e[:3] if not eL: print 'M', ' '.join(m[:3]) print ' ', '\n '.join(m[3:])

EMBOSS: first explorations

This is a quick update on my exploration of EMBOSS. It's not a set of detailed notes to follow or anything like that. I'm just trying to get a feeling for how to use it, and recording my impressions.

As mentioned last time, I can launch a window with the JEMBOSS GUI like so:

$ /usr/local/share/EMBOSS/jemboss/runJemboss.sh

On Mac OS X, it installed into /usr/local/share... I found a tutorial here:

/usr/local/share/EMBOSS/doc/tutorials/emboss_tut.tar.gz

Unpack and grab emboss_tutorial.pdf. Nice short tutorial. We get info about a program like say wossname by one of these:

wossname -opt wossname -help wossname -outfile ~/Desktop/file.txt wossname > help.txt

Even I can figure out that what wossname does is to search for info on a specific topic, say "restriction":

$ wossname Finds programs by keywords in their short description Text to search for, or blank to list all programs: restriction SEARCH FOR 'RESTRICTION' rebaseextract Process the REBASE database for use by restriction enzyme applications recoder Find restriction sites to remove (mutate) with no translation change redata Retrieve information from REBASE restriction enzyme database remap Display restriction enzyme binding sites in a nucleotide sequence restover Find restriction enzymes producing a specific overhang restrict Report restriction enzyme cleavage sites in a nucleotide sequence showseq Displays sequences with features in pretty format silent Find restriction sites to insert (mutate) with no translation change

I put a FASTA-formatted sequence file on my desktop and do:

$ remap example.fasta Display restriction enzyme binding sites in a nucleotide sequence Comma separated enzyme list [all]: Minimum recognition site length [4]: Output file [da247.remap]: EMBOSS An error in remap.c at line 244: Cannot locate enzyme file. Run REBASEEXTRACT

$ REBASEEXTRACT Process the REBASE database for use by restriction enzyme applications REBASE database withrefm file: Error: Input file is required

Google led me to REBASE. I was confused because the example input files shown on a help page here don't match the three files I downloaded. But then I noticed that (some of) the example output files match. REBASE/embossre.enz looks like link_emboss_e.txt.

Actually reading the docs: the input file must be the "withrefm" file of a REBASE distribution. For example, the withrefm file for REBASE version 005 is at: ftp://ftp.neb.com/pub/rebase/withrefm.005. But the link is dead.

Google: site:rebase.neb.com withrefm

No matches. So, I can't find the required input files on the REBASE site.

[UPDATE: the files are there, see the end of the post]

I found what look like equivalents for two of the four output files including REBASE/embossre.enz. The reference file format looks like it's changed. And the fourth file REBASE/embossre.equ isn't there.

Check out the files in

$ ls /usr/local/share/EMBOSS/data/REBASE/ dummyfile embossre.enz embossre.equ embossre.ref embossre.sup

Some look good, like embossre.equ, but some are dummy files, e.g. embossre.enz.

Try this:

sudo mv embossre.enz embossre.enz.old sudo cp ~/Desktop/link_emboss_e.txt /usr/local/share/EMBOSS/data/REBASE/embossre.enz $ remap -sequence example.fasta Display restriction enzyme binding sites in a nucleotide sequence Comma separated enzyme list [all]: Minimum recognition site length [4]: 6 Output file [da247.remap]: results.txt

The result looks like things are working:

DA247 NheI | Cac8I Cac8I | | BmtI | StuI NspI XmnI \ \ \ \ \ \ \ GATGAACGCTAGCGGCAGGCCTAACACATGCAAGTCGAGGGAGAAGCCCTTCGGGGCGGA 10 20 30 40 50 60 ----:----|----:----|----:----|----:----|----:----|----:----| CTACTTGCGATCGCCGTCCGGATTGTGTACGTTCAGCTCCCTCTTCGGGAAGCCCCGCCT / / / / / / / | | NheI | StuI NspI XmnI | Cac8I Cac8I BmtI D E R * R Q A * H M Q V E G E A L R G G M N A S G R P N T C K S R E K P F G A E * T L A A G L T H A S R G R S P S G R K ----:----|----:----|----:----|----:----|----:----|----:----| S S R * R C A * C M C T S P S A R R P P X H V S A A P R V C A L R P L L G E P R I F A L P L G L V H L D L S F G K P A S

So, what I might do is look harder for withrefm, but it seems to work, at least. Not crazy about the output format, they should do something like:

. . . . . 60 GATGAACGCTAGCGGCAGGCCTAACACATGCAAGTCGAGGGAGAAGCCCTTCGGGGCGGA CTACTTGCGATCGCCGTCCGGATTGTGTACGTTCAGCTCCCTCTTCGGGAAGCCCCGCCT

[UPDATE: Should've looked harder. The files I need are right there on the REBASE downloads page (#5 and #31).

$ sudo rebaseextract Password: Process the REBASE database for use by restriction enzyme applications REBASE database withrefm file: ~/Desktop/link_withrefm.txt REBASE database proto file: ~/Desktop/link_proto.txt c-98-236-78-154:Desktop telliott_admin$ cat /usr/local/share/EMBOSS/data/REBASE/embossre.enz ..

Repeat remap and it looks good. It would be nice if the output of enzymes with the number of sites listed the position of each cut.

# Enzymes that cut Frequency Isoschizomers Acc65I 1 Asp718I AflIII 2 ApoI 1 AcsI,XapI BbvII* 1 BpiI,BpuAI,BstV2I,BbsI BmtI 1 BspOI BsaAI 1 BstBAI,Ppu21I ..

EMBOSS

I made a stab at building EMBOSS but failed and succeeded, I think.

According to the adminguide, it requires zlib and libpng as well as gd. Since EMBOSS is set up for Linux, I will probably have to edit the config file to show where these libraries are on my machine.

I installed zlib and libpng to get matplotlib up and running (according to Gavin Huttley's wiki instructions).

Let's get gd. Downloaded gd-2.0.35.tar.gz from this page. According to the readme, it's the usual:

./configure make install

(using sudo). And it has a problem:

gdft.c:1403:35: error: fontconfig/fontconfig.h: No such file or directory gdft.c:1466: error: expected '=', ',', ';', 'asm' or '__attribute__' before '*' token gdft.c:1479: error: expected ')' before '*' token gdft.c: In function 'font_pattern': gdft.c:1516: error: 'FcPattern' undeclared (first use in this function) gdft.c:1516: error: (Each undeclared identifier is reported only once gdft.c:1516: error: for each function it appears in.) gdft.c:1516: error: 'font' undeclared (first use in this function) gdft.c:1517: error: 'FcChar8' undeclared (first use in this function) gdft.c:1517: error: 'file' undeclared (first use in this function) gdft.c:1518: error: 'pattern' undeclared (first use in this function) gdft.c:1525: warning: assignment makes pointer from integer without a cast gdft.c:1531: error: expected ')' before 'FcChar8' gdft.c:1538: error: 'FC_FILE' undeclared (first use in this function) gdft.c:1538: error: 'FcResultMatch' undeclared (first use in this function) make[1]: *** [gdft.lo] Error 1 make: *** [install-recursive] Error 1

It seems to be related to a library called Fontconfig.

Get fontconfig-2.8.0.tar.bz2 from here. According to INSTALL do:

./configure --sysconfdir=/etc --prefix=/usr --mandir=/usr/share/man

but it fails because there is no configure, why not? From a web search it seems that

autoconf processes configure.in to produce a configure script

so I try this:

$ autoconf configure.in:36: error: possibly undefined macro: AM_INIT_AUTOMAKE If this token and others are legitimate, please use m4_pattern_allow. See the Autoconf documentation. configure.in:38: error: possibly undefined macro: AM_MAINTAINER_MODE configure.in:59: error: possibly undefined macro: AM_CONFIG_HEADER configure.in:64: error: possibly undefined macro: AC_LIBTOOL_WIN32_DLL configure.in:65: error: possibly undefined macro: AM_PROG_LIBTOOL configure.in:78: error: possibly undefined macro: AM_CONDITIONAL

and then do as before:

$ sudo ./configure --sysconfdir=/etc --prefix=/usr --mandir=/usr/share/man ./configure: line 1728: syntax error near unexpected token `fontconfig,' ./configure: line 1728: `AM_INIT_AUTOMAKE(fontconfig, 2.8.0)'

It looks like I'm supposed to use autoconf, but first I have to fix the above error.

After another web search, try adding this to configure.in

m4_pattern_allow([AM_INIT_AUTOMAKE]) m4_pattern_allow([AM_MAINTAINER_MODE]) m4_pattern_allow([AM_CONFIG_HEADER]) m4_pattern_allow([AC_LIBTOOL_WIN32_DLL]) m4_pattern_allow([AM_PROG_LIBTOOL]) m4_pattern_allow([AM_CONDITIONAL])

autoconf goes fine. Now:

$ sudo ./configure --sysconfdir=/etc --prefix=/usr --mandir=/usr/share/man ./configure: line 1736: syntax error near unexpected token `fontconfig,' ./configure: line 1736: `AM_INIT_AUTOMAKE(fontconfig, 2.8.0)'

And I'm stuck. You can't feed fontconfig to AM_INIT_AUTOMAKE.No EMBOSS for me, it seems.

[UPDATE: I found instructions for installing GD here. It turns out I already had /usr/X11R6/include/fontconfig. The test of GD looks fine.

Followed simple instructions for EMBOSS and did:

./configure make make install

And after searching around for 10 minutes I finally ran:

/usr/local/share/EMBOSS/jemboss/runJemboss.sh

Now I have a window to play with (graphic at top of post)! ]