Burroughs-Wheeler Transform-BWT (3)

I'm exploring the Burroughs-Wheeler Transform, which has become important for fast methods of matching short reads to large genomes. The first two posts are here and here.

In the first we saw what BWT is and showed how to reverse or invert it. In the second, we saw how to use BWT to do a search. Here, I show a fairly efficient way to generate the BWT for a sequence.

Last time, we used a script which does a rotation on a deque followed by sorting to generate the BWT of a sequence:

GAATTCAAGCTTGGATCCGGAAAGATCTGATC

and obtained

GACGAAGGGATTCTGTG^GATACTTAAACTACC

The result matches what the visualization tool on this page generates. However, our script does a lot more work (and uses more memory) than is really necessary. If we look again at the first few lines of output from bwt.py:

20 AAAGATCTGATC^GAATTCAAGCTTGGATCCGG 21 AAGATCTGATC^GAATTCAAGCTTGGATCCGGA 6 AAGCTTGGATCCGGAAAGATCTGATC^GAATTC

Each number is the index into the sequence where the row (the suffix) begins (0-based indexing). The array of numbers from top to bottom of the output is the suffix array.

As you can see from the last part of the script below, the BWT is easily generated from the suffix array, it's just text[j-1] for each j in the array. That's because the BWT is the last column and the sequences wrap, so each character in the BWT is -1 with respect to the character in the first column for any row.

The Python script shown here sorts the rows but only looks at as many of the characters as are needed to do the sort. We have a list L to contain the growing suffix array. This consists of indexes j into the text, one corresponding to each row, which are in sorted order in the suffix array.

We move an index i along the text and for the suffix starting at each i, find the place where it belongs in the suffix array, then insert it.

I'm sure there are even faster ways to generate the BWT for a text. The most glaring problem with this attempt is that repeated insertion into a list is expensive. We need a data structure that combines ease of insertion with ease of traversal (like a linked list).

Output:

> python fast_bwt3.py GACGAAGGGATTCTGTGGATACTTAAACTACC 20 21 6 1 22 7 14 24 29 2 5 16 17 26 9 31 19 13 23 28 8 18 12 4 15 25 30 27 11 3 10 32

If you compare this with what we got last time you'll see it matches.

fast_bwt.py

#text = 'BANANA' text = 'GAATTCAAGCTTGGATCCGGAAAGATCTGATC' t = text + '^' N = len(t) # L is the suffix array L = list() # i: index to walk current suffix along text # j: indexes for sorted suffix already in suffix array # k: enumeration of L, to give insertion point def doOne(i,L): for k,j in enumerate(L): n = 0 while t[i+n] == t[j+n]: n += 1 if t[i+n] < t[j+n]: L.insert(k,i) return L.append(i) for i in range(1,N): doOne(i,L) bwt = ''.join([t[j-1] for j in L]) print bwt, '\n' for k,j in enumerate(L): if k and not k % 10: print print '%3d' % j,

Burroughs-Wheeler Transform-BWT (2)

I'm exploring the Burroughs-Wheeler Transform, which has become very important for fast methods of matching short reads to large genomes. The first post is here.

I found a nice page about it here, which is to be part of a series. There is a great applet(?), anyhow a tool for visualization is on that page, that can walk you through the steps of the search. I found it very helpful for understanding what's going on. I hope my explanation is also useful.

With the output from the Python script from yesterday, we'll try using BWT to search for the query GGATC in the text:

GAATTCAAGCTTGGATCCGGAAAGATCTGATC

We'll work from the end of the query backwords. The rows we're working on can be found in the output from the script, which is listed at the end of this post.

Q1 = C.
Begin by selecting all the rows that start with C:

5 CAAGCTTGGATCCGGAAAGATCTGATC^GAATT 16 CCGGAAAGATCTGATC^GAATTCAAGCTTGGAT 17 CGGAAAGATCTGATC^GAATTCAAGCTTGGATC 26 CTGATC^GAATTCAAGCTTGGATCCGGAAAGAT 9 CTTGGATCCGGAAAGATCTGATC^GAATTCAAG 31 C^GAATTCAAGCTTGGATCCGGAAAGATCTGAT

Q2 = T
So next, restricted to this selection,
find all the positions in the BWT (last column) that contain T.
The first one is

5 CAAGCTTGGATCCGGAAAGATCTGATC^GAATT

Within the BWT as a whole, this is the first T.
There are four T's in the BWT in the range selected.
So our range is now 1-4.

Next, go to the rows starting with T,
that are of rank 1-4 with respect to T:

4 TCAAGCTTGGATCCGGAAAGATCTGATC^GAAT 15 TCCGGAAAGATCTGATC^GAATTCAAGCTTGGA 25 TCTGATC^GAATTCAAGCTTGGATCCGGAAAGA 30 TC^GAATTCAAGCTTGGATCCGGAAAGATCTGA

Q3 = A.
Find all the rows with A in the BWT for this range.
The first one is:

15 TCCGGAAAGATCTGATC^GAATTCAAGCTTGGA

Within the BWT as a whole, this is the seventh A.
There are three A's in the BWT in the range selected.
So our range is now 7-9.

Next, go to the rows starting with A,
that are of rank 7-9 with respect to A:

14 ATCCGGAAAGATCTGATC^GAATTCAAGCTTGG 24 ATCTGATC^GAATTCAAGCTTGGATCCGGAAAG 29 ATC^GAATTCAAGCTTGGATCCGGAAAGATCTG

Q4 = G.

Find all the rows with G in the BWT for this range.
That would be all of them.

The first one is:

14 ATCCGGAAAGATCTGATC^GAATTCAAGCTTGG

Within the BWT as a whole, this is the 3rd G.
So our range is now 3-5.

Finally, go to the rows starting with G,
that are of rank 3-5 with respect to G:

13 GATCCGGAAAGATCTGATC^GAATTCAAGCTTG 23 GATCTGATC^GAATTCAAGCTTGGATCCGGAAA 28 GATC^GAATTCAAGCTTGGATCCGGAAAGATCT

Q5 is G.
There is one row with a match.
Since is the last letter of the query, we're done.

Our match has index 13.

That's the position of the query GATCC in the text:

GAATTCAAGCTTGGATCCGGAAAGATCTGATC +

using zero-based indexing.

Notice that we don't use the whole table, and we don't use the suffix array (the indexes) until the end. Instead, all we use is the BWT and the sorted characters (actually, the index of the first A, C, etc.):

txt GAATTCAAGCTTGGATCCGGAAAGATCTGATC^ BWT GACGAAGGGATTCTGTG^GATACTTAAACTACC chr AAAAAAAAAACCCCCCGGGGGGGGTTTTTTTT^

One more thing. If we have a selected range in the BWT and are finding say, all the G's from index 6-10, we need to find quickly how many G's come before that position in the BWT.

Not using the whole table is important, because the BWT has to be memory-efficient. It would not be useful if we required memory proportional to the square of the genome size.

I think it might be worth it to work the whole example again, using just the BWT and the sorted characters:

[ NOTE: Blogger keeps screwing with the formatting. Please let me know if it doesn't seems correct later on. ]

BWT GACGAAGGGATTCTGTG^GATACTTAAACTACC chr AAAAAAAAAACCCCCCGGGGGGGGTTTTTTTT^

Q1 = C Q2 = T + GACGAAGGGATTCTGTG^GATACTTAAACTACC AAAAAAAAAACCCCCCGGGGGGGGTTTTTTTT^ ------

The indicated T is T1 in the BWT
R = 1-4

Go to T1-T4

Q3 = A + ++ + + + + BWT GACGAAGGGATTCTGTG^GATACTTAAACTACC chr AAAAAAAAAACCCCCCGGGGGGGGTTTTTTTT^ ----

The indicated A is A7 in the BWT
R = 7-9

Go to A7-A9

Q4 = G + + + BWT GACGAAGGGATTCTGTG^GATACTTAAACTACC chr AAAAAAAAAACCCCCCGGGGGGGGTTTTTTTT^ ---

The indicated G is G3 in the BWT
R = 3-5

Go to G3-G5

Q5 = G + + +++ + + + BWT GACGAAGGGATTCTGTG^GATACTTAAACTACC chr AAAAAAAAAACCCCCCGGGGGGGGTTTTTTTT^ ---

The indicated G is G8 in the BWT
We're done.

The index (see below) is 13 for the match.

output from bwt.py

seq GAATTCAAGCTTGGATCCGGAAAGATCTGATC^ 20 AAAGATCTGATC^GAATTCAAGCTTGGATCCGG 21 AAGATCTGATC^GAATTCAAGCTTGGATCCGGA 6 AAGCTTGGATCCGGAAAGATCTGATC^GAATTC 1 AATTCAAGCTTGGATCCGGAAAGATCTGATC^G 22 AGATCTGATC^GAATTCAAGCTTGGATCCGGAA 7 AGCTTGGATCCGGAAAGATCTGATC^GAATTCA 14 ATCCGGAAAGATCTGATC^GAATTCAAGCTTGG 24 ATCTGATC^GAATTCAAGCTTGGATCCGGAAAG 29 ATC^GAATTCAAGCTTGGATCCGGAAAGATCTG 2 ATTCAAGCTTGGATCCGGAAAGATCTGATC^GA 5 CAAGCTTGGATCCGGAAAGATCTGATC^GAATT 16 CCGGAAAGATCTGATC^GAATTCAAGCTTGGAT 17 CGGAAAGATCTGATC^GAATTCAAGCTTGGATC 26 CTGATC^GAATTCAAGCTTGGATCCGGAAAGAT 9 CTTGGATCCGGAAAGATCTGATC^GAATTCAAG 31 C^GAATTCAAGCTTGGATCCGGAAAGATCTGAT 19 GAAAGATCTGATC^GAATTCAAGCTTGGATCCG 0 GAATTCAAGCTTGGATCCGGAAAGATCTGATC^ 13 GATCCGGAAAGATCTGATC^GAATTCAAGCTTG 23 GATCTGATC^GAATTCAAGCTTGGATCCGGAAA 28 GATC^GAATTCAAGCTTGGATCCGGAAAGATCT 8 GCTTGGATCCGGAAAGATCTGATC^GAATTCAA 18 GGAAAGATCTGATC^GAATTCAAGCTTGGATCC 12 GGATCCGGAAAGATCTGATC^GAATTCAAGCTT 4 TCAAGCTTGGATCCGGAAAGATCTGATC^GAAT 15 TCCGGAAAGATCTGATC^GAATTCAAGCTTGGA 25 TCTGATC^GAATTCAAGCTTGGATCCGGAAAGA 30 TC^GAATTCAAGCTTGGATCCGGAAAGATCTGA 27 TGATC^GAATTCAAGCTTGGATCCGGAAAGATC 11 TGGATCCGGAAAGATCTGATC^GAATTCAAGCT 3 TTCAAGCTTGGATCCGGAAAGATCTGATC^GAA 10 TTGGATCCGGAAAGATCTGATC^GAATTCAAGC 32 ^GAATTCAAGCTTGGATCCGGAAAGATCTGATC

Burroughs-Wheeler Transform-BWT (1)

Some time ago I was introduced to bowtie (see this post).

As the front page says:

Bowtie is an ultrafast, memory-efficient short read aligner. It aligns short DNA sequences (reads) to the human genome at a rate of over 25 million 35-bp reads per hour. Bowtie indexes the genome with a Burrows-Wheeler index to keep its memory footprint small: typically about 2.2 GB for the human genome (2.9 GB for paired-end).

So the next question is: what is the Burroughs-Wheeler Transform (BWT) and how does it speed up this process? It turns out that the discovery of BWT was not all that long ago (1994 according to wikipedia), and its application to this bioinformatics problem is more recent than that. BWT has been applied to the problem of data compression in the well known bzip tool. We'll talk about the compression part today, and deal with search in another post.

I wrote a simple Python script to explore the BWT; the listing is at the end of the post. The basic transformation couldn't be easier. Take a short text, like 'BANANA', and add a marker to the end, giving 'BANANA^'.

[UPDATE: I chose '^' because it sorts after the uppercase letters by the default sort. This may grate a little for regex pros who expect '^' to mark the beginning of something. Just define your own cmp function if you are one of these. ]

Generate all the rotational permutations of this text:

BANANA^
^BANANA
A^BANAN
NA^BANA
ANA^BAN
NANA^BA
ANANA^B

Now sort the list:

ANANA^B
ANA^BAN
A^BANAN
BANANA^
NANA^BA
NA^BANA
^BANANA

That's it! The last column, 'BNN^AAA' is the BWT of 'BANANA'. I think you can see that the runs of identical symbols would allow opportunity for efficient compression.

But here's the main point: by itself, the last column contains all of the information that was in the original text. That is, the transformation is reversible or invertible. The inversion depends on the cyclic rotation.

Here are the steps. Take the BWT and sort it:

A
A
A
B
N
N
^

By the way, that's what we already have in column 1 (or 0, for you Pythonistas).

Because of the rotation, the BWT plus column 1 gives all pairs of symbols in the original text. The symbol pairs are (BWT first, then column 1):

BA
NA
NA
^B
AN
AN
A^

Now iteratively do this: sort whatever we generated in the previous step, and then tack on the BWT in front. Here are the triplets:

BAN
NAN
NA^
^BA
ANA
ANA
A^B

Rinse, lather and repeat:

BANA
NANA
NA^B
^BAN
ANAN
ANA^
A^BA

BANAN
NANA^
NA^BA
^BANA
ANANA
ANA^B
A^BAN

BANANA
NANA^B
NA^BAN
^BANAN
ANANA^
ANA^BA
A^BANA

Although 'BANANA' is a nice simple example to begin with, in some respects it is too simple. For example:

> python bwt.py ACGTGAATTCGAAACCGGAA 11 AAACCGGAA^ACGTGAATTCG 12 AACCGGAA^ACGTGAATTCGA 5 AATTCGAAACCGGAA^ACGTG 18 AA^ACGTGAATTCGAAACCGG 13 ACCGGAA^ACGTGAATTCGAA 0 ACGTGAATTCGAAACCGGAA^ 6 ATTCGAAACCGGAA^ACGTGA 19 A^ACGTGAATTCGAAACCGGA 14 CCGGAA^ACGTGAATTCGAAA 9 CGAAACCGGAA^ACGTGAATT 15 CGGAA^ACGTGAATTCGAAAC 1 CGTGAATTCGAAACCGGAA^A 10 GAAACCGGAA^ACGTGAATTC 4 GAATTCGAAACCGGAA^ACGT 17 GAA^ACGTGAATTCGAAACCG 16 GGAA^ACGTGAATTCGAAACC 2 GTGAATTCGAAACCGGAA^AC 8 TCGAAACCGGAA^ACGTGAAT 3 TGAATTCGAAACCGGAA^ACG 7 TTCGAAACCGGAA^ACGTGAA 20 ^ACGTGAATTCGAAACCGGAA

The BWT (last column) does not group all identical letters together. And in the inversion step, the reconstructed text is not necessarily in the first row.

In the output above, I've also listed the "suffix array", which is the index into the original text where we start each row. The suffix array will become important for the application to search, next time.

bwt.py

import sys from collections import deque try: text = sys.argv[1] except IndexError: text = 'BANANA' text += '^' n = len(text) #--------------------------------------------- # BWT rL = list() dq = deque(text) for i in range(n): rL.append(''.join(list(dq))) dq.rotate(1) rL.sort() bwt = [s[-1] for s in rL] sfx = [(n - e.index('^') - 1) for e in rL] for s,e in zip(sfx,rL): print '%2d' % s, print e #--------------------------------------------- # reverse # take the last col and sort to generate col 0 c0 = sorted(bwt) print '\n'.join(c0) # get pairs as a demo: c1 = [bwt[i]+c0[i] for i in range(n)] c1.sort() print '\n'.join(c1) # iterate starting with bwt again: L = sorted(bwt) for j in range(n-2): L.sort() L = [bwt[i]+L[i] for i in range(n)] print '\n'.join(L) + '\n' L = [e for e in L if not '^' in e] assert ''.join(L[0]) == text[:-1]

Julia

Julia is a new "high-level, high-performance dynamic programming language for technical computing." Website here and rationale here. Interview with Stefan Karpinski.

I followed the instructions here to download and build it. The only issue is that I failed to read the platform-specific notes before starting, so I had the wrong (old) gfortran. There are links in the docs, but I followed the notes in this thread and got it from (download link).

[ UPDATE: re the first comment, some discussion here. ]

DNA binding sites 7

Continuing in the same vein as a number of recent posts (some links here), the Dropbox link below leads to my version of a site analysis program written in C. It's not very user friendly, for example, paths to the input files are hard-coded. A substantial part of the effort required went into ScoreList.c, which implements a function for constructing a list of scores from a list of nucleotide counts for each position, and a second function which reads such a score list in from disk.

The logic of the program was explained in the previous post. I won't spend much time explaining the details here.

Let's just say that if this blog has a motto it is: "learn by doing."

So... get in there and try writing your own version. If you get stuck, you can see how (or if) I handled the problem.

A few other notes: I tested the program by using a list of counts for crp sites obtained as described (here). The output shows the position, score and sequence of crp sites in the E. coli genome. The first few (threshold = 12) are:

$ ./test 18968 14.33 tattgtgaactatcgcaaagaa 42067 19.16 ttctgtgattggtatcacattt 42187 12.44 attggtgatccataaaacaata 49633 14.48 aagagtgacgtaaatcacactt 70157 15.02 aagtgtgacgccgtgcaaataa 141283 18.95 atgtgtgatcgtcatcacaatt 141663 12.43 tgatgtgaaaatcctcaaagat 184120 13.67 aattgtgcttattttagcattt 223495 13.60 ggatgtgaatcacttcacacaa 243785 12.34 atttctgacgttagtcatattt 268498 13.40 taatgtgaacatgatcaacgaa 281362 18.40 atatgtgatccagcttaaattt 304272 12.45 tgttgttattcactacacgttt 312539 17.25 ttttttgacatgtatcacaaat 365617 14.65 atgagtgagctaactcacatta ..

After implementing the algorithm (which handles the sequence one character at a time) I realized that it would be good to have the site sequence available at the time of printing (as shown above). This required remembering the current site's sequence, which I grafted onto the program at the very end. Luckily it didn't change the running time by much. I think now that while the circular linked list seemed like an elegant approach, it caused as many problems as it solved, and I probably wouldn't do it again.

You can tell just from looking at the patterns that we are getting good matches to the crp consensus [T/A]₃TGTGAN₆TCACA[T/A]₃. Also, the last site is the lac 1 (lacZ) site---reversed. So I think it's working OK.

Finally, it runs in less than 2 seconds! That's a huge improvement on the previous time. Zipped files on Dropbox (here).

DNA binding sites 6

A few weeks ago I had some posts about finding binding sites in DNA using a simple PSSM (position specific scoring matrix; last post here). The time required to evaluate all the potential sites in a bacterial genome of ≈ 5 million bp is prohibitively long (> 1 min on my machines), so I wanted to explore ways to do it faster. I started looking at Cython (here), but then realized that I need to brush up my C skills first (here).

I also had an idea that I thought at first would make the code faster, though probably it doesn't. Then I thought it would make the code clearer, though I realized after actual implementation that it doesn't do that either. What it does do is make our pass through the sequence simpler in logical terms, and it's the basis of my C-program to evaluate sites or motifs in a DNA sequence. We construct a circularly linked list, as illustrated schematically in the graphic. Each item (node) in the list holds the current value for an accumulating score that will become the score for an individual site in the sequence.



Probably it's best to illustrate with an example. Let's say we have a scoring system used to evaluate potential sites. For example, a site with a in position 1 (we'll use 1-based indexing here) contributes score a1, c in position 2 adds c2, etc.

Suppose we're looking for sites of length N = 4, and we have this sequence: acgtg. We construct a list like the following (shown here in rows to make the layout simple). After the initial priming, we end up with this:

-> 1 2 g1 3 c1 g2 4 a1 c2 g3

I set up a toy example with scores of
0.01 .. 0.04 for a1 .. a4; 0.11 to 0.14 for c1 .. c4, etc.
Output from the priming phase looks like this:

nt = a item = 4 before: 0.00, add: 0.01 after: 0.01 nt = c item = 3 before: 0.00, add: 0.11 after: 0.11 item = 4 before: 0.01, add: 0.12 after: 0.13 nt = g item = 2 before: 0.00, add: 0.21 after: 0.21 item = 3 before: 0.11, add: 0.22 after: 0.33 item = 4 before: 0.13, add: 0.23 after: 0.36 1 0.000000 2 0.210000 3 0.330000 4 0.360000 end priming

Now we consider the next nucleotide in the sequence: t. We add the appropriate scores for t to each item in the list, then read the score for the item that contains a total of N scores, and finally, zero that item. The score we read is the score for the site that has sequence acgt. Schematically:

t 1 t1 2 g1 t2 3 c1 g2 t3 -> 4 a1 c2 g3 t4 1 t1 2 g1 t2 3 c1 g2 t3 -> 4

Output looks like this:

t item = 1 before: 0.00 add: 0.31 after: 0.31 item = 2 before: 0.21 add: 0.32 after: 0.53 item = 3 before: 0.33 add: 0.33 after: 0.66 item = 4 before: 0.36 add: 0.34 after: 0.70 zeroing item = 4 final score = 0.70

The next nucleotide is g:

g 1 t1 g2 2 g1 t2 g3 -> 3 c1 g2 t3 g4 4 g1 1 t1 g2 2 g1 t2 g3 -> 3 4 g1

g item = 4 before: 0.00 add: 0.21 after: 0.21 item = 1 before: 0.31 add: 0.22 after: 0.53 item = 2 before: 0.53 add: 0.23 after: 0.76 item = 3 before: 0.66 add: 0.24 after: 0.90 zeroing item = 3 final score = 0.90

The site acgt has score:
a1 + c2 + g3 + t4 = 0.01 + 0.12 + 0.23 + 0.34 = 0.70


The site cgtg has score:
c1 + g2 + t3 + g4 = 0.11 + 0.22 + 0.33 + 0.24 = 0.90

Zipped files on Dropbox (here). These will change in the next few days as I debug and test the project more...

The C code for the linked list looks like this:

item *make_circular_linked_list(int N) { int i; item *current, *first, *previous; for (i=0; i<N; i++) { current = (item *) malloc(sizeof(item)); current->num = i+1; current->score = 0; if (i==0) { first = current; } else { previous->next = current; } previous = current; } current->next = first; return first; }

DNA binding sites 5

Here is an example that shows how difficult the problem of searching for sites is, at least by the method we've used. I include two scripts in the zipped files: search.EC.py and plot.py. The site we want to search for is hard-coded in the file as crp.

As its name implies, the first script searches the E. coli sequence. You'll need to get the sequence first (and save it in the right place---see the script) before this will run. We slide a window over the length of the sequence, shifting one base at a time, and score what's showing in the window under the site's scoring scheme. Scores were multiplied by 10 (to convert to ints) and then saved to disk (22 MB or so).

In the second phase, the search.EC.py script also randomizes the sequence and re-runs the same search. The plot.py script filters the results for those above a threshhold and plots a histogram of the results (see the above graphic). The red bars are the results obtained with the authentic sequence and the yellow bars are with the randomized sequence.

The point is that although the extreme high values are clearly higher in the real sequence, for a scoring range like 110 - 120 (that is 11.0 - 12.0 in the original scheme), the ratio of the likelihoods for the two models (real E. coli v. random) is not even greater than 2. So, upon observing a site with a value of 11.5, say, its significance isn't clear.

One idea that would improve the significance is to note that the genome is subject to selection. For the Crp system to work properly, there has likely been selection against randomly placed sites, so the random model is not really the appropriate one to test against.

Again, zipped files here.

DNA binding sites 4

I grabbed the data for crp and purR sites in E. coli from George Church's server (here). This is the first part of the crp data set:

>aldB -18->4 attcgtgatagctgtcgtaaag >ansB 103->125 ttttgttacctgcctctaactt >araB1 109->131 aagtgtgacgccgtgcaaataa >araB2 147->169 tgccgtgattatagacactttt

The graphics above are representations of the information analysis for crp and purR produced by site_score.py. Notice the different patterns. The important bases for crp are two short pentamers separated by one turn of the helix. Not so for purR, which suggests a different mode of binding.

According to (Schumacher 1994 PMID 7973627):

The DNA-binding domain contains a helix-turn-helix motif that makes base-specific contacts in the major groove of the DNA. Base contacts are also made by residues of symmetry-related alpha helices, the "hinge" helices, which bind deeply in the minor groove. Critical to hinge helix-minor groove binding is the intercalation of the side chains of Leu54 and its symmetry-related mate, Leu54', into the central CpG-base pair step. These residues thereby act as "leucine levers" to pry open the minor groove and kink the purF operator by 45 degrees.

It's that minor groove interaction that is giving the strong signal in the "middle" of the site.

Here is what we calculate for site scores for crp. Notice that the lac site is a relatively poor one:

$ python site_utils.py crp tnaL 20.2 nupG2 18.7 lac 17.6 cdd 17.3 deoP2 16.7 malT 16.5 .. cya 12.8 .. crp 12.0 .. lac 8.8 .. avg for 100000 random seqs: -15.64 13.16 12.18 11.87 11.51 10.78 10.66 10.58 10.56 10.2 9.87 158 sites in random seq for cutoff = 5

DNA binding sites 3

Continuing with binding site analysis, first two posts here and here.

Tom Schneider also invented Sequence Logos, which display the information for binding sites in an intuitive, graphical way (Schneider 1990 PMID 2172928). In that paper, (following Shannon) they define an uncertainty measure for each position in an alignment:



where H(l) equals minus the sum over the four nucleotides of the frequency of each base b at that position times the log₂(freq). Then, the information is:



where e(n) is a small sample correction factor. Thus, uncertainty plus information is constant, and approximately equal to 2 (bits). It's no coincidence that lacking any information about which nucleotide is present at some position in a sequence, you need to ask me two yes-no questions to obtain the identity. For example: is it a purine? Yes. Then is is adenine? Yes. Two questions, two bits.

The script site_score.py does this calculation for the fis sites example, and we plot our home-grown version of the logo as the graphic below.


You can compare that output to what is in the paper:



The colors are switched for the central position because the values for A and T are exactly equal, and we sorted to plot T on top, while Schneider did the reverse.

There is also a site on the web for making logos. To use that, we need to strip the names out of the sequence file.

FH = open('fis.sites.txt','r') data = FH.read() FH.close() for line in data.strip().split('\n')[1:]: print line.split()[1]

The only significant difference is at the middle position (11). We didn't use the reversed sequences, so we see mainly A at that position. This is an artifact of the web site's approach.



Zipped project files on Dropbox (here).

DNA binding sites 2

The first few sequences in fis.sites.txt are:

#format Schneider 1 tttgccgattatttacgcaaa 2 agtgactaaaatttacactca 3 gtggtgcgataattactcata 4 attgcatttaaaatgagcgtg 5 attggtcaaagtttggccttt

The first few lines of output from site_utils.py:

ttatgtacaaatagtaagaaatgtctgaga.. [45, 10, 26, 39] 0.567 -1.603 -0.2245 0.3605 ttggtatatatactatacacctatatttga.. [28, 22, 21, 49] -0.1175 -0.4655 -0.5326 0.6898

For Schneider's approach each site is present in the alignment in both forward and reverse complement orientation. In the output we are looking at each column of the alignment. The first base 't' of the top line of output is from sequence #1, the second base 't' is from the reverse of sequence #1, the third base 'a' is from sequence #2, and so on. The counts are given next, and sum to N = 120 = 60 * 2. The scores are calculated as given previously:

2 + log2(freq) - 0.018

The scores for the first 11 sequences are:

1 12.2 2 11.8 3 9.0 4 6.5 5 12.2 6 8.5 7 8.4 8 4.6 9 12.0 10 5.3 11 10.4

If you look at the graphic from last time (or the paper) you'll see that we match. So I think we're doing things correctly.

The last thing we do when running this basic site_utils.py file as __main__ is to look at some random sequences. For starters, we calculate

avg for 100000 random seqs: -10.85

which is roughly -0.5 for each position in the 21 nucleotide alignment. According to my understanding the purpose of the correction term was to make this be zero, and I'm not sure why it isn't.

The top 10 sites from the random sequences had scores:

11.72 10.73 10.73 10.52 10.48 10.28 10.08 10.08 10.07 9.86

There were a total of 9 sequences with scores > 10, while only 14 of the 60 authentic sites scored that high. So depending on where our cutoff is for an authentic site we'll get a lot of false positives. For this example, with a cutoff of 5, I found 267 of them.

DNA binding sites 1

I haven't written much about classification of DNA binding sites or motifs on the blog. (I did a series on Gibbs sampling to find new motifs starting here).

But a couple of years ago, before the blog started, I did some work on this and there is a modest writeup on my pages at mac.com (here). At that time, I put a lot of effort into providing a Quartz GUI for the basic functionality, but I haven't tested it recently and I'd be pretty surprised if it still works. Still, that discussion (and see here) is reasonable. Let's see if we can improve upon it.

The reason the subject came up is in considering Cython (here). One use case that's definitely appealing is to predict new members of the family in a bacterial genome (given a set of known binding sites), where there are (just about) as many sites to be tested as base pairs. That can take a substantial time using Python (a few minutes).

To start working on this problem, we need some data. A repository of alignments for various transcription factors of E. coli is available from George Church's lab (here), but I'm going to use Tom Schneider's method (website) for constructing the scoring system, so I'll use the example from his paper about fis (Hengen 1997 PMID 9396807). Here is a graphic from that paper showing a few of the 60 sites they analyzed and the bit scores for each site. We're going try to to recreate this. I've got a file with the 60 sequences; you can get them and the code we'll be using from Dropbox (here).



As usual for these problems, we consider each position in the alignment to be independent, and add up the scores for all the columns to obtain a total score. The score for a single column of the alignment depends on the the count of nucleotides at that position. Suppose we have N = 20 sequences and this distribution:

A:5 C:5 G:5 T:5

We calculate the frequency as 5/20 = 0.25 for each base, and then the score is 2 + log₂(freq) - correction. The correction is for the small sample size of known sites (see here to begin) and it's about 0.018 in this particular case, but we'll neglect it for this brief discussion. The method essentially computes a log odds score for the competing hypotheses of a real binding site versus random sequence.

The score for each base in this alignment, which seems to be random, is calculated as 2 log₂(0.25) = 0. On the other hand, suppose we have:

A:18 C:1 G:1 T:0

Then the scores are:

A = 2 + log2(0.9) = 2 - 0.152 = 1.848 C = G = 2 + log2(0.05) = 2 - 4.322 = -2.322

We assign a pseudocount of 1 for T, even though no site had T, so the score for T is the same as for C and G. In this way we end up with an m x n matrix of floats, where m = 4 (ACGT) and n = the length of the alignment.

To score a candidate site, we observe each nucleotide in turn, and retrieve the corresponding score for that nucleotide and position in the alignment. For example, in the non-random case if the site contains A at that position we add 1.848 to the score, while if it has C we subtract 2.322 from the score. Here we see the great strength of this approach: we don't just reward sequences for being close to the "consensus", we penalize them for having a nucleotide that is rarely observed at the corresponding position in authentic sites.

That's the simple, basic idea. As I say, the files will be on Dropbox (here). Next time we'll show some results.

Likelihood revisited

As I've said many times, I like Higgs & Attwood a lot. But I ran into one section of the book that's giving me some trouble, and I'd like to work through it here. It's about likelihood, on pp. 229-230.

Let's go back to the beginning, and start by making the standard application of Bayes rule. The joint probability of two events A and B is equal to:

P(A,B) = P(A|B) P(B) = P(B|A) P(A)

Now, what we really want to talk about is not just any two events A and B, but data and models. In the figure, we sketch out an event space of all possible outcomes before seeing any data, which is divided into rectangles of different sizes (10 horizontal steps and 20 vertical steps). Each horizontal step is the same size (1/10), as is each vertical step (1/20).



Let's say we partition the event space using horizontal steps and assign each partition to a model: A through J.

Consider the purple rectangle at the extreme lower right. In this part of the event space, data z was observed and model J is correct.

The size of this rectangle---the probability P(z,J)---is one-half the size of the purple rectangle to its left, labeled 1 (for its size, which is 1% of the total). The total area of the extreme lower right rectangle: P(z,J) is 1/20 x 1/10 = 1/200 = 0.005 (i.e. 0.5%). In symbols, this calculation is:

P(z,J) = P(z|J) P(J) = 1/20 x 1/10

The probability of data z, given that model J is correct, times the probability that model J is correct, is equal to the probability of both events: that we will observe the data z and that model J is correct.

Now, suppose we back up and think about all the models. We observe data z and ask, what is P(J|z)? That is, rather than asking how often model J will generate data z, we know the data already, and we want to learn the probability that model J is correct.

Using Bayes rule and rearranging terms we just do this:

P(J|z) = P(z,J) / P(z) = P(z|J) P(J) / P(z)

where

P(z) = (5 + 4.5 + 4 + 3.5 + 3 + 2.5 + 2 + 1.5 + 1 + 0.5)/100 = 0.275 P(J|z) = 0.005 / 0.275 ≈ 0.2

Having observed data z, we can compute P(J|z). The three terms we need to do the calculation are the:

prior: P(J) likelihood: P(z|J) evidence: P(z) = P(z|A) + P(z|B) + .. P(z|J)

What's nice about this figure is that it's easy to see the influence of our prior evaluation of P(J). In this case, we've taken an agnostic (uninformative) prior and assumed that each model is equally probable (the horizontal width of each rectangle is the same). We could accomodate prior information by adjusting these widths.

For any well-specified model we should be able to calculate P(data|model). But by the time we're doing the calculation, we've already seen the data. So then it seems contrived to talk about the probabiity. A new term is used, the likelihood, and we talk about the likelihood L of model J given what we just observed (data z):

L = P(z|J)

But what Higgs and Attwood say is different:

.. the likelihood of the data according to each of the models. This is written L(D|M_k) and is read as "the likelihood of the data, given model k." These likelihoods are defined as the probability distributions over all possible sets of data that could be described by the model. If we sum the likelihoods over all possible data for any given model, we must get one (i.e. the likelihood is normalized).

They are saying that something must have happened, and in the event space where model J is true, given that normalization (only J is true) then

P(z) + P(y) + P(x) = 1.

They are summing vertically inside the box that is model J.



In contrast, Felsenstein uses the phrase "Likelihood Prob(D|H)", and gives the example of estimating the value of parameter p for a binomial distribution, given data of 5 heads and 6 tails. Given that the model parameter p is any particular value, we can calculate P(data|p). The maximum likelihood estimate of p works out to 5/11.



Felsenstein says:

Note that although this looks rather like a distribution, it is not. It plots the probabilities of the same data D for different values of p. Thus it does not show the probabilities of different mutually exclusive outcomes, and the area under the curve need not be 1.

The equivalent observation in our case is that the bar-graph of purple rectangles in the figure shows the probabilities of the same data z for different models. These are not mutually exclusive outcomes (we could have obtained other data), and the total area of the bars need not be 1.

I'm with Felsenstein. Talk about the "likelihood of the model having observed some data", and don't assume that the likelihoods of different models given the same data should sum to 1. They don't.

Gene Ontology: showing the graph 2

A brief final update on the GO project. I modified GOUtils.py to remember both the child and parent in moving up the tree: handle_request(D,target,debug=False) returns a list of pairs which are of the form

('GO:0006139', 'GO:0034641').

Having the GO data in memory, it's easy to find the corresponding names.

Now what we need is some graphing code that can obtain the name from the go_id and make the plot. Without worrying about it too much yet, I found some old code that works, and used that. I remember that it was hard to parse everything to make the labels, so I wasn't eager to revisit the logic behind it yet. Here is the graph for a transcription factor, I think it was GAL4 but I've lost track. Zipped project files here. Hope this helps someone.

Gene Ontology: showing the graph 1

The graphing program I used is Graphviz, from here. On the laptop I'm using, I have Graphviz.app (2.20), which was installed about 2 years ago, but they are currently up to version 2.26, and 2.27 is in development. The Mac OS X .app version seems to not be upgraded anymore, but my copy still works. You can also get Graphviz from MacPorts.

I can launch the command line version like so (or whatever the path to the MacPorts version is):
/Applications/Graphviz.app/Contents/MacOS/Graphviz ~/Desktop/example1.dot

but I normally double-click (or just drop a file with commands on top of the application icon). A very simple example gives the graphic above:

digraph G { A -> B; A -> C; B -> C; C -> D; D[color="cyan", style=filled] }

Graphviz does the heavy lifting for us. I like that!

Code for the example from last time:

digraph G { node[style="filled", fontsize="22",color="cyan"]; A[label="nuclear\npart"]; B[label="intracellular\norganelle\npart"]; A -> B; C[label="protein\ncomplex"]; D[label="macromolecular\ncomplex"]; C -> D; E[label="intracellular\npart"]; F[label="cell\npart"]; E -> F; G[label="organelle\npart"]; B -> G; B -> E; H[label="cellular_component"]; F -> H; I[label="mediator\ncomplex"]; I -> C; J[label="nucleoplasm\npart"]; I -> J; J -> A; D -> H; G -> H; }

Gene Ontology continued

This is the continuation of a project using the Gene Ontology (first post here). For this part, you'll need to get the annotations associated with the yeast genome---at least that's what I used (here). In the project files (link below) you'll find a short script that loads the data from this file. It expects to find the file in the db folder.

Another short script useGO.py just exercises things a bit. We load the GO data and the yeast annotations. Given a target list (in this case ['pheromone']), then we look for all the yeast genes containing that word in the description field (at index 9 of the original yeast db file). We recover these GO ids and print all the applicable GO terms, obtained using the recursive code from the first post.

Sample output shows a single one of the genes found:

MFA1 Mating pheromone a-factor, made by a cells ['GO:0000750'] ['pheromone-dependent signal transduction involved in conjugation with cellular fusion'] GO:0000750 pheromone-dependent signal transduction involved in conjugation with cellular fusion GO:0007186 G-protein coupled receptor protein signaling pathway GO:0007166 cell surface receptor linked signaling pathway GO:0023033 signaling pathway GO:0023052 signaling GO:0008150 biological_process GO:0032005 regulation of conjugation with cellular fusion by signal transduction GO:0007165 signal transduction GO:0050794 regulation of cellular process GO:0050789 regulation of biological process GO:0065007 biological regulation GO:0031137 regulation of conjugation with cellular fusion GO:0046999 regulation of conjugation GO:0043900 regulation of multi-organism process GO:0048610 cellular process involved in reproduction GO:0009987 cellular process GO:0022414 reproductive process

I think you can see what GO is supposed to be about. We gradually progress to more and more general categories as we work our way up the tree.

What's not obvious in the approach I used so far is that these chains of terms end with one of three different major categories. These are:

GO:0003674 ['molecular_function'] GO:0008150 ['biological_process'] GO:0005575 ['cellular_component']

The other thing is that I've obscured the branching, but I have a modification to the code that gives this information. And I have a graph that plots it. More in a later post. Zipped project files here.

GO (Gene Ontology)

According to the main page:

The project provides a controlled vocabulary of terms for describing gene product characteristics and gene product annotation data from GO Consortium members, as well as tools to access and process this data.

The idea is for people to browse from the internet, but you can download the ontology relationships (download link). The database consists of entries separated by double newlines which look like:

[Term] id: GO:0008883 name: glutamyl-tRNA reductase activity namespace: molecular_function def: "Catalysis of the reaction: (S)-4-amino-5-oxopentanoate + NADP(+) + tRNA(Glu) = L-glutamyl-tRNA(Glu) + H(+) + NADPH." [EC:1.2.1.70, RHEA:12347] subset: gosubset_prok synonym: "L-glutamate-semialdehyde: NADP+ oxidoreductase (L-glutamyl-tRNAGlu-forming)" EXACT [EC:1.2.1.70] xref: EC:1.2.1.70 xref: KEGG:R04109 xref: MetaCyc:GLUTRNAREDUCT-RXN xref: RHEA:12347 is_a: GO:0016620 ! oxidoreductase activity, acting on the aldehyde or oxo group of donors, NAD or NADP as acceptor

and so on... There's a single header (which doesn't have [ ] in the first line), and a few entries at the end that start with [Typedef]. A format guide is here. (I didn't read it).

What we want to do is follow all the 'is_a' links. When I looked at GO previously, I found some cyclical references. So I put in a test for whether a particular item has been seen before, not unlike the tree-traversal code here. Anyway, the following is the output for this particular target.

We print the details for this item, and then we follow a chain of 'is_a' all the way up to GO:0003674 molecular_function.

GO:0008883 def "Catalysis of the reaction: (S)-4-amino-5-oxopenta .. id GO:0008883 is_a GO:0016620 ! oxidoreductase activity, acting on th .. name glutamyl-tRNA reductase activity namespace molecular_function subset gosubset_prok synonym "L-glutamate-semialdehyde: NADP+ oxidoreductase (L .. xref EC:1.2.1.70 KEGG:R04109 MetaCyc:GLUTRNAREDUCT-RXN RHEA:12347 GO:0008883 glutamyl-tRNA reductase activity GO:0016620 oxidoreductase activity, acting on the aldehyde or oxo group of donors, NAD or NADP as acceptor GO:0016903 oxidoreductase activity, acting on the aldehyde or oxo group of donors GO:0016491 oxidoreductase activity GO:0003824 catalytic activity GO:0003674 molecular_function

I'll be doing more with this, so I put the zipped project files up on Dropbox (here). I will update the zip at that link as I work more on the project. You will need to add the GO database to the /db folder for it to run. It's about 20 MB.

Code to find restriction sites

Having downloaded a list of restriction enzymes from REBASE earlier today for another post, I couldn't resist writing a short program to find restriction sites in a DNA sequence. The input enzyme data looks like this:

AarI CACCTGC (4/8) AatII GACGT^C AbsI CC^TCGAGG AccI GT^MKAC AceIII CAGCTC (7/11) AciI CCGC (-3/-1)

We ignore the numbers (4/8) etc. To make it more interesting, I used a simple form of Python's regular expressions (docs here), which are pre-compiled.

The results can be limited to non-degenerate six-cutters (or longer) and sorted by the index of the match. It took a bit more than an hour, and was great fun. Naturally, the code hasn't been tested carefully. (Note: 0-based indexing).

output:

$ python script.py Eco47III AGCGCT 11 AGCGCT HaeII RGCGCY 11 AGCGCT TsoI TARCCA 23 TAACCA HgiCI GGYRCC 35 GGCACC AflIII ACRYGT 56 ACGCGT MluI ACGCGT 56 ACGCGT AclI AACGTT 62 AACGTT BclI TGATCA 83 TGATCA BspGI CTGGAC 93 CTGGAC MstI TGCGCA 107 TGCGCA BsgI GTGCAG 120 GTGCAG HindII GTYRAC 140 GTCAAC BspMI ACCTGC 190 ACCTGC

code listing:

import re sample_dna = ''' ATGACCCTTTTAGCGCTCGGTATTAACCATAAAACGGCACCTGTATCGCT GCGAGAACGCGTAACGTTTTCGCCGGACACGCTTGATCAGGCGCTGGACA GCCTGCTTGCGCAGCCAATGGTGCAGGGCGGGGTCGTGCTGTCAACCTGT AACCGTACAGAGCTGTATCTGAGCGTGGAAGAGCAGGATAACCTGCAAGA''' # load data from downloaded file # http://rebase.neb.com/rebase/link_proto def load_data(fn = 'link_proto.txt'): FH = open(fn,'r') data = FH.read() FH.close() L = data.strip().split('\n\n') return L # parse data into names and sites def preprocess(data): L = data.strip().split('\n') names = list() sites = list() for e in L: n,s = e.split(' ',1) names.append(n) words = s.strip().split() if words[0][0] == '(': sites.append(words[1]) else: sites.append(words[0]) return names,sites # codes for degenerate positions # http://www.bioinformatics.org/sms/iupac.html def get_pattern(s): D = { 'A':'A','C':'C','G':'G','T':'T', 'R':'[AG]','Y':'[CT]','N':'.', 'S':'[GC]','W':'[AT]','K':'[GT]', 'M':'[AC]','B':'[CGT]','D':'[AGT]', 'H':'[ACT]','V':'[ACG]', '^':''} rL = [D[c] for c in s] return ''.join(rL) # dictionary of pre-compiled regexps def make_dict(data): names,sites = preprocess(data) D = dict() for n,s in zip(names,sites): p = get_pattern(s) p = re.compile(p) i = s.find('^') s = s.replace('^','') rD = { 'name':n,'site':s, 'pattern':p,'i':i } D[n] = rD return D def search(dna,D,minlength=4,ignore_ambig=True): N = max([len(n) for n in D.keys()]) M = max([len(D[n]['site']) for n in D]) rL = list() for n in D: rD = D[n] n,p,s = rD['name'], rD['pattern'],rD['site'] if len(s) < minlength: continue if ignore_ambig and 'N' in s: continue m = p.search(dna) if m: i = m.start() j = i + len(rD['site']) e = [n.ljust(N),s.ljust(M)] e += [i,dna[i:j]] rL.append(e) rL.sort() return rL def show(result): def f(s): return s[2] # index of match result = sorted(result, key=f) for line in result: # left i as int, so convert to str i = line[2] line[2] = str(i).rjust(4) print ' '.join(line) if __name__ == '__main__': dna = ''.join([c for c in sample_dna if c in 'ACGT']) data = load_data() L = data[3].strip() # type II enzymes only D = make_dict(L) result = search(dna,D,minlength=6) show(result)

Unifrac analysis: simulating sequences

I want to show an example of phylogenetic analysis using Unifrac (website (old), PMID 16893466). Before we do that, we need to get some sequences to work with. I am going to do that here, and I'll talk about Unifrac next time.

We're going to simulate sequences from two different populations. The simulation is based on 4 authentic rRNA sequences that we get from Genbank:

AY005045.1 Streptococcus_mitis_bv2 AB302401.1 Pseudomonas_cinnamophila L14639.1 Capnocytophaga_gingivalis AF411020.1 Achromobacter_xylosoxidans_AU1011

We will simulate the two samples by picking sequences with some distribution. To start with, I'm going to try this one:

A_distr = {0:1,1:3,2:5,3:1} B_distr = {0:0,1:3,2:1,3:5}

So, sample A will have 1 sequence from S. mitis, 3 from P. cinnamophila, etc. As I work through the Unifrac analysis, I might go back and change this to make the results more interesting, but this is a reasonable start.

Also, we'll mutagenize each sequence. The script to set this up is at the end of the post. Note that you must have PyCogent installed to follow along at home. You will also need muscle and FastTree installed somewhere accessible on your $PATH.

Here is a tree based on the sequences as printed out at the Terminal by the script.

/-A3 | | /-A4 | /0.042---| | | \-B2 |-0.744---| | | /-B1 ---------| \0.623---| | \-B3 | | /-A2 | | | | /-A1 | | | | | | /-A7 \0.578---| /0.966---| /0.644---| | | | | \-A6 | | | | | | \1.000---| /-A9 | | | /0.389---| | | | | \-B4 \0.996---| \0.225---| | | /-A5 | \0.822---| | \-A8 | | /-A10 | | \1.000---| /-B8 | | \0.591---| /-B5 | /0.950---| | | \-B6 \0.618---| | /-B7 \0.742---| \-B9

import random, os, sys random.seed(137) from cogent import LoadSeqs, DNA from cogent.db.ncbi import EFetch from cogent.app.muscle import align_unaligned_seqs from cogent.app.fasttree import build_tree_from_alignment def fetch_ncbi_data(ofile,s): # get the seqs from Genbank input = [e.split() for e in s.strip().split('\n')] id_list = [t[0] for t in input] names = [t[1] for t in input] ef = EFetch(id=','.join(id_list), rettype='fasta') data = ef.read().strip() # title lines are too long, replace by genus_species rL = list() for i,e in enumerate(data.split('\n\n')): old_title, seq = e.strip().split('\n',1) new_title = '>' + names[i] seq = seq[:500] rL.append('\n'.join([new_title,seq])) FH = open(ofile,'w') FH.write('\n\n'.join(rL)) FH.close() def mutagenize(seq, percent=5): L = list(seq) D = { 'A':'CGT', 'C':'AGT', 'G':'ACT', 'T':'ACG' } N = int(percent / 100.0 * len(seq)) X = len(seq) for i in range(N): j = random.choice(range(X)) nt = L[j] if not nt in 'ACGT': continue L[j] = random.choice(D[nt]) return ''.join(L) def distribute_seqs(ifile,ofile,mut_freq=10): # set up our two samples FH = open(ifile,'r') data = FH.read().strip().split('\n\n') FH.close() A_distr = {0:1,1:3,2:5,3:1} B_distr = {0:0,1:3,2:1,3:5} A = list() B = list() x = y = 0 for i,e in enumerate(data): title,seq = e.split('\n',1) seq = ''.join(seq.split()) for j in range(A_distr[i]): x += 1 copy = mutagenize(seq[:],mut_freq) new_seq = DNA.makeSequence(copy,'A' + str(x)) A.append(new_seq) for j in range(B_distr[i]): y += 1 copy = mutagenize(seq[:],mut_freq) new_seq = DNA.makeSequence(copy,'B' + str(y)) B.append(new_seq) FH = open(ofile,'w') L = [seq.toFasta() for seq in A + B] FH.write('\n\n'.join(L)) FH.close() def align_seqs(ifile,ofile): seqs = LoadSeqs(ifile, moltype=DNA, aligned=False) aln = align_unaligned_seqs(seqs, DNA) aln.writeToFile(ofile) return aln def get_tree(ifile): aln = LoadSeqs(ifile, moltype=DNA, aligned=True) tr = build_tree_from_alignment(aln,moltype=DNA) return tr #=============================================== s = ''' AY005045.1 Streptococcus_mitis_bv2 AB302401.1 Pseudomonas_cinnamophila L14639.1 Capnocytophaga_gingivalis AF411020.1 Achromobacter_xylosoxidans_AU1011 ''' fn1 = 'rRNA_gb.fasta' fn2 = 'samples.fasta' fn3 = 'samples.aln.fasta' fn4 = 'samples.tree' if not os.path.exists(fn1): fetch_ncbi_data(fn1,s) if not os.path.exists(fn2): distribute_seqs(fn1,fn2) if not os.path.exists(fn3): aln = align_seqs(fn2,fn3) tr = get_tree(fn3) print tr.asciiArt() tree_str = tr.getNewick(with_distances=True) FH = open(fn4,'w') FH.write(tree_str + '\n') FH.close()

Matplotlib in OS X: Heat Maps

I've been exploring matplotlib. It's complex enough that I probably need to get a book. I decided to revisit the problem of producing heat maps. In the last couple of posts (here and here), I mentioned a couple of methods that might be used: (scatter, pcolormesh, and pcolor).

But, today I'm going to use matplotlib.patches.Rectangle. Rather than post the code (140 lines in two scripts), I made a zipped Archive that includes:

• a directory containing three data files
(counts, topo.colors, colors for the rownames)
• Helper.py and HeatMapper.py
• example.pdf

Above is a part of the figure the code draws The whole thing is up on Dropbox here.

The code to draw the boxes is:

p = Rectangle((x,y),width=dx,height=dy, facecolor=box_color,edgecolor='w') ax.add_patch(p) if n == 0: continue t = ax.text(x+dx/2,y+dy/2,str(n), color=label_color,fontsize=12, fontname='Helvetica', ha='center',va='center')

It took significantly less time to write this than the Cocoa app that does the same thing, and I think it gains a lot in flexibility.

Here is an example in black and white. It uses the result from binCounts() as the color (i.e. gray).

Matplotlib: topo.colors

I looked around a bit on the web for topo colors for matplotlib, but didn't see anything. So, I used my new understanding of how to make a LinearSegmentedColormap, and the previous analysis of R's topo.colors from here. Reprinting the figure from that post



we see how the RGB components vary individually over the course of the range 0-1000 (or 0.000 to 1.000). Examining the values for the colors obtained from R at those points (in R do: topo.colors(1000), I estimated these values:

1 4c00ff 84 0000ff 334 00e5ff 335 00ff4d 418 00ff00 667 e6ff00 668 ffff00 850 ffdd62 1000 ffe0b3

I converted the hex values for the colors at the breakpoints to fractions of 256 in Python:

L = ['4c','b3','dd','e0','e5'] >>> for e in L: ... print e, int(e,16), int(e,16)*1.0/256 ... 4c 76 0.296875 b3 179 0.69921875 dd 221 0.86328125 e0 224 0.875 e5 229 0.89453125

And that allowed me to construct the Colormap. In filling out the map, I came to appreciate the unusual data structure. It's natural to extend something like this:

x y0 y1 x y0 x

by adding the next y1 to the middle row (the value of the color at the beginning of the interval), and the next y0 to the last row (the value at the end of the interval). Here's what I got. It looks pretty good to me:


Here's the code:

import matplotlib.pyplot as plt import matplotlib import numpy as np # topo colors r = ((0.000, 0.000, 0.297), (0.084, 0.000, 0.000), (0.418, 0.000, 0.000), (0.667, 0.895, 0.895), (0.668, 1.000, 1.000), (1.000, 1.000, 1.000)) g = ((0.000, 0.000, 0.000), (0.084, 0.000, 0.000), (0.334, 0.895, 0.895), (0.335, 1.000, 1.000), (0.667, 1.000, 1.000), (0.850, 0.863, 0.863), (1.000, 0.875, 1.000)) b = ((0.000, 0.000, 1.000), (0.334, 1.000, 1.000), (0.335, 0.297, 0.297), (0.418, 0.000, 0.000), (0.667, 0.000, 0.000), (1.000, 0.700, 1.000)) colors = {'red':r, 'green':g, 'blue':b} f = matplotlib.colors.LinearSegmentedColormap m = f('my_color_map', colors, 256) A = np.array(range(100)) A.shape = (10,10) plt.pcolor(A,cmap=m) plt.colorbar() plt.savefig('example.png')