PyCogent 13: Tree with colored edges

I'm exploring PyCogent (docs, download page)---first post here.

Continuing from the last post, we use PyCogent to plot the tree. By looking at the source code, I figured out how to color the edges, but I haven't seen how to color the tip labels. It would also be useful to be able to edit the labels.

from cogent import LoadTree from cogent.draw import dendrogram fn = 'temp.tree' t = LoadTree(fn) def f(node): L = ['DQ450530','AJ278451','AF411020'] if node.isTip(): name = node.getTipNames()[0] for e in L: if e in name: return 'blue' return 'red' return 'black' height, width = 500,500 np = dendrogram.ContemporaneousDendrogram(t) np.drawToPDF('temp.pdf', width, height, font_size=14, edge_color_callback=f)

PyCogent 12: simple Phylogenetic Tree

I'm exploring PyCogent (docs, download page)---first post here.

I want to use the sequences we downloaded last time to make a simple phylogenetic tree. Before we do that, however, let's modify the title lines since they are way too long:

>gi|15384333|gb|AF411019.1| Achromobacter xylosoxidans strain AU0665 16S ribosomal RNA gene, partial sequence

Since I don't know how to do that with an Alignment (or even a Sequence Collection) object yet, we'll modify the original file of downloaded sequences and repeat the alignment.

fn = 'temp.fasta' def short(s): stuff,genus,species = s.split()[:3] gb = stuff.split('|')[3].split('.')[0] return '_'.join((gb,genus[0],species)) def modifyNames(fn): FH = open(fn,'r') data = FH.read().strip() FH.close() L = data.split('\n\n') pL = list() for e in L: title,seq = e.strip().split('\n',1) title = short(title) pL.append('>' + title + '\n' + seq) FH = open(fn,'w') FH.write('\n\n'.join(pL)) FH.close() modifyNames(fn)

The next step is to do the alignment. (I showed this code in the other post). Given the alignment, we make the tree like this:

from cogent import LoadSeqs, DNA from cogent.evolve.models import GTR,HKY85 from cogent.phylo import distance, nj def prelim(): fn = 'temp.aln.fasta' seqs = LoadSeqs(fn,format='fasta') model = HKY85() dcalc = distance.EstimateDistances(seqs,model) dcalc.run(show_progress=False) d = dcalc.getPairwiseDistances() tree = nj.nj(d) print '\n' + tree.asciiArt() + '\n' tree.writeToFile('temp.tree') print d.items()[0][0] print d.items()[0][1] prelim()

The output is (the progress indicator is turned off):

/-DQ450530_A_bacterium /edge.1--| | \-AJ278451_A_xylosoxidans /edge.2--| | | /-AJ002809_A_sp. | \edge.0--| -root----| \-EU373389_A_xylosoxidans | |--AF411020_A_xylosoxidans | \-AF411019_A_xylosoxidans ('AF411019_A_xylosoxidans', 'DQ450530_A_bacterium') 0.00403969203526

The last two lines give the first of the pairwise distances. This is a dictionary keyed by a tuple of the the two sequence titles. We'll explore using PyCogent to draw the tree next time. Having written it to disk, we can use our old friend R:

setwd('Desktop') library(ape) t = read.tree('temp.tree') plot(t) plot(t,cex=1.5) axisPhylo() plot(t,cex=1.3) axisPhylo()

The result is shown at the top of the post. There is a slight problem that I haven't bothered to fix. The sequence titles saved in the Tree are quoted, and R prints them as is.

PyCogent 11: start w/multiple sequences

I'm exploring PyCogent (docs, download page)---first post here.

Next up are multiple sequences. In the first part we grab some rRNA sequences from NCBI.

import sys from cogent.db.ncbi import EFetch L = ['AF411020.1','EU373389.1', 'DQ450530.1','AJ278451.1', 'AF411019.1','AJ002809.1'] fn = 'temp.fasta' def test1(): e = EFetch(db='nucleotide', rettype='fasta', id=','.join(L)) s = e.read() #print s FH = open(fn,'w') FH.write(s) FH.close() test1()

In the second part, we use muscle to align them (like before):

from cogent import LoadSeqs, DNA from cogent.app.muscle import align_unaligned_seqs as malign def test2(): seqs = LoadSeqs(filename=fn, aligned=False, format='fasta') print type(seqs) aln = malign(seqs,DNA) aln.writeToFile('temp.aln.fasta') test2()

<class 'cogent.core.alignment.SequenceCollection'>

In the third part we just load the aligned sequences back into Python.

def test3(): fn = 'temp.aln.fasta' sc = LoadSeqs(fn,format='fasta') print sc[:40] print type(sc) test3()

>gi|21436540|emb|AJ278451.1| Achromobacter xylosoxidans subsp. denitrificans partial 16S rRNA gene, strain DSM 30026 (T) ----AGAGTTTGATCATGGCTCAGATTGAACGCTAGCGGG >gi|15384333|gb|AF411019.1| Achromobacter xylosoxidans strain AU0665 16S ribosomal RNA gene, partial sequence ------AGTTTGATCCTGGCTCAGATTGAACGCTAGCGGG >gi|2832590|emb|AJ002809.1| Alcaligenes sp. 16S rRNA gene, isolate R6 ------------------------ATTGAACGCTAGCGGG >gi|92919431|gb|DQ450530.1| Alcaligenaceae bacterium LBM 16S ribosomal RNA gene, partial sequence -ATTAGAGTTTGATCCTGGCTCAGATTGAACGCTAGCGGG >gi|171191189|gb|EU373389.1| Achromobacter xylosoxidans strain TPL14 16S ribosomal RNA gene, partial sequence TCGGAGAGTTTGATCCTGGCTCAGATTGAACGCTAGCGGG >gi|15384334|gb|AF411020.1| Achromobacter xylosoxidans strain AU1011 16S ribosomal RNA gene, partial sequence ------AGTTTGATCCTGGCTCAGATTGAACGCTAGCGGG <class 'cogent.core.alignment.Alignment'>

Notice that when we load unaligned sequences, we have a cogent.core.alignment.SequenceCollection, while an aligned set is a cogent.core.alignment.Alignment. We'll use these going forward. And yes, those are mighty long title lines. And we can do slices on the alignment.

PyCogent 10: basic DNA methods

I'm exploring PyCogent (docs, download page)---first post here.

We'll use this method from the last post (and you'll also need to load the data from your favorite FASTA-formatted sequence file), as shown in that example.

from cogent import Sequence, DNA # method 1: from a string seq = Sequence(moltype=DNA, name=n, seq=s) print '1', seq[:24] print type(seq) print seq[:36].toFasta()

<class 'cogent.core.sequence.DnaSequence'> >SThemA ATGACCCTTTTAGCGCTCGGTATTAACCATAAAACG

print seq.Name print seq[:36] print seq[:36].complement() print seq.MW() # this strand only

SThemA ATGACCCTTTTAGCGCTCGGTATTAACCATAAAACG TACTGGGAAAATCGCGAGCCATAATTGGTATTTTGC 365443.89

print seq[:36].rc() # reverse complement print seq.canPair(seq.rc())

CGTTTTATGGTTAATACCGAGCGCTAAAAGGGTCAT True

print seq[-36:].toRna() print seq.hasTerminalStop()

CUGAAUAUUCUGCGCGACAGCCUCGGGCUGGAGUAG True

s = seq.withoutTerminalStopCodon() print s[-33:] print s.getTranslation()[-11:] print s.getOrfPositions()

CTGAATATTCTGCGCGACAGCCTCGGGCTGGAG LNILRDSLGLE [(0, 1254)]

s = seq.toFasta().split('\n',1)[1] print s[:36] print type(s)

ATGACCCTTTTAGCGCTCGGTATTAACCATAAAACG <type 'str'>

mseq = seq.shuffle() print seq[:36] print mseq[:36] print seq.count('g') == mseq.count('g') print all([seq.count(n) == mseq.count(n) for n in 'acgt'])

ATGACCCTTTTAGCGCTCGGTATTAACCATAAAACG GGTCCAATGTCCCAGCGACTGAGCTGTGAACGTTAA True True

PyCogent 9: getting serious with Sequence

I'm exploring PyCogent (docs, download page)---first post here.

The most central and basic object type in PyCogent is Sequence. There are several ways to construct a DNA Sequence object. For the demo, I need to load a FASTA-formatted sequence, and recover the title and DNA sequence separately. (Substitute your own to follow along).

from cogent import LoadSeqs, DNA, Sequence def loadData(fn): FH = open(fn, 'r') data = FH.read() FH.close() return data fn = 'cogent/SThemA.fasta' data = loadData(fn) n,s = data.strip().split('\n',1) n = n[1:]

Now, here are five ways to create a DNA Sequence object:

# method 1: from a string seq = Sequence(moltype=DNA, name=n, seq=s) print '1', seq[:24] # method 2: from a text file tfn = 'cogent/SThemA.txt' seq = Sequence(moltype=DNA, filename=tfn, format='fasta') print '2', seq[:24] # method 3: from a file # guess format from extension seq = Sequence(moltype=DNA, filename=fn) print '3', seq[:24] # method 4: a method in DNA seq = DNA.makeSequence(s) seq.Name = n print '4', seq[:24] # method 5: by joining two DNA sequences a = DNA.makeSequence("AGTACACTGGT") seq2 = seq[:6] + a print '5', seq2 # method 6: from an NCBI record (via str) from cogent.db.ncbi import EFetch e = EFetch(db='nucleotide', rettype='fasta', id='154102') data = e.read() n,s = data.strip().split('\n',1) # go to method 1 seq = Sequence(moltype=DNA, name=n, seq=s) print '6', seq[:24]

output:

1 ATGACCCTTTTAGCGCTCGGTATT 2 ATGACCCTTTTAGCGCTCGGTATT 3 ATGACCCTTTTAGCGCTCGGTATT 4 ATGACCCTTTTAGCGCTCGGTATT 5 ATGACCAGTACACTGGT 6 ATGACCCTTTTAGCGCTCGGTATT

Here is one more that is a little trickier:

# method 7: using a parser from cogent.parse.fasta import FastaParser # must modify for lowercase f=open(fn) g = FastaParser(f) n,s = g.next() s = str(s).upper() # go to method 1 seq = Sequence(moltype=DNA, name=n, seq=s) print '7', seq[:24] # simpler use is uppercase sequence fn2 = '.temp' FH = open(fn2,'w') FH.write('>' + n + '\n' + s.upper()) FH.close() f=open(fn2) g = FastaParser(f) n,s = g.next() # go to method 1 seq = Sequence(moltype=DNA, name=n, seq=s) print '8', seq[:24]

7 ATGACCCTTTTAGCGCTCGGTATT 8 ATGACCCTTTTAGCGCTCGGTATT

We'll look at the most basic methods of DNA Sequence objects next.

elementary $PATH stuff

As I mentioned the other day, it seems like PyCogent likes to use the $PATH variable rather than Python's sys.path to look for executables like clustalw and muscle. While I was updating my clustalw (actually clustalw2) installation at work to version 2.0.12, I worked on this quite elementary BASH shell issue. (I got clustalw from EMBL-EBI, see FTP download link here).

There is a directory in my home directory ~/bin where my bioinformatics binaries live. I copied the new folder clustalw-2.0.12-macosx from the download into there, as usual. Then, I just made a symbolic link to the executable like so:

ln -s ~/bin/clustalw-2.0.12-macosx/clustalw2 ~/bin/clustalw

I used clustalw for the link since that's what Py-Cogent expects. I tested it by doing ~/bin/clustalw from the Desktop. Now I needed to put ~/bin on my PATH. You can look at the $PATH variable as follows:

echo $PATH

output before the modification:

/usr/bin:/bin:/usr/sbin:/sbin:/usr/local/bin:/usr/X11/bin

To modify it, I made a text file called on my Desktop containing this line:

PATH=$PATH:$HOME/bin/ ;export PATH

I put it in the right place by doing:

cp ~/Desktop/x.txt ~/.bash_profile

and tested by first doing

echo $PATH

which gives /Users/te/bin/ appended to the previous output. The final test (for now, anyway) is to do clustalw from the Desktop.

$ clustalw ************************************************************** ******** CLUSTAL 2.0.12 Multiple Sequence Alignments ******** ************************************************************** 1. Sequence Input From Disc 2. Multiple Alignments 3. Profile / Structure Alignments 4. Phylogenetic trees S. Execute a system command H. HELP X. EXIT (leave program)

PyCogent 8: getting matplotlib

I'm exploring PyCogent (docs, download page)---first post here.

As I mentioned the other day, I was having a little trouble with the drawing code for case study #2 from the paper (PMID 17708774, link). (Note that this part of the example is not in the text, but in supplementary file #4, link).

Although PyCogent initially used ReportLab for at least some drawing, I'm told that they have switched over completely to matplotlib. It's a good thing I didn't know that at first. I've struggled with matplotlib in the past. They have a special page about installing on OS X (the "it sucks to be you" kind of page), where they strongly recommend that you install a separate Python.

The senior author on the PyCogent paper is Gavin Huttley, whose group is at the John Curtin School of Medical Reearch at ANU in Canberra, Australia. I posted to the PyCogent help forum at sourceforge, and Gavin very patiently walked me through the solution to the issues I encountered. Most important, he worked up a wiki page explaining in detail how to build and install matplotlib for OS X Snow Leopard.

Another problem that slowed me down was that I had installed PyCogent first with pip, and subsequently with subversion, and somehow screwed things up. I ripped PyCogent out and rebuilt it and it works fine now. See the help thread for details.

Now, in the directory holding all the files I got from the supplementary data for the paper, I do:

cd src python rRNA_display_tree.py

where the script has been modified slightly as described in the thread. A pdf is written to the figs directory. Here is a screenshot of the upper-left hand corner.



As described in the paper:

We display low G+C% to high G+C% on a spectrum from yellow to blue...We note that, in general, closely related taxa exhibit a similar color. However, certain lineages appear to have evolved to low G+C%, quickly raising questions about environmental and/or changes in DNA metabolism that may distinguish these organisms from their sister taxa. Another feature apparent from this tree is a general trend for earlier diverging lineages to be intermediate in G+C%, and for sud- den changes toward low G+C% to be more common than sud- den changes to high G+C% (there are more yellow branches surrounded by blue or green neighbors than blue branches surrounded by yellow or green neighbors).

I don't get much out of knowing the identity of the organism that has gone to low GC recently. But here is a paper (PMID 19001264) from Dan Andersson's group that may be relevant (ultimately) as to mechanism.

Python package installations

This post is really for me, in the sense of "letting google organize my head." I'm getting into PyCogent (first post here). They're switching over to using matplotlib, which has complicated installation instructions for OS X. This has always scared me away before. But, I want PyCogent to work. Before I get into that, I need to re-learn what Python installation tools are about.

I'm sure any Pythonista knows all this already, but here goes...

What is easy_install?

Easy Install comes from PEAK (Python Enterprise Application Kit). The docs are here. How did it get on my Snow Leopard system? According to the PyCogent docs front page, it comes with the stock Python install on Snow Leopard. I couldn't find it by searching the system. My site-packages directory /Library/Python/2.6/site-packages has a file named easy-install.pth but that's it. easy_install must be somewhere because this works:

easy_install --help

In particular the output of that command includes:

--upgrade (-U) force upgrade (searches PyPI for latest versions) usage: easy_install [options] requirement_or_url ...

According to the docs easy_install is related to setuptools (see below), like a "wrapper" for them. According to the Pylons docs, it is actually part of the setuptools module. What about that easy-install.pth file?

According to the same page:

Easy Install works by maintaining a file in your Python installation’s lib/site-packages directory called easy_install.pth. Python looks in .pth files when it is starting up and adds any module ZIP files or paths specified to its internal sys.path. Easy Install simply maintains the easy_install.pth file so that Python can find the packages you have installed with Easy Install.

The first four lines of my easy-install.pth are:

import sys; sys.__plen = len(sys.path) ./numpy-1.3.0-py2.6-macosx-10.6-universal.egg ./pip-0.6.1-py2.6.egg ./Sphinx-0.6.3-py2.6.egg

Now, I think numpy-1.3 is stock on Snow Leopard. The time stamp on the referenced file is Friday, September 4, 2009 8:38 AM. I think that is the date of my upgrade to Snow Leopard, so it checks. easy_install probably comes with Snow Leopard and is used to do the numpy install.

What is setuptools?

Python code to "build, install, upgrade and uninstall Python package - easily!" It comes from here. It is what we're using when we do:

python setup.py build sudo python setup.py install

So what is PyPi? It's the Python Package Index, "a repository of software for the Python programming language…currently 8369 packages." It either is or used to be named CheeseShop (a Monty Python thing). According to this tutorial, easy_install "gives you a quick and painless way to install packages remotely." Ah... now I get it.

What is pip?

The PyCogent Quick Installation instructions mention this. According to PyPi pip is intended as a replacement for easy_install. Among other things, pip has explicit requirement files. For example, cogent-requirements.txt:

cogent Sphinx>=0.6 numpy>=1.3.0 http://www.reportlab.org/ftp/ReportLab_2_3.tar.gz MySQL-python>=1.2.2 SQLAlchemy>=0.5 http://jcsmr.anu.edu.au/org/genbio/compgen/PyxMPI.tar.gz

Note that Apple has provided a special place for /site-packages to live. It is not where you might expect:
/System/Library/Frameworks/Python.framework/Versions/2.6/lib/python2.6

(although there is a lot there) but rather here:
/Library/Python/2.6/site-packages
You'll can see this in sys.path.

PyCogent 7: Muscle

More on PyCogent here. The objective for this time is to use the Application Controller for Muscle, and construct a phylogenetic tree.

As with the previous example using clustal, it's important to be sure the the binary is on my path. I did this from the Desktop:

ln -s /Users/te/bin/muscle3.6_src/muscle ~/bin/muscle export PATH=/Users/te/bin:$PATH

Now, we expect to run muscle as we did clustal, but there is a slight problem:

from cogent.app.muscle import Muscle def test1(): app = Muscle() result = app('opuntia.fasta') print result for k in ['StdErr','StdOut']: print k print result[k].read() test1()

{'ExitStatus': 0, 'StdErr': , 'StdOut': } StdErr MUSCLE v3.6 by Robert C. Edgar http://www.drive5.com/muscle This software is donated to the public domain. Please cite: Edgar, R.C. Nucleic Acids Res 32(5), 1792-97. opuntia 7 seqs, max length 902, avg length 896 00:00:00 10 MB(2%) Iter 1 100.00% K-mer dist pass 1 00:00:00 10 MB(2%) Iter 1 100.00% K-mer dist pass 2 00:00:01 10 MB(2%) Iter 1 100.00% Align node 00:00:01 10 MB(2%) Iter 1 100.00% Root alignment 00:00:01 10 MB(2%) Iter 2 100.00% Refine tree 00:00:01 10 MB(2%) Iter 2 100.00% Root alignment 00:00:01 10 MB(2%) Iter 2 100.00% Root alignment 00:00:02 10 MB(2%) Iter 3 100.00% Refine biparts 00:00:03 10 MB(2%) Iter 4 100.00% Refine biparts StdOut

The result of the app() call is supposed to be a dictionary-like object with a key 'Align' but that key is not present. Of the two file-like objects that are present, 'StdErr' shows that the program ran, but there aren't any results. This is probably an error, since we're supposed to have a "consistent API."

To work around this issue, we go to an example in the online docs, and try the following instead. I've changed the sequence file to test.fasta, which contains five bacterial 16S rRNA sequences.

from cogent import LoadSeqs, DNA from cogent.app.muscle import align_unaligned_seqs as malign def test2(): seqs = LoadSeqs(filename='test.fasta', aligned=False) #print seqs aln = malign(seqs,DNA) #print aln aln.writeToFile('test.aln.fasta') test2()

The file ends up where it should be, with the first two lines:

>Actinomyces_naeslundii ----------TTGATCCTGGCTCAGGACGAACGCTGGCGGCGTGCTTAACACATGCAAGT

The next step is to make a tree.

from cogent.evolve.models import GTR from cogent.phylo import distance, nj def test3(): seqs = LoadSeqs(filename='test.aln.fasta', aligned=True) model = GTR() dcalc = distance.EstimateDistances(seqs,model) dcalc.run(show_progress=True) d = dcalc.getPairwiseDistances() tree = nj.nj(d) print tree.asciiArt() tree.writeToFile('test.tree') test3()

Doing [1/10]: Actinomyces_naeslundii <-> Achromobacter_xylosoxidans Doing [2/10]: Actinomyces_naeslundii <-> Abiotrophia_defectiva Doing [3/10]: Actinomyces_naeslundii <-> Aeromonas_hydrophila Doing [4/10]: Actinomyces_naeslundii <-> Acinetobacter_haemolyticus Doing [5/10]: Achromobacter_xylosoxidans <-> Abiotrophia_defectiva Doing [6/10]: Achromobacter_xylosoxidans <-> Aeromonas_hydrophila Doing [7/10]: Achromobacter_xylosoxidans <-> Acinetobacter_haemolyticus Doing [8/10]: Abiotrophia_defectiva <-> Aeromonas_hydrophila Doing [9/10]: Abiotrophia_defectiva <-> Acinetobacter_haemolyticus Doing [10/10]: Aeromonas_hydrophila <-> Acinetobacter_haemolyticus /-Actinomyces_naeslundii /edge.0--| | \-Abiotrophia_defectiva | -root----| /-Aeromonas_hydrophila |-edge.1--| | \-Acinetobacter_haemolyticus | \-Achromobacter_xylosoxidans

For the last part, I downloaded the code from the supplementary files for the PyCogent paper. I won't show the code because, unfortunately, it doesn't work.

from matplotlib.path import Path ImportError: No module named path

Although the paper says it uses ReportLab, and that is what I installed the other day, the code is looking for matplotlib, which I don't have. But still, we're making progress.

PyCogent 6: NCBI

More on PyCogent here. Here are three examples from the cookbook.

• EFetch to obtain a nucleotide sequence from the id
• EUtils to obtain a medline report from Pubmed
• Using a parser to process a Genbank file from disk

import sys from cogent.db.ncbi import EFetch from cogent.db.ncbi import EUtils from cogent.parse.genbank import parse def test1(): e = EFetch(db='nucleotide', rettype='fasta', id='154102') result = e.read() L = result.split('\n') print L[0][:40] print L[1][:40] test1()

>gi|154102|gb|J04243.1|STYHEMAPRF S.typh GGATCCACTGCCGCAGGCTGTTTAACGGAATCGGCATCCC

def test2(): e = EUtils(db='pubmed', rettype='medline') item = e['2544564'] # PMID result = item.read() print result.split('\n')[0] test2()

PMID- 2544564

This requires a BioPython file from here.

def test3(): gb = parse(open('ls_orchid.gbk')) print gb.Name #print gb.Info L = gb.Info['features'] L = [e for e in L if e['type'] == 'gene'] print L[0]['gene'][0]

Z78533 5.8S rRNA

Uncomment the line to see all of what is available.

PyCogent 5: Alignment

More on PyCogent here. This is short and sweet.

Two alignments by Needleman-Wunsch (global) and Smith-Waterman (local).

from cogent.align.algorithm import nw_align from cogent.align.algorithm import sw_align seq1 = 'TTCATA' seq2 = 'TGCTCGTA' print nw_align(seq1,seq2) print sw_align(seq1,seq2)

('T--TCATA', 'TGCTCGTA') ('TCATA', 'TCGTA')

PyCogent 4: Clustalw

I'm continuing with an exploration of PyCogent (first post here). The goal for this post is to run clustalw using an "Application Controller" from PyCogent. The example from the cookbook uses an input file from Biopython (download link), opuntia.fasta, which is on my Desktop.

I ran into two issues that were solved pretty easily. The first is that I don't usually manipulate the $PATH variable from the shell, but rather specify needed paths explicitly or add them to sys.path from within Python. However, PyCogent does not seem to look at sys.path:

import sys p = '/Users/te/bin/clustalw-2.0.10-macosx' sys.path.insert(0,p)

When I run my script (below), I get this:

cogent.app.util.ApplicationNotFoundError

So, I did this from the command line before running my test script:

export PATH=/Users/te/bin/clustalw-2.0.10-macosx:$PATH

The other issue is that the Clustal app controller expects clustalw. I have clustalw2. I don't know how to fix the App Controller, so I made a symbolic link in that directory:

ln -s clustalw2 clustalw

Come to think of it, I could have just done that from the Desktop and solved both issues at once. I held my breath, and did this:

from cogent.app.clustalw import Clustalw app = Clustalw() result = app('opuntia.fasta')['Align'] print result print ''.join(result.readlines()[:12])

<open file 'opuntia.aln', mode 'r' at 0x101471690> CLUSTAL 2.0.10 multiple sequence alignment gi|6273285|gb|AF191659.1|AF191 TATACATTAAAGAAGGGGGATGCGGATAAATGGAAAGGCGAAAGAAAGAA gi|6273284|gb|AF191658.1|AF191 TATACATTAAAGAAGGGGGATGCGGATAAATGGAAAGGCGAAAGAAAGAA gi|6273287|gb|AF191661.1|AF191 TATACATTAAAGAAGGGGGATGCGGATAAATGGAAAGGCGAAAGAAAGAA gi|6273286|gb|AF191660.1|AF191 TATACATAAAAGAAGGGGGATGCGGATAAATGGAAAGGCGAAAGAAAGAA gi|6273290|gb|AF191664.1|AF191 TATACATTAAAGGAGGGGGATGCGGATAAATGGAAAGGCGAAAGAAAGAA gi|6273289|gb|AF191663.1|AF191 TATACATTAAAGGAGGGGGATGCGGATAAATGGAAAGGCGAAAGAAAGAA gi|6273291|gb|AF191665.1|AF191 TATACATTAAAGGAGGGGGATGCGGATAAATGGAAAGGCGAAAGAAAGAA ******* **** *************************************

Note: as usual, all code is executed from the command line by doing
python script.py.

PyCogent 3: Translation

I'm continuing with an exploration of PyCogent (first post here). I found a draft cookbook on the web. That will be a big help!

The problem that I had in following the first tutorial is this: the function LoadSeqs creates an object of this class:

<class 'cogent.core.alignment.SequenceCollection'>

What we want is a single DNA sequence. We can get the individual sequences back from the collection like so:

from cogent import LoadSeqs, DNA fn = 'SThemA.fasta' s = LoadSeqs(fn,moltype=DNA,aligned=False) name,seq = s.items()[0]

seq is now a

<class 'cogent.core.sequence.DnaSequence'>

and it supports slicing. We can also construct a similar object directly:

seq2 = DNA.makeSequence("AGTACACTGGT") seq2.Name = 'X' print seq2.toFasta()

>X AGTACACTGGT

We still can't translate a stop codon:

pep = seq.getTranslation()

cogent.core.alphabet.AlphabetError: TAG

But slicing makes it easy to fix that:

pep = seq[:-3].getTranslation() print type(pep) print pep[:24]

<class 'cogent.core.sequence.ProteinSequence'> MTLLALGINHKTAPVSLRERVTFS

Sequences are complex objects:

L = dir(seq) L = [e for e in L if not e[0] == '_'] for e in L[:10]: print e print len(L), 'items'

CodonAlphabet Info LineWrap MW MolType Name PROTEIN addAnnotation addFeature annotateFromGff 81 items

Going through the dir reveals the intended method for translating with a terminal stop:

seq = DNA.makeSequence("ATGACACTGGTGTAG") seq.Name = 'X' print seq.toFasta() print seq.hasTerminalStop() seq2 = seq.withoutTerminalStopCodon() print seq2.getTranslation() # reverse complement print seq.rc().canPair(seq)

>X ATGACACTGGTGTAG True MTLV True

PyCogent (2)

I've started playing with PyCogent. The first simple tutorial is here. We start by importing stuff:

from cogent import LoadSeqs, DNA fn = 'SThemA.txt' s = LoadSeqs(fn)

This fails, even though the file is in FASTA format.

cogent.parse.record.FileFormatError: Unsupported file format txt

We must change the extension.

fn = 'SThemA.fasta' s = LoadSeqs(fn) print s.MolType

Now, this kinda sorta works, but there is a slight problem:

MolType(('a', 'b', 'c', 'd', 'e', 'f', 'g', 'h', 'i', 'j', 'k', 'l', 'm', 'n', 'o', 'p', 'q', 'r', 's', 't', 'u', 'v', 'w', 'x', 'y', 'z', 'A', 'B', 'C', 'D', 'E', 'F', 'G', 'H', 'I', 'J', 'K', 'L', 'M', 'N', 'O', 'P', 'Q', 'R', 'S', 'T', 'U', 'V', 'W', 'X', 'Y', 'Z'))

We must tell it what type of characters we're dealing with.

fn = 'SThemA.fasta' s = LoadSeqs(fn,moltype=DNA,aligned=False) print s.MolType print s.SeqLen

That's much better:

MolType(('T', 'C', 'A', 'G')) 1257

Let's print the sequence. The commented versions don't work:

# print s[:24] # print s.toDna()[:24] print s.toFasta()[:24]

>ST hemA ATGACCCTTTTAGCG

The error with the first line is:

TypeError: 'SequenceCollection' object is unsubscriptable

Finally, let's translate:

print s.getTranslation()

>ST hemA MTLLALGINHKTAPVSLRERVTFSPDTLDQALDSLLAQPMVQGGVVLSTCNRTELYLSVEEQDNLQEALIRWLCDYHNLNEDDLRNSLYWHQDNDAVSHLMRVASGLDSLVLGEPQILGQVKKAFADSQKGHLNASALERMFQKSFSVAKRVRTETDIGASAVSVAFAACTLARQIFESLSTVTVLLVGAGETIELVARHLREHKVQKMIIANRTRERAQALADEVGAEVISLSDIDARLQDADIIISSTASPLPIIGKGMVERALKSRRNQPMLLVDIAVPRDVEPEVGKLANAYLYSVDDLQSIISHNLAQRQAAAVEAETIVEQEASEFMAWLRAQGASETIREYRSQSEQIRDELTTKALSALQQGGDAQAILQDLAWKLTNRLIHAPTKSLQQAARDGDDERLNILRDSLGLE

But, it only works if I trim off the stop codon in the file. I have no idea how to do this from code, since the objects don't support slicing.

print s.__dict__

{'Info': {'Refs': {}}, '_seqs': [DnaSequence(ATGACCC... 1254)], 'Name': None, 'NamedSeqs': {'ST hemA': DnaSequence(ATGACCC... 1254)}, 'Alphabet': ('T', 'C', 'A', 'G', '-', 'B', 'D', 'H', 'K', 'M', 'N', 'S', 'R', 'W', 'V', 'Y', '?'), 'SeqLen': 1254, 'Names': ['ST hemA'], 'SeqData': [DnaSequence(ATGACCC... 1254)], 'MolType': MolType(('T', 'C', 'A', 'G'))}

Any ideas?

What we can do is manipulate the data before making the object. Now, we can even use a text file!

FH = open('SThemA.txt') data = FH.read().strip() data = data.split('\n',1)[1][:-3] seqs = { 'SThemA':data } s = LoadSeqs(data=seqs, moltype=DNA,aligned=False) print s.getTranslation().toFasta()[:24]

>ST hemA MTLLALGINHKTAPV

We had to split off the title line first, then lose the last three bases. That is good enough for now.

PyCogent: baby steps

Rob Knight (at the University of Colorado) has been on my radar for a while now because of his involvement in a number of metagenomics projects, like with Jeff Gordon. Besides a solid statistics and phylogenetics background, he brings an enthusiasm for Python to the mix, and I find that encouraging since I can at least hope to figure out what he's up to. I have already used UniFrac for a bacterial census project I'm involved in, and I mean to post about that sometime later---like when we have actually written the paper!

Rob is also a developer of PyCogent (PMID 17708774), billed as a "toolkit for making sense from sequence." I read the paper a year ago, but didn't get into it at the time because it seemed like a pretty complex undertaking.

Now, my desire has become more pressing because of a new paper that is coming out from Rob which integrates the PyCogent framework with NAST, called PyNast (PMID 19914921). I want to learn how to use template alignments.

I followed the instructions for "Quick installation" of PyCogent here. I ran into a problem with MySQL that I'm still working on. And there were other apparent errors in the pip-log file that led me to go back and try to re-install each component one-at-a-time. This seemed to work (except for MySQL, which needs mysql_config), but when I looked for tests they were not present. It seems that these don't end up in the site-packages directories. Since I wasn't confident that I'd actually installed PyCogent correctly, I used subversion to grab the source and do it again:

svn co https://pycogent.svn.sourceforge.net/svnroot/pycogent/trunk PyCogent

and followed the instructions in the README. That seems to have worked since this:

sh run_tests

gave:

[snip] Can't find blastall executable: skipping test [snip] ... [snip] ... --- [snip] --- Ran 3403 tests in 185.362s OK

Obviously, I will have to let PyCogent know where executables live, but otherwise, that looks reassuring. However, the data file for the very first tutorial is not present in what I have on my machine. Seems there are a few rough spots to work through. I'll keep you posted.

Politically incorrect, I suppose

I promise that we won't do much of this, but as long as we're being personal for the moment (see the last post), I can't help writing a bit about where I am politically. My favorite blogs are:

Steve Benen
Ezra Klein
Matthew Iglesias
Kevin Drum
James Fallows
Tom Ricks
Brad DeLong
Stephen Walt
Bernard Avishai

I watched President Obama's speech on Afghanistan the other night on TV. Politically, I think he has to do this extra something, though he doesn't want to, and I don't either. If he failed to do it, the right wing would cream him.

I get discouraged because I think what we are looking at here, ultimately, is a "lost decade" or even more. The correct response to 9/11 should have been, in order:

• lock the cockpit doors on commercial aircraft and alter our response to threatened hijackings

• think about arming marshals and putting them aboard

• do bad things to Osama what'hoosits and his buddy KSM---these would involve helicopters at an altitude of 1350 over the WTC site if I had my way...

Instead we spent 1 trillion dollars (or is it 2?) on an unnecessary war and the aftermath. I know we're a rich country and all, but this is totally crazy. Don't these guys read history?

As Bronowski says near the "Ascent of Man" (a favorite):

"It seems very pessimistic to talk about western civilization with a sense of retreat. I have been so optimitistic about the ascent of man; am I going to give up at this moment? Of course not. The ascent of man will go on. But do not assume that it will go on carried by western civiliation as we know it. We are being weighed in the balance at this moment. If we give up, the next step will be taken---but not by us. We have not been given any guarantee that Assyria and Egypt and Rome were not given. We are waiting to be somebody's past too, and not necessarily that of our future."

Oh, the horror!

Maybe this should be more like a real blog. At least once in a while.

One of the formative experiences of my youth was seeing a live performance of The Rocky Horror Show, in a West L.A. venue that is lost to time. It was fantastic! The stage show had all of the good things in the movie, especially the songs, and none of the crap. (Though no Susan Sarandon). Tim Curry was overpowering. Why did life not give him more of a career?

In that spirit, check out Richard O'Brien doing the intro---"Science Fiction, Double Feature" on youtube. Or buy it from iTunes, it's definitely worth it.

And you can explore all the references from the song at wikipedia.

Man, those were the days. I suppose it's just a coincidence that I also downloaded "All the young dudes" today.

So it's kind of weird. Today, my generation look so... so... (bourgeois isn't quite what I'm looking for), but all that stuff is still inside my head. And the next question is, what was in the elder generation's head when Pete Townshend sang "hope I die before I get old" ?? Were they wild and crazy too? I thought we were special.

In that spirit, here is a friend of mine.

Boggle

I found a fun text problem here, called Boggle. It is even an app for the iPhone. The basic rules are that we have a 4 x 4 matrix containing letters

asta
ifue
mtnd
aaos

Given a puzzle, the object is to form words, starting at any position but using two rules:

• sequential letters must be neighbors (diagonals are OK)
• no position can be used twice

(duplicate letters are valid puzzles, as shown in the example)

To make life simpler, we can use a dictionary of valid words, for example (on my Snow Leopard system):

/usr/share/dict/words

Rather than choose letters for the puzzle equally often, I chose letters at the frequency they occur in the dictionary. Here is one run:

234936 asta ifue mtnd aaos 954 amandus dentist donatist itonama mantodea mistend situate situated student tamandu tamandua tsunami

There are almost 235,000 words in the dictionary. 954 of them pass the first filter (below). We only print words longer than a threshold T = 5.

There are lots of solutions shown, in different languages (link).

My version uses an idea adapted from BLAST. Make a list of all k-words in the table (adjacent letters, k=2), and filter the dictionary against the list. This reduces the search space by about 200- to 1000-fold.

# get k-words (k=2), by position def getKWords(puzzle): L = list() for i in range(16): temp = list() c1 = puzzle[i] for j in near[i]: temp.append(c1 + puzzle[j]) L.append(temp) # convert to simple list kL = list() for eL in L: kL.extend(eL) return list(set(kL))

The second phase tests whether a word is a valid solution to the puzzle by checking its "path." I won't post that, but you can find my full solution here.

Grouping list elements

I was browsing Stack Overflow again. I really like the site, and I'm using it to improve my knowledge of Python. The rating system for questions and answers (and answerers) is really helpful.

But sometimes I have to wonder... For example, here is a question from yesterday about how to group elements from a sequence (and since it's from yesterday, it is naturally a duplicate). The answer that people like has me scratching my head. Here is how I do it:

from string import uppercase as u L = list(u) def groupby(L,N=3): R = range(0,len(L),N) return [L[i:i+N] for i in R] L2 = groupby(L) print L2[:1] # prints: [['A', 'B', 'C']]

I like this because: (a) it works, and (b) it's simple enough to understand at a glance. For very long sequences you might want to use xrange within the list comprehension. Here is the popular answer:

>>> from itertools import izip_longest >>> L=[1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11] >>> list(izip_longest(*[iter(L)]*3)) [(1, 2, 3), (4, 5, 6), (7, 8, 9), (10, 11, None)]

This is just crazy. It is not at all clear how this works or what it does. It would have to be commented in code. And, how does it work? The docs are here. We first make a list of iter(L) objects:

>>> [iter(L)]*2 [<listiterator object at 0x10048b810>, <listiterator object at 0x10048b810>]

But when I try to unpack it with *, I get:

>>> *[iter(L)*3] File "", line 1 *[iter(L)*3] ^ SyntaxError: invalid syntax

So the bottom line is that I don't really know how it works because I can't take it apart.

What izip_longest does is take a list of lists and pop items off each one in turn to put into the groups. From the docs:

izip_longest('ABCD', 'xy', fillvalue='-')
--> Ax By C- D-

It seems like a lot of extra work is being done here, to go along with the obfuscation. These guys should go back to Perl!

Roman numerals - update

After I posted about converting between decimal and Roman numerals here, I came upon a simpler way to handle the forward direction, from decimal to Roman.

Rather than divide sequentially by 1000, 100, and finally 10, we include 900, 500, 400 etc. in our list of divisors, and modify the symbols to be saved accordingly. This idea is included in the Python Cookbook. I'm not sure who thought it up first, but it's a nice trick.

Here is the code:

def decimalToRoman2(n): divs = [1000, 900, 500, 400, 100, 90, 50, 40, 10, 9, 5, 4, 1] symbols = ['M', 'CM', 'D', 'CD', 'C', 'XC', 'L', 'XL', 'X', 'IX', 'V', 'IV', 'I'] vals = list() for i,d in enumerate(divs): x = n/d if x: vals.append(symbols[i]*x) n = n%d return ''.join(vals)