nothing in biology makes sense except in the light of evolution

I've been reading Mike Yarus's book: Life from an RNA world (Amazon). It's a very readable account of evolution from the perspective of sequences and the RNA world. Mike's a highly intelligent guy, and that and his wit inform every page. In one chapter, he gets Rob Knight and Steve Freeland to do an evolution simulation.

we let a computer write out a random string, mutate 1 in 100 characters in each generation, and select changes only if they match Dobzhansky

In my version:

in each generation we pick one position

if it already matches Dobzhansky, continue to the next generation

mutate to a random choice from the set of symbols

This isn't a very good model of evolution. It was just fun to spend a few minutes coding it.

Here is the beginning, middle and end of one run:

> python evolve.py 0 LVSTUknwIfGIRHFgHESZ M zthWtNQTk qgtoeMvjeJAzOidKWEZO ZwNjDyCvq 100 LVSTUknwIfGIRHFgHESZ M zthetNQTk qgtoeMvjeJAzOidKWEZO ZwNjDyCvq 200 LVSTUknwIfGIRHFgHESZ M zthetNQTk qgtoeMvjeJAzOidhWEZO ZwNjDyCvq 300 LVSTiknwIfGIRHFgHESZ M zthetNQ k qgtoeMvjeJAzOidhWEZO ZwNjuyivq 6100 notTing if biHlogy makes sense except iv the Oight of Zvolution 6200 notTing if biHlogy makes sense except iv the Oight of Zvolution 6300 nothing if biHlogy makes sense except iv the Oight of Zvolution 6400 nothing if biHlogy makes sense except iv the Oight of Zvolution 13800 nothing in biHlogy makes sense except in the light of evolution 13900 nothing in biHlogy makes sense except in the light of evolution 14000 nothing in biHlogy makes sense except in the light of evolution 14100 nothing in biHlogy makes sense except in the light of evolution 14167 nothing in biology makes sense except in the light of evolution 14167 nothing in biology makes sense except in the light of evolution

import random s = 'nothing in biology makes sense' s += ' except in the light of evolution' N = len(s) symbols = 'abcdefghijklmnopqrstuvwxyz ' symbols += symbols.upper() L = [random.choice(symbols) for i in range(N)] print ' 0 ' + ''.join(L) i = 0 while L != list(s): i += 1 v = i and not i % 100 if v: print str(i).rjust(5), j = random.choice(range(N)) if L[j] == s[j]: if v: print ''.join(L) continue c = random.choice(symbols) if c == s[j]: L[j] = c if v: print ''.join(L) print print str(i).rjust(5), s print str(i).rjust(5), ''.join(L)

Mutual Information (3)

We're working on a paper from Michael Laub's lab at MIT (Skerker et al 2009 PMID 18555780). Previous posts here and here.

Now it's time to analyze the data. The self comparisons (a single position in the alignment, analyzed for mutual information against itself) yield info values ranging from 4.1 to 0.01.

90 4.102 260 4.025 122 3.993 209 3.976 234 3.944 101 3.939 235 3.938 152 3.935 63 3.91 222 3.887 159 0.191 250 0.138 132 0.133 291 0.086 199 0.067 242 0.052 16 0.012 106 0.012 156 0.012 158 0.012

As you can see in the screenshot, the residue in column 16 (1-based indexing), which has a very low score, is the conserved (catalytic) histidine.



We'll filter out the self comparisons for the plots. Here is the histogram of information values that I got. It's a bit different from the paper, but not much.



In the next part, we'll try to match up the residues which look like they might interact (that have high mutual information) and see if that makes sense in terms of the protein structures.

I picked a familiar HK/RR pair to do this: EnvZ and OmpR. This screenshot of part of Fig 2 shows a section of each protein.



Searching in the alignment file (the alignments don't have titles that I recognize), I recovered a sequence that I think is probably the right one. It matches the Figure as far as I checked:

>26250004-26250005/1-457 KQLADDRTLLMAGVSHDLRTPLTRI RLATEMMSAESINKDIEECNAIIEQ FIDYLRTGMADLNAVLGEVIA--AE SGYEREIETAL-YVKMHPLSIKRAV ANMVVNAARY-GNGWIKVSSGTEAW FQVEDDGPGIAPEQRKHLFQPFVRG DISGTGLGLAIVQRIVDNHNGMLEL GTSERGGLSIRAWLPNYKILVVDDD MRLRALLERYLTEQGFQVRSVANAE QMDRLLTRESFHLMVLDLMLPGEDG LSICRRLRSQSPMPIIMVTAKGEEV DRIVGLEIGADDYIPKPFNPRELLA RIRAVLRR

I put the newlines in to help me count. The alignment is 308 residues total. From Fig 2,

EnvZ (the HK) sequence starts with:

AGVKQLADDRTLLMAGVSHDL KQLADDRTLLMAGVSHDL

The second line above is from the alignment. The sequence doesn't start at residue 1, however. By my calculation, residue 0 in the alignment (the K) is residue 15 in the protein (1-based index). So we'll adjust the indexes we obtain by adding 15 to compare with the actual protein sequence.

OmpR (the RR) sequence starts with

MQENYKILVVDDDMRLRALLER NYKILVVDDDMRLRALLER

The second line above is from the alignment, where there's an N at position 0 in this fragment. By my calculation, that N is residue 3 of OmpR in Fig 2 (1-based index), and residue 190 of the alignment, so we subtract 187 for values >= 190.

Going back to the plotting script, we filter the data for pairs in which one comes from the HK and one from the RR, and print the top 20. We do the math as outlined above. As a check, we grab the putative EnvZ/OmpR sequence from the alignment file and print the sequence starting at the position we've identified. Here are the results for the top 20 (in each pair the RR is printed first and the HK second).

[Note: to keep the code simple, I ignored the situation where a column contains '-' a gap in the alignment for some of the sequences. That's what's giving us the two strange results below at position 16 and 20.]

> python plot.py 14 RLRAL 42 ATEMM 0.822 18 LLERY 42 ATEMM 0.816 22 YLTEQ 42 ATEMM 0.769 15 LRALL 37 TRIRL 0.730 15 LRALL 42 ATEMM 0.694 56 MLPGE 54 DIEEC 0.688 14 RLRAL 54 DIEEC 0.678 21 RYLTE 42 ATEMM 0.677 14 RLRAL 38 RIRLA 0.668 83 KGEEV 22 TLLMA 0.667 22 YLTEQ 21 RTLLM 0.663 4 YKILV 42 ATEMM 0.658 83 KGEEV 54 DIEEC 0.657 22 YLTEQ 18 ADDRT 0.650 22 YLTEQ 54 DIEEC 0.646 18 LLERY 86 --AES 0.644 15 LRALL 54 DIEEC 0.643 18 LLERY 38 RIRLA 0.636 22 YLTEQ 45 MMSAE 0.633 22 YLTEQ 86 --AES 0.631

The two residues with highest mutual information are OmpR residue 14 and EnvZ residue 42. It looks pretty good to me.

[ UPDATE: The heatmap looks pretty boring, so I'm going to skip it.
But I plotted the top 20 interactions (graphic at the top of the post), the repetition indicates we're on the right track.. ]

import sys from utils import load_data import matplotlib.pyplot as plt data = load_data('results.2.txt') data = data.strip().split('\n') data = [e.split() for e in data] data = [(int(t[0]), int(t[1]), float(t[2])) for t in data] def f(t): return float(t[2]) def part1(): L = [t for t in data if t[0] == t[1]] L = sorted(L,key=f,reverse=True) for t in L[:10]: print str(t[0]+1).rjust(3),round(t[2],3) print for t in L[-10:]: print str(t[0]+1).rjust(3),round(t[2],3) sys.exit() #part1() L = [t[2] for t in data if t[0] != t[1]] X = 1.0 plt.hist(L,bins=X*50) ax = plt.axes() ax.set_xlim(0,X) plt.savefig('example.png') # t[0] always > t[1] N = 190 data = [t for t in data if t[0] >= N and t[1] < N] aln = load_data('cell3925mmc4.fa') aln = aln.strip().split('>')[1:] aln = [e for e in aln if e.startswith('26250004-26250005/1-457')] envZ_ompR = aln[0].split('\n')[1] for t in sorted(data,key=f,reverse=True)[:20]: i = t[0] rr = i - 187 j = t[1] hk = j + 15 print str(rr).rjust(3), envZ_ompR[i:i+5], print str(hk).rjust(3), envZ_ompR[j:j+5], print '%3.3f' % t[2]

Mutual information (2)

We're working on a paper from Michael Laub's lab at MIT (Skerker et al 2009 PMID 18555780). The first post is here.

In this part we'll load the alignment (supplementary data file S4---the annotation on the page is incorrect), and crunch the numbers. I just write the results to disk.

We'll do the analysis in another post.

python info.py > results.txt

import sys from utils import load_data import info_helper as ih #aln = 'AASSASSTTT\nNMWWNTTKKS\nGTSNTYRSTA\nGGGGGGGGGG' fn = 'cell3925mmc4.fa' data = load_data(fn) data = data.strip().split('>')[1:] data = [e.split('\n')[1].strip() for e in data] def show(data): print 'starting:', len(data) for i in range(7): print i, L = [e for e in data if e.count('-') <= i] print len(L) sys.exit() #show(data) def transpose(L): R = range(len(L[0])) rL = list() for i in R: rL.append(''.join([item[i] for item in L])) return rL data = [e for e in data if e.count('-') <= 4] #data = data[:100] cols = transpose(data) pD = ih.make_prob_dict(cols) info = dict() for i in range(len(cols)): for j in range(i+1): info[(i,j)] = ih.get_info(i,j,cols,pD,v=False) for i,j in sorted(info.keys()): print i,j,round(info[(i,j)],3)

Mutual information

I want to talk about a really nice paper from Michael Laub's lab at MIT (Skerker et al 2009 PMID 18555780). It'll give us an opportunity to exercise our matplotlib skills.

We're going to try to recreate Fig 1, which is visible in the PubMed page, or you can get the original paper from the link to Cell.

Two-component signal transduction systems are ubiquitous in bacteria (wikipedia). The canonical design consists of a membrane-bound sensor (histidine) kinase (HK) and a cytoplasmic response regulator (RR). E. coli contains about 30 such pairs. The members of each pair have substantial specificity. The HK of the ntr system has specificity for its own RR, and likewise in the phoBR system, phoB has specificity for phoR. We may speak of a HK and its cognate RR.

For our purposes the important thing is that each system comprises (in the simplest design) two protein partners with complementary surfaces. These systems (a pair of proteins) are the products of ancient gene duplication events, and have since diverged over time. Amino acids at interacting sites are constrained to co-evolve in each pair.

If this sounds too vague or too complicated, consider an even simpler example: a stem of paired RNA residues in rRNA named H15.

Here is the H15 sequence in 1D (the parentheses indicate residues involved in pairing---see the link above for details):

(((( ((((( ))))) )))) TGCACAATGGGCGCAAGCCTGATGCA

And here is the inner stem drawn in 2D to show the base-pairing more directly:

TGGGC GTCCG

The base-pairing of this stem is more important to rRNA function than the identities of the bases. The result is that in some bacteria the identities of the central bases have been switched:

original co-evolved TGGGC TGCGC GTCCG GTGCG

Presumably this happened in 2 discrete steps, but I don't know of any examples where the intermediate state has been preserved. Maybe we should look for some, and it's undoubtedly been studied.

To quantify this kind of coevolution, we'll draw on a concept (and mathematical definition) called mutual information. The steps in the calculation will be:

make a multiple sequence alignment

compare column X and column Y

total number of sequences (length of each column) = c

for each residue x in column X calculate px

for each residue y in column Y calculate py

px is the probability of residue x in column X

We'll write the columns horizontally to save space.
Suppose column X and Y are:

X: AASSASSTTT Y: NMWWNTTKKS

For column X we have:
pA = 0.3 (3 A out of a total of 10 residues)
pS = 0.4
pT = 0.3

For column Y:

pK = 0.2
pM = 0.1
pN = 0.2
pS = 0.1
pT = 0.2
pW = 0.2

We pre-calculate these values for each column. When we calculate the information, we'll refer to the probabilities for column Y as q rather than p, to keep them straight from the p's for column X.

Now, we consider each pair of residues, one from column X and one from column Y. This pair is made up of residues in two interacting protein surfaces or rRNA chains, that may have co-evolved.

pxy is the number of sequences with x in column X and y in column Y

divided by c, the total number of pairs:

X: AASSASSTTT Y: NMWWNTTKKS


pAM = 0.1
pAN = 0.2
pST = 0.2
pSW = 0.2
pTK = 0.2
pTS = 0.1

Finally, to compute the mutual information for this pair of columns, we do this calculation for each individual pair of residues and then sum:

pAM * log (pAM / pAqM) = 0.1 * log (0.1 / (0.3 * 0.1)) = 0.0523
I would have used log₂, but Skerker et al used log₁₀, so I matched them.

Here is part of the output of the script below:

NMWWNTTKKS AASSASSTTT pKT 0.2 pK 0.2 qT 0.3 temp 3.33 final 0.1 pMA 0.1 pM 0.1 qA 0.3 temp 3.33 final 0.05 pNA 0.2 pN 0.2 qA 0.3 temp 3.33 final 0.1 pST 0.1 pS 0.1 qT 0.3 temp 3.33 final 0.05 pTS 0.2 pT 0.2 qS 0.4 temp 2.50 final 0.08 pWS 0.2 pW 0.2 qS 0.4 temp 2.50 final 0.08 info 0.47

temp is the result of the calculation inside the parentheses, above. Next time we'll apply this method to the data from Skerker.

from math import log def log2(n): return log(n)*1.0/log(2) def log10(n): return log(n)*1.0/log(10) # cache character probabilities for each column def make_prob_dict(cols): # input is a list of columns # c is the number of sequences in the alignment c = len(cols[0]) pD = list() for col in cols: char_kinds = list(set(col)) values = [col.count(k)*1.0/c for k in char_kinds] pD.append(dict(zip(char_kinds,values))) return pD def get_info(i,j,cols,pD,v=False): col1, col2 = cols[i], cols[j] if v: print col1 + '\n' + col2 # as before, c is the number of sequences c = len(col1) assert c == len(col2) info = 0 pairs = [col1[k] + col2[k] for k in range(c)] pL = sorted(list(set(pairs))) for p in pL: pXY = pairs.count(p)*1.0/c pX = pD[i][p[0]] pY = pD[j][p[1]] inside = (pXY * 1.0) / (pX * pY) if v: print 'p' + p, pXY, if v: print 'p' + p[0], pX, if v: print 'q' + p[1], pY, if v: print 'temp', '%3.2f' % round(inside, 2), outside = pXY * log10(inside) if v: print 'final', round(outside,2) info += outside if v: print 'info', round(info,2) return info if __name__ == '__main__': aln = 'AASSASSTTT\nNMWWNTTKKS\nGTSNTYRSTA\nGGGGGGGGGG' cols = aln.split('\n') pD = make_prob_dict(cols) info = dict() for i in range(len(cols)): for j in range(i): info[(i,j)] = get_info(i,j,cols,pD,v=True)

16S rRNA V regions, continued

As promised last time (here), this is the code to plot V regions from the sequences in the Enterobacteriales, obtained from RDP. If you compare to this post from the other day, you'll see we match pretty well.

The figure looks better with a wider window, but the limits are established more accurately with a smaller window. The threshold T needs to be adjusted based on the window size.

One detail: I did the redirect below to save the extreme values and then parsed them with the code at the end of the post.

python script.py > extreme_values.txt

script.py

import utils as ut from info import shannon import matplotlib.pyplot as plt data = ut.load_data('entero.txt') data = data.strip().split('>')[1:] EC = [e for e in data if 'X80725' in e][0] EC = ut.clean_fasta(EC)[1] data = [ut.clean_fasta(e)[1] for e in data] L = ut.make_count_list(data) pos = 0 R = range(1,1451) iL = list() for i,c in enumerate(EC): if c == '-': continue pos += 1 if not pos in R: continue cD = L[i] e = shannon(cD,'ACGT') iL.append(e) #print str(pos).ljust(3), c + ' ', #print ''.join([str(cD[k]).ljust(5) for k in 'ACGT']), #print round(e,2) aL = list() w = 20 T = 1.8 for i in range(len(iL)): j, k = i - w, i + w + 1 if j < 0: j = 0 if k > len(iL): k = len(iL) m = ut.mean(iL[j:k]) if m < T: print i+1 aL.append(m) plt.plot((1,1451),(T,T),lw=2,color='r',zorder=0) plt.scatter(R,aL) ax = plt.axes() ax.set_xlim(-5,1455) ax.set_ylim(0.8,2.05) plt.savefig('example.pdf')

analyze.py

import utils as ut data = ut.load_data('extreme_values.txt') L = [int(n) for n in data.strip().split('\n')] current = L.pop(0) print current, '-', while L: next = L.pop(0) if next != current + 1: print current print next, '-', current = next print next

16S rRNA V regions

I'm exploring ways to visualize the sequence diversity in 16S rRNA. I finally realized that to see the "V" regions clearly, I need an alignment in which the V regions can align, that is, the sequences need to be closely enough related. To get such a set, I went to RDP (Browsers), and grabbed 187 type sequences from the Order Enterobacteriales. In the next post, I'm actually going to generate the graphic above (it's a tease). The figure shows the information content of the sequence set computed for a window sliding across the length of the gene. The troughs are V regions.

But first, I thought I would post some utility code separately. Here is a script that imports the utility functions shown in the last part, below, and exercises them on the enteric sequences:

import utils as ut data = ut.load_data('entero.txt') data = data.strip().split('>')[1:] seqs = [ut.clean_fasta(e)[1] for e in data] L = ut.make_count_list(seqs,kL='ACGT') for D in L[75:80]: print D

Here is what it prints:

> python rdp.py {'A': 174, 'C': 2, 'T': 3, 'G': 1} {'A': 51, 'C': 2, 'T': 3, 'G': 123} {'A': 3, 'C': 172, 'T': 2, 'G': 3} {'A': 127, 'C': 5, 'T': 1, 'G': 47} {'A': 1, 'C': 76, 'T': 3, 'G': 99}

The first module contains a function to return the Shannon entropy for a distribution:

info.py

from math import log def f_logf(f): if f == 0: return 0 return f*log(f)*1.0/log(2) def shannon(cD,kL=None): if not kL: kL = cD.keys() L = [cD[k] for k in kL] S = sum(L) if S == 0: raise ValueError('Empty list') L = [n*1.0/S for n in L] L = [f_logf(f) for f in L] return 2 + sum(L)

The second module contains some functions that I use all the time. I put this in my site-packages directory:

utils.py

def load_data(fn): FH = open(fn,'r') data = FH.read().strip() if '\r\n' in data: data = s.replace('\r\n', '\n') FH.close() return data def reverse_complement(seq): import string.maketrans tt = maketrans('ACGT','TGCA') return seq[::-1].translate(tt) def write_data(fn,s): FH = open(fn,'w') FH.write(s) FH.close() def clean_fasta(s): title, seq = s.strip().split('\n',1) seq = ''.join(seq.strip().split()) return title, seq def make_count_list(L,kL=None,upper=True): # L is a list of strings (seqs) # each str the same length # count all chars assert len(set([len(s) for s in L])) == 1 iL = L[:] if upper: iL = [e.upper() for e in iL] if not kL: kL = sorted(list(set(''.join(iL)))) rL = list() R = range(len(iL[0])) for i in R: temp = [e[i] for e in iL] counts = [temp.count(k) for k in kL] rL.append(dict(zip(kL,counts))) return rL

Human Oral Microbiome Database

Bruce Paster and colleagues at the Forsyth Institute pioneered the molecular classification of oral microbes. There is a curated database called HOMD (Human Oral Microbiome Database) and also a BLAST server. I thought we could take a quick look at the database, and exercise it by comparing the taxonomic assignments to those we get using the RDP classifier accessed through QIIME.

The standard HOMD 16S rRNA RefSeq database (version 10.1, 2010-02-08) can be downloaded from this page. The filename is oral16S_20100208_9up.txt. It is formatted for Windows. (Always something to look out for). Just do:

data = data.replace('\r\n','\n')

or split on '>' and then strip(), ignoring the first (empty) value. It contains 755 sequences.

The phylogenetic assignments are in files linked here. We'll just use the first one, called homd_taxonomy_table.txt.

The taxonomy table has 625 entries in 27 columns. The zeroth column is labeled HOT_ID (001 through 850)---which look like integers and appear to be unique. The next 6 columns (L[1:7]) contain the phylognetic classification. Then comes a column (L[7]) labeled Species which sometimes looks like "sp. oral taxon 123", and sometimes looks like "invisus". Finally, four columns over is one labeled 16S_rRNA, which may have one or more Genbank IDs (separated by '|'), or perhaps something like "To be submitted".

The sequences have title lines like:

>001A28SC| Bartonella sp. | Oral Taxon 001 | Strain A28SC | GQ422708 | 0 | U

where the 001 is appears to be the oral taxon number, followed by the strain or clone number, although sometimes it is part of the Genbank ID, and sometimes it's the first few letters of the genus and species. Sometimes, an underscore follows the oral taxon number, but usually not.

>002_3433| Caulobacter sp. | Oral Taxon 002 | Strain TWE165 | DQ493433 | 2 | U

It appears that the rule was to make the title (the part before the first '|') contain 8 characters, and if the strain / clone identifier naturally contains 5, they dropped the underscore.

It seems simple enough to take the first 3 characters of the sequence's title line to map the sequences to the taxonomic information, and see how that works out. A possible difficulty is that a substantial fraction of these are not unique. There are four sequences labeled 058:

>058_8632| Streptococcus sp. | Oral Taxon 058 | Clone C3MLM097 | AY278632 | 251 | U >058BM035| Streptococcus sp. | Oral Taxon 058 | Clone BM035 | AY005043 | 251 | U >058BW009| Streptococcus sp. | Oral Taxon 058 | Clone BW009 | AY005042 | 251 | U >058CH016| Streptococcus sp. | Oral Taxon 058 | Clone CH016 | AY005044 | 251 | U

The corresponding entry of the taxonomy table is

058 Bacteria Firmicutes Bacilli Lactobacillales Streptococcaceae Streptococcus sp. oral taxon 058 AY005042|AY005043|AY005044|AY278632 251 0.72 26 Oral Clone BM035|Oral Clone BW009|Oral Clone C3MLM097|Oral Clone CH016 0 251 0 SEQF1704

Since this is just a demo, I'm going to remove duplicates from the sequences. That should leave us with 618. In the process we'll collapse the sequence title down to just those three characters, and rename the file 'homd.txt'.

The next step is to parse the phylogenetic information. That's also elementary so I won't show it either. We grab columns 0, 2 through 7, and 11.

The next thing to do is run the RDP classifier (using QIIME's script, see my post here).

> assign_taxonomy.py -i homd.txt -m rdp -o tax

Those results look like this:

001 Root;Bacteria;Proteobacteria;Alphaproteobacteria;Rhizobiales;Bartonellaceae;Bartonella 1.000

where the taxonomy is given in 7 levels starting with "Root". Some designations are quoted: "Lactobacillales", for example. The last value is a measure of the confidence of the assignment. I will strip the first 2 levels from the RDP taxonomy (just as a I stripped "Bacteria" from the HOMD taxonomy).

Now we just load the two data sets. Look at the first few entries:

001 homd: Proteobacteria:Alphaproteobacteria:Rhizobiales:Bartonellaceae:Bartonella rdp: Proteobacteria:Alphaproteobacteria:Rhizobiales:Bartonellaceae:Bartonella 002 homd: Proteobacteria:Alphaproteobacteria:Caulobacterales:Caulobacteraceae:Caulobacter rdp: Proteobacteria:Alphaproteobacteria:Caulobacterales:Caulobacteraceae 003 homd: Proteobacteria:Alphaproteobacteria:Sphingomonadales:Sphingomonadaceae:Sphingomonas rdp: Proteobacteria:Alphaproteobacteria:Sphingomonadales:Sphingomonadaceae:Sphingomonas

Now we look for differences. There seems to be an issue with the TM7, so we'll skip those (there are only about 10 anyway). Of 618 examined, there were 31 with discrepancies (more than I expected). I thought perhaps these would be ones that RDP was uncertain about, but that doesn't seems to be the case generally. Here are a few:

318 homd: Bacteroidetes:Flavobacteria:Flavobacteriales:Bacteroidetes[F-7]:Bacteroidetes[G-7] rdp: Bacteroidetes:Sphingobacteria:Sphingobacteriales:Sphingobacteriaceae:Pedobacter 504 homd: Firmicutes:Mollicutes:Acholeplasmatales:Acholeplasmatales[F-1]:Acholeplasmatales[G-1] rdp: o 507 homd: Bacteroidetes:Bacteroides:Bacteroidales:Bacteroidetes[F-5]:Bacteroidetes[G-5] rdp: Bacteroidetes:o 606 homd: Firmicutes:Mollicutes:Mycoplasmatales:Mycoplasmataceae:Mycoplasma rdp: Tenericutes:Mollicutes:Mycoplasmatales:Mycoplasmataceae:Mycoplasma:o

This is all such basic use of Python that I won't post the code. The next step is to use the HOMD data to explore map variation along the 16S rRNA sequence.

Mapping 16S rRNA structure to sequence

There has been a lot of discussion recently about which region(s) of bacterial 16S rRNA are the best (most informative, easily amplified, and short) for phylogenetic inference, particularly when using Illumina reads (e.g. Schloss 2010 PMID 20628621).

I have two different pictures in my head when thinking about the 16S rRNA sequence, the linear one of the sequence, and the two-dimensional one of secondary structure. I thought it might be useful to map the secondary structure onto the primary sequence. I've used the E. coli rrnB gene sequence as the reference. Error-correction is welcome.

The colors have no special meaning except that blue and red were chosen for short-range interactions.

The authoritative reference seems to be Neefs 1990 PMID 1692117. Here is a birds-eye-view. The secondary structure shown in these figures is from Case et al 2007 PMID 17071787.

I found a nice introductory post about the topic (here).

QIIME (1)

Yesterday and this morning, I've been working on installing QIIME and doing the tutorial. It stands for: Quantitative Insights into Microbial Ecology (link). It's yet another major software project out of Rob Knight's lab.

I'd like to post on the issues/questions/answers I run into as I work through it. It'll also give me an opportunity to get back to Bioinformatics (e.g. multiple sequence alignments).

The first thing is that they recommend you use a virtual machine (Virtual Box) and install the QIIME Virtual Box on it, which is a huge download that packages Ubuntu Linux with all of the QIIME dependencies, correctly configured. Sounds like a great idea.

I started by trying this. I had to download overnight on my machine at work (> 2 GB), then buy a new thumb drive to bring the file home, since I don't have the recommended amount of memory for the VM (1 GB) on that machine. But I ran into trouble---basically much aggravation dealing with the different keys on Linux, but also the Virtual Box and extra work trying to figure out copy/paste, moving files over, etc. The killer was when the VM prompted for an administrator's password, and of course I don't have it. There shouldn't be one needed..

[ After a search of the Forum the mystery word is revealed to be: qiime. I shoulda guessed it. ]

So I decided instead to deal with the dependencies for QIIME. I got almost all of them (like 24 or so), but for the "default pipeline" you really don't need so many. Particularly if you already have PyCogent and matplotlib installed, as I do, I would recommend just going down the list. It's not difficult, and it's better to be in an environment where you're comfortable working.

There was really no trouble configuring them (at least the essentials).

#1: Python 2.6

I have it, since I'm on OS X. I ran into a little trouble because (from playing with MacPorts) I still had /opt/local/bin on my path as changed in .bash_profile. The QIIME install checks this variable and modifies its scripts (the she-bang thing), so I had to fix that later on.

#2: Numpy 1.3.0

I have this too.

>>> import numpy >>> print numpy.__version__ 1.3.0

Not sure what the default OS X version is, but easy_install will fix you up quickly.

#3: PyCogent

Of course, this has its own dependencies. But we've been through that before (here, here, here).

#4: MUSCLE

This is not one of the primary dependencies, but it's an old friend. I upgraded to the latest and greatest version, but it caused an error in the QIIME tests. Luckily I still had muscle3.6_src around. As usual, I put it in ~/Software and do:

ln -s ~/Software/muscle3.6_src/muscle ~/bin/muscle

and I put ~/bin on my $PATH in .bash_profile:

export PATH=$HOME/bin/:$PATH

#5: MAFFT

It's an installer. Piece of cake. Check it:

> mafft

#6 uclust

See the QIIME install notes for this. This is the work of Robert Edgar (also the creator of MUSCLE), who is a talented programmer, but also a businessman. Hence, no source code, and no 64-bit for some things. But there's a download link on the QIIME site. The only problem is that Safari put a .txt extension on the download and I thought there was a problem, went off on a wild goose chase, and got a version that is not new enough. Enough said. Put it in ~/Software and link as usual

> uclust uclust v1.2.21q

#7 PyNAST

(download)

unpack and move to Software

python setup.py install cd tests python all_tests.py

After installing MUSCLE and MAFFT passes all tests except:

AssertionError: Got (DnaSequence(ACGTACG... 23), DnaSequence(ACGTACG... 23)), but expected (DnaSequence(ACGTACG... 23), DnaSequence(ACGTACG... 23))

That's pretty silly.

#8 Greengenes files

greengenes core set data file (fasta)
greengenes alignment lanemask file (txt)

Where to put them? I put them in qiime_tutorial, see below, but will eventually want them in Qiime somewhere.

#9 FastTree

fasttree 2.1.0 (src)

gcc -Wall -O3 -finline-functions -funroll-loops -o FastTree -lm FastTree-2.1.0.c

move to Software and link

#10 Java

I don't need to add the Java runtime. I got it from Apple a while ago (here):

> java -version java version "1.6.0_22" Java(TM) SE Runtime Environment (build 1.6.0_22-b04-307-10M3261) Java HotSpot(TM) 64-Bit Server VM (build 17.1-b03-307, mixed mode)

#11 RDP Classifier

rdp_classifier-2.0.1 (src)

in .bash_profile:

export RDP_JAR_PATH=$HOME/Software/rdp_classifier/rdp_classifier-2.0.jar

Check for mismatch between the name of the classifier file. What I actually got was named 2.0 not 2.0.1.

#12 QIIME

That's it for the default pipeline..

svn co https://qiime.svn.sourceforge.net/svnroot/qiime/trunk

I followed their advice and put the install scripts in a special place (but I wouldn't do that again).

python setup.py install --install-scripts=~/bin/qiime/bin/ cd tests python all_tests.py

The problem with the special directory is you will need a .qiime_config file (see the docs)

Make sure to put it in ~/.qiime_config.

Also not needed yet but I'll list them here anyway:

#13 BLAST

blast-2.2.22

We've done that one before. The BLASTMAT variable must point to the NCBI data directory. In .bash_profile:

export BLASTMAT=$HOME/Software/blast-2.2.22/data

#14 Infernal

(download)

./configure make make check sudo make install

#15 R

Run R and do:

install.packages('randomForest')

I upgraded to R 2.12.1 from 2.10.0, but it has some issue. Trying to get the ape package from CRAN hangs the app. I backed off to 2.11.1 (see here).

Still to come:

10 more bits of software.

But they're not needed for the first part of the tutorial, so we should just do that first. Oh..I also have Sphinx installed (from the PyCogent instructions), so the documentation got built. But I've just been working from the web version anyway. On to the fun stuff.

Exploring mitochondrial genomes 2

Continuing with the problem from the last post (here), we want to construct a plot similar to what the ape book has on p. 60 (Fig 3.5), which shows the nucleotide composition of the mitochondrial genes from the species we're looking at. The code shown below is in the same script with the listing from last time, and depends on it.

In the first part, we organize the sequence data into a dict keyed by gene name, where the value is a list of all the sequences for that gene in the data. In the second part we go through that dict, pull out the sequences for each of our target genes and count all the A, C, G, T. The data structure is a nested one: a dict of dicts. The third part is somewhat more complicated. We want to plot not the total counts for each nt, but the fraction of each. The first line does the calculation:

values = [nD[nt]*1.0/S for nt in 'atcg'] L = [sum(values[:j]) for j in range(1,5)] for z,c,f in zip([4,3,2,1],colors,L): b = plt.bar(i+1,f,width=1,color=c,zorder=z)

I chose the ATCG because then we can see at a glance what the GC content is for each gene. If we wanted to get slightly fancy we could sort on that. We're going to do a bar plot, so the value for each bar is not its fraction, but its fraction plus the sum of the fractions of all the nucleotides plotted previously to it. The plt.bar call looks a little convoluted: the purpose is to plot the bars in the correct z-order. In reality, we plot the bars one on top of the other, but we don't see what's underneath. We also tweak matplotlib into plotting the gene names as the tick labels.

The figure could use a key relating the colors to the nucleotides, but I'll leave that as an exercise for the reader.

Code listing:

D = dict() targets = ['ATP6','ATP8','COX1','COX2','COX3', 'CYTB','ND1','ND2','ND3','ND4','ND4L', 'ND5','ND6'] for gene in targets: D[gene] = [] for item in rest.strip().split('\n\n'): if 'Sequence does not exist' in item: continue title,seq = item.strip().split('\n',1) i = title.find('(') j = title.find(')') gene = title[i+1:j] D[gene].append(seq) #--------------------------------- rD = dict() for k in targets: nD = dict(zip('acgt',[0]*4)) for seq in D[k]: for nt in 'acgt': nD[nt] += seq.count(nt) rD[k] = nD #--------------------------------- colors = ['0.2','0.4','0.6','0.8'] iL = list() for i,k in enumerate(targets): iL.append(i+1) nD = rD[k] S = sum(nD.values()) values = [nD[nt]*1.0/S for nt in 'acgt'] L = [sum(values[:j]) for j in range(1,5)] for z,c,f in zip([4,3,2,1],colors,L): b = plt.bar(i+1,f,width=1,color=c,zorder=z) ax = plt.axes() ax.set_xlim(1,len(targets)+1 ) ax.xaxis.set_ticks([i+0.5 for i in iL]) ax.xaxis.set_ticklabels(targets) plt.savefig('example.png')

Exploring mitochondrial genomes

Following the ape book (Phylogenetics in R; starting on p. 55), they use a data set in a file named mammal_mtGenome.fasta. It's not part of the datasets with the ape package (too big), but a Google search led me to versions of Ch. 3 on the web, so you can follow along even if you don't have the book.

The example uses a database called OGRe---Organellar Genome Retrieval system (here). It has 1244 organisms! What's nice is that they've organized the data by genes and they have the nucleotide sequences.

If we go to select species from taxonomy and expand the phylogeny: Chordata > Vertebrata > Mammalia, we can get a checkbox for all mammals, then at the bottom of the page do display sequences and from the table at the very bottom choose these genes: rRNAs plus protein encoding genes: ATP6 .. ND6. I saved the download in 'sequences.fasta' and it's 2.7 MB.

The data file contains a header, followed by a section with entries like:

#HOMSAPMIT : _Homo sapiens_ (human) : MITOCHONDRION

after that we find the sequences of the genes, grouped by gene (all the ATP6 genes first, then ATP8 all the way to ND6 and beyond). Those entries look like this:

>GALPYRMIT(ATP6) atgaac..

or

>ORYCUNMIT(12SrRNA)Sequence does not exist

A double newline ('\n\n') is found between successive "does not exist" entries, but '\n\n\n' between other entries. There's a total of 233 species. Parsing the data file is a fairly trivial Python exercise. We use the fact that the organism code name is the first element, that the genus contains '_' as its first letter while the species contains '_' as its last letter, and that the common name is enclosed in parentheses. We sort and print the results. Next time, we'll do something with the data. Here we just print the first few lines:

Acinonyx jubatus ACIJUBMIT cheetah Ailuropoda melanoleuca AILMELMIT giant panda Ailurus styani FULSTYMIT red panda Ammotragus lervia AMMLERMIT aoudad Anomalurus GP-2005 ANOSP1MIT scaly-tailed squirrel Arctocephalus forsteri ARCFORMIT New Zealand fur seal Arctocephalus pusillus ARCPUSMIT cape fur seal Arctocephalus townsendi ARCTOWMIT Guadalupe fur seal Artibeus jamaicensis ARTJAMMIT Jamaican fruit-eating bat

Code listing:

import matplotlib.pyplot as plt fn = 'sequences.fasta' FH = open(fn,'r') data = FH.read().strip() FH.close() results = list() header, defs, rest = data.split('\n\n',2) for line in defs.strip().split('\n'): L = line.split() code = L[0][1:] genus = [w for w in L if w[0] == '_'][0][1:] species = [w for w in L if w[-1] == '_'][0][:-1] i = line.find('(') j = line.find(')') name = line[i+1:j] results.append((genus,species,code,name)) def test(): print len(results) results.sort() for g,s,c,n in results: print (g + ' ' + s).ljust(30), c, ' ', n test()

Likelihood in Seven Pages

I just re-read a "gentle introduction" to likelihood by Peter Foster. He is here, and I got the paper here. It's called

The Idiot's Guide to the Zen of Likelihood in a Nutshell in Seven Days for Dummies, Unleashed

Now that I've got the idea about Q and P matrices (e.g. here) and a bit of practice with Maximum Likelihood for phylogenetic models (here), it makes perfect sense. Feels like I'm making progress.

Highly recommended.

Tracking evolution of coding sequences 2

I grabbed some sequences (based on the figure in Page & Holmes): ATP synthase alpha subunit, the medium wavelength opsin, terminal transferase, beta globin, and IL6. I got both the human and mouse CDS for each. I did it by hand, and include the files for each in the zip. I used MUSCLE to do the alignment, and then fixed two obvious problems with the IL6 alignment manually.

Finally, I ran the script evaluate.py which does what we said we were going to do last time. We count non-synonymous sites, and synonomous and non-synonomous differences between each pair of sequences.

There isn't a lot to say about how it works. Except of course that the whole point of the exercise is to notice the huge difference in non-synonymous changes allowed in IL6, say, compared to the others.

Here are the results:

gene syn non non/syn ATP_synthase 0.418 0.013 0.031 opsin 0.35 0.066 0.189 TdT 0.539 0.122 0.226 beta_globin 0.379 0.142 0.375 IL6 0.683 0.519 0.76

Tracking evolution of coding sequences

I have two posts exploring calculation of the ratio of synonymous to non-synonymous changes when comparing coding sequences. It's more complicated than I thought. This part sets up the utilities we'll need to process the sequences next time. Our goal is to be able to calculate data for a plot like the graphic above, which is from Page & Holmes.

The GeneticCode class constructs the "universal" genetic code in a few lines using list comprehensions and a hard-coded string of one-letter abbreviations for the amino acids.

aa = 'FFLLSSSSYY**CC*WLLLLPPPPHHQQRRRR' +\ 'IIIMTTTTNNKKSSRRVVVVAAAADDEEGGGG'

Objects of this class also contain a reverse dict mapping amino acids to codons, and a third dictionary with codons as keys and a list of synonymous codons as the value. It uses a special cmp function to sort codons in order 'TCAG', i.e. TTT to GGG. I added a function hamming which calculates the Hamming distance between two codons, and a function to construct a set of utility dictionaries. These contain all the codons at a specified hamming distance from each codon. Like this:

nD1 CTT ['TTT', 'CTC', 'CTA', 'CTG', 'CCT', 'CAT'] .. TAG ['TTG', 'TCG', 'TAT', 'TAC', 'TAA', 'TGG'] .. .. nD2 CTT ['TTC', 'TTA', 'TTG', 'TCT', 'TAT', 'TGT'] .. TAG ['TTT', 'TTC', 'TTA', 'TCT', 'TCC', 'TCA'] .. .. nD3 CTT ['TCC', 'TCA', 'TCG', 'TAC', 'TAA', 'TAG'] .. TAG ['CTT', 'CTC', 'CTA', 'CCT', 'CCC', 'CCA'] .. ..

The calculation of synonymous and non-synonymous subtitutions proceeds in two parts (Yang). The code to do this is in the module SN_utils.py. The first part is to assess, for each codon one at a time, what fraction of the 9 codons that have a hamming distance of 1 to it, are either synonymous or not. Nonsense codons are discarded before doing the calculation.

For example, CTT has 9 neighbors, none of them stop codons, and 3 of them synonymous. We calculate (3/9) * 3 = 1.0 and record a tuple containing that value as well as 3 minus the value => (1.0, 2.0).

On the other hand, two of TAT's neighbors are stop codons, and one neighbor is synonymous, so we record a value of (1/7) * 3 = 0.429 and 3 minus that => (0.429, 2.571). The values are saved in a dict keyed by each codon.

When we come to evaluate sequences we'll use this info as follows: for each sequence we accumulate the sums of the values for the list of codons (the two elements of the tuple, S and N, are summed separately). Then we average those sums over the two sequences that we're comparing. The two averages will be the denominators of the final ratios for S and N.

The second part of the calculation involves comparisons between two codons. There are four cases, corresponding to a Hamming distance of 0, 1, 2, or 3.

0:
(The same codon), so there is nothing to do, as we just ignore these pairs.

1:
Also simple, we record the tuple (1.0, 0.0) if they are synonomous and (0.0, 1.0) if not.

2:
For this case, there are two positions which need to be changed to convert one codon into the other (at least in the simplest pathway, which is all that we look at). We consider each of the single changes as candidate intermediates. This results in 2 possible paths, with 2 changes for each: C1 => I1 => C2, and C1 => I2 => C2. For each of the 4 single changes, then, we record (0.5, 0.0) if they are synonymous and (0.0, 0.5) otherwise. We sum the tuples for the final value. For example, for the pair CTT CAA, there are two paths (asterisks mark the base that just changed):

* * * * CTT => CAT => CAA CTT => CTA => CAA Leu His Gln Leu Leu Gln

One of the four changes (CTT => CTA) is synonymous, while the other 3 are not. We record the final value (0.5, 1.5).

3:
This is the more complex case. Here there 3 differences and hence 6 possible paths. For the pair GGC and CTT, the paths are:

* * * * * * GGC => CGC => CTC => CTT GGC => CGC => CGT => CTT Gly Arg Leu Leu Gly Arg Arg Leu * * * * * * GGC => GTC => CTC => CTT GGC => GTC => GTT => CTT Gly Val Leu Leu Gly Val Val Leu * * * * * * GGC => GGT => CGT => CTT GGC => GGT => GTT => CTT Gly Gly Arg Leu Gly Gly Val Leu

Of the 18 changes, 6 are synonymous. We record the final value (3.0, 6.0). The debug output shows the program processing this pair:

next pair: GGC CTT d = 3 test: CGC process: d = 1 GGC CGC value = (0.0, 1.0) process: d = 2 CGC CTT test intermediates: CTC syn CTT count 1 CGT syn CGC count 2 value = (1.0, 1.0) test: GTC process: d = 1 GGC GTC value = (0.0, 1.0) process: d = 2 GTC CTT test intermediates: CTC syn CTT count 1 GTT syn GTC count 2 value = (1.0, 1.0) test: GGT process: d = 1 GGC GGT value = (1.0, 0.0) process: d = 2 GGT CTT test intermediates: GTT count 0 CGT count 0 value = (0.0, 2.0) 3.0 6.0

The values for each pair are saved in a second dict. We'll exercise that next time. The zipped project files are on Dropbox (here).

Models of sequence evolution 4

Just to finish up with this topic. I wrapped the linear algebra routines into a function, then (i) generate a TN93 Q matrix for a particular choice of rates and π_N, find P for that Q at a particular short t (0.1), multiply P by itself many times to see the convergence for the π values, and finally compare the P we got by matrix stuff to the P we obtain using the equations from Yang. It looks good to me. Time to move on.

Q [[-0.7 0.6 0.04 0.06] [ 0.4 -0.5 0.04 0.06] [ 0.04 0.06 -0.4 0.3 ] [ 0.04 0.06 0.2 -0.3 ]] evals [ -1.10000000e+00 -2.00000000e-01 -1.62026873e-17 -6.00000000e-01] L [[ 0.89583414 0. 0. 0. ] [ 0. 0.98019867 0. 0. ] [ 0. 0. 1. 0. ] [ 0. 0. 0. 0.94176453]] U [[ -8.32050294e-01 5.00000000e-01 -5.00000000e-01 -6.12373075e-17] [ 5.54700196e-01 5.00000000e-01 -5.00000000e-01 -4.08248716e-17] [ -2.98507646e-18 -5.00000000e-01 -5.00000000e-01 -8.32050294e-01] [ 4.00977456e-18 -5.00000000e-01 -5.00000000e-01 5.54700196e-01]] iU [[ -7.21110255e-01 7.21110255e-01 -1.36437676e-17 -5.30723331e-17] [ 4.00000000e-01 6.00000000e-01 -4.00000000e-01 -6.00000000e-01] [ -4.00000000e-01 -6.00000000e-01 -4.00000000e-01 -6.00000000e-01] [ -1.60118642e-16 -2.40177963e-16 -7.21110255e-01 7.21110255e-01]] P [[ 0.93354022 0.05655912 0.00396027 0.0059404 ] [ 0.03770608 0.95239326 0.00396027 0.0059404 ] [ 0.00396027 0.0059404 0.96109845 0.02900088] [ 0.00396027 0.0059404 0.01933392 0.97076542]] converges to: [[ 0.2 0.3 0.2 0.3] [ 0.2 0.3 0.2 0.3] [ 0.2 0.3 0.2 0.3] [ 0.2 0.3 0.2 0.3]] explicit P: [[ 0.93354022 0.05655912 0.00396027 0.0059404 ] [ 0.03770608 0.95239326 0.00396027 0.0059404 ] [ 0.00396027 0.0059404 0.96109845 0.02900088] [ 0.00396027 0.0059404 0.01933392 0.97076542]]

code:

from __future__ import division import math import numpy as np def JC69_Q(): a = 0.01/3 Q = np.asarray( [ -3*a, a, a, a, a, -3*a, a, a, a, a, -3*a, a, a, a, a, -3*a ]) Q.shape = (4,4) return Q def TN93_values(): pi = [0.2,0.3,0.2,0.3] T,C,A,G = pi Y = T + C R = A + G a1 = 2 # Y transitions a2 = 1 # R transitions b = 0.2 # transversions return T,C,A,G,Y,R,a1,a2,b def TN93_Q(): T,C,A,G,Y,R,a1,a2,b = TN93_values() Q = np.asarray( [ -(a1*C + b*R), a1*C, b*A, b*G, a1*T, -(a1*T + b*R), b*A, b*G, b*T, b*C, -(a2*G + b*Y), a2*G, b*T, b*C, a2*A, -(a2*A + b*Y) ]) Q.shape = (4,4) return Q def TN93_P(t): # order TCAG T,C,A,G,Y,R,a1,a2,b = TN93_values() e2 = math.e**(-b*t) e3 = math.e**(-(R*a2 + Y*b)*t) e4 = math.e**(-(Y*a1 + R*b)*t) TY,CY,AR,GR = T/Y,C/Y,A/R,G/R M = [ T + TY*R*e2 + CY*e4, C + CY*R*e2 - CY*e4, A*(1 - e2), G*(1 - e2), T + TY*R*e2 - TY*e4, C + CY*R*e2 + TY*e4, A*(1 - e2), G*(1 - e2), T*(1 - e2), C*(1 - e2), A + AR*Y*e2 + GR*e3, G + GR*Y*e2 - GR*e3, T*(1 - e2), C*(1 - e2), A + AR*Y*e2 - AR*e3, G + GR*Y*e2 + AR*e3 ] P = np.array(M) P.shape = (4,4) return P def convert(Q,t,debug=False): if debug: print 'Q', '\n', Q evals, evecs = np.linalg.eig(Q) if debug: print 'evals', '\n', evals L = np.diag([math.e**(k*t) for k in evals]) if debug: print 'L', '\n', L U = evecs if debug: print 'U', '\n', U iU = np.linalg.inv(U) if debug: print 'iU', '\n', iU P = np.dot(U,np.dot(L,iU)) if debug: print 'P', '\n', P return P def long_time(P,N=5000): M = P for i in range(N): M = np.dot(P,M) return M if __name__ == '__main__': #Q = JC69_Q() Q = TN93_Q() t = 0.1 P = convert(Q,t,debug=True) print 'converges to:' print long_time(P) print 'explicit P:' print TN93_P(t)

Models of sequence evolution 3

I wrote a couple of long posts about models of sequence evolution (here and here), in which (and it may not be obvious) we solved the problem of how to obtain:

P = eQt

for any t, by getting eigenvalues and eigenectors, using both R and Python. That's outstanding.

If you look carefully, you might notice something. There are just two terms in Q for the JC69 model:

-0.01 (on diagonal) and 0.003333 (off diagonal)

P has only two terms as well and they turn out to be (for t = 0.1):

0.999 and 0.000333

That is, the off-diagonal terms are apparently just multiplied by t, and naturally the diagonal is 1 - 3 times the other term. There was really no need for all the exponentiation of matrices, etc.

How does this happen? The formulas for the two terms in JC69 are:

PNN(t) = 1/4 + 3/4 exp(-4λt) PNM(t) = 1/4 - 1/4 exp(-4λt)

But going back to the infinite series expansion for e (or just remembering the approximation for small x):

ex = 1 + x1/1! + x2/2! + .. ex = 1 + x

So the exponential is just 1 - 4λt and we obtain:

PNN(t) = 1/4 + 3/4 (1 - 4λt) = 1 - 3λt PNM(t) = 1/4 - 1/4 (1 - 4λt) = λt

Let's take another look at the formulas for these terms in the TN93 model. Look at only the first row of P. We first consider the two exponential terms:

e2 = exp(-βt) e4 = exp[-(πYα1 + πRβ)t]

Approximating for small t:

e2 = 1 - βt e4 = 1 - (πYα1 + πRβ)t

Look what happens to P13 and P14:

P13 = &piA(1 - e2) = &piAβt

This is t times Q13. Exactly the same thing happens for P14. P12 is a bit more complicated, but it's just algebra:

P12 = πC + (πCπR/πY)e2 - (πC/πY)e4

The second term is:

πCπR/πY (1 - βt) πCπR/πY - (πCπR/πY)βt

The third term is:

πC/πY (-1 + πYα1t + πRβt)

We can see that (considering both the second and the third term) there are two copies of: (πCπR/πY)βt with opposite sign, leaving:

P12 = πC + πCπR/πY - πC/πY + πCα1t

The last term is t times Q12. So I'm betting that the rest will go away. We have:

πC + πCπR/πY - πC/πY

but

πY + πR = 1 πR - 1 = -πY

so

πC + πCπR/πY - πC/πY = πC + πC/πY(πR - 1) = πC - πC = 0

P12 = πCα1t = t Q12

How about that!

Models of sequence evolution 2

We have an equation to solve:

P = eQt

It looks easy enough, except that here Q is a 4 x 4 matrix. But linear algebra teaches us that if we can find the eigenvalues and eigenvectors of Q, we can form a succession of 3 matrices such that

Q = U Λ U-1

where &Lambda is a diagonal matrix. The upshot of that is that any function we want to do on Q is equal to:

f(Q) = U f(Λ) U-1

U and U^-1 are the right and left eigenvectors of Q, and moreover as the notation implies:

U U-1 = I

So, let's see where that takes us. Long ago, we had this matrix:

A = [ 5 2 ] [ 9 2 ]

Using the characteristic equation we found eigenvalues

λ = [8,-1]

We also found eigenvectors v, one for each eigenvalue, such that:

Av = λv [ 5 2 ] [ x ] = λ [ x ] [ 9 2 ] [ y ] [ y ]

For example:

[ 5 2 ] [ 2 ] = 8 [ 2 ] [ 9 2 ] [ 3 ] [ 3 ]

Checking our multiplication, rows x columns:

[ 16 ] = [ 16 ] [ 24 ] [ 24 ]

We scale v to be a unit vector:

[ 2/√13 ] = [ 0.5547 ] [ 3/√13 ] [ 0.8321 ]

And the second eigenvector is:

[ 1/√10 ] = [ 0.3162 ] [ -3/√10 ] [ -0.9487 ]

In R:

M = matrix(c(5,2,9,2),nrow=2,byrow=T) eigen(M) $values [1] 8 -1 $vectors [,1] [,2] [1,] 0.5547002 -0.3162278 [2,] 0.8320503 0.9486833

Any linear multiple of v is still an eigenvector, so it's OK that R has our v₂ x -1.

So far, so good. Now, Yang says that what we need to do is "decompose" our matrix into right and left eigenvectors. We have a right eigenvector, but now we need a new left eigenvector for each eigenvalue such that:

vA = vλ = λv

Luckily, we can use inverses. He says that we want:

Q = U Λ U-1

Where U is the right eigenvector of Q and the "rows of " U^-1 are the left eigenvector of Q. &Lambda is a diagonal vector with the eigenvalues λ. Here is the code separately (makes it easier for you to run it):

result = eigen(M) U = as.matrix(result$vectors,nrow=2) iU = solve(U) U %*% iU L = diag(result$values) L U %*% L %*% iU

and what it prints:

> result = eigen(M) > U = as.matrix(result$vectors,nrow=2) > iU = solve(U) > U %*% iU [,1] [,2] [1,] 1 0 [2,] 0 1 > L = diag(result$values) > L [,1] [,2] [1,] 8 0 [2,] 0 -1 > U %*% L %*% iU [,1] [,2] [1,] 5 2 [2,] 9 2

I would just note in passing that although svd sounds like what we want (singular value decomposition), it doesn't give the values we're looking for. I'm not sure why yet.

Let's try a phylogenetics example using a simple matrix like Jukes-Cantor's Q:

a = 0.01/3 m = c(-3*a, a, a, a, a, -3*a, a, a, a, a, -3*a, a, a, a, a, -3*a) Q = matrix(m,nrow=4)

result = eigen(Q) U = as.matrix(result$vectors,nrow=4) U iU = solve(U) iU U %*% iU L = diag(result$values) L U %*% L %*% iU

> result = eigen(Q) > U = as.matrix(result$vectors,nrow=4) > U [,1] [,2] [,3] [,4] [1,] -0.5 0.8660254 0.0000000 0.0000000 [2,] -0.5 -0.2886751 -0.1601282 0.8006408 [3,] -0.5 -0.2886751 -0.6133112 -0.5389954 [4,] -0.5 -0.2886751 0.7734393 -0.2616453 > iU = solve(U) > iU [,1] [,2] [,3] [,4] [1,] -5.000000e-01 -0.5000000 -0.5000000 -0.5000000 [2,] 8.660254e-01 -0.2886751 -0.2886751 -0.2886751 [3,] -7.698012e-18 -0.1601282 -0.6133112 0.7734393 [4,] 3.849006e-17 0.8006408 -0.5389954 -0.2616453 > U %*% iU [,1] [,2] [,3] [,4] [1,] 1.000000e+00 1.110223e-16 -5.551115e-17 1.110223e-16 [2,] 4.293802e-18 1.000000e+00 0.000000e+00 0.000000e+00 [3,] 1.173089e-17 0.000000e+00 1.000000e+00 0.000000e+00 [4,] 1.173089e-17 -2.775558e-17 2.775558e-17 1.000000e+00 > L = diag(result$values) > L [,1] [,2] [,3] [,4] [1,] -1.734723e-18 0.00000000 0.00000000 0.00000000 [2,] 0.000000e+00 -0.01333333 0.00000000 0.00000000 [3,] 0.000000e+00 0.00000000 -0.01333333 0.00000000 [4,] 0.000000e+00 0.00000000 0.00000000 -0.01333333 > U %*% L %*% iU [,1] [,2] [,3] [,4] [1,] -0.010000000 0.003333333 0.003333333 0.003333333 [2,] 0.003333333 -0.010000000 0.003333333 0.003333333 [3,] 0.003333333 0.003333333 -0.010000000 0.003333333 [4,] 0.003333333 0.003333333 0.003333333 -0.010000000

Check.

Now, what we need to do is to exponentiate L. That is, we do exp { λt } for each eigenvalue, and place the result in a diagonal matrix. The code:

t = 0.1 result$values k = exp(t*result$values) D = diag(k) P = U %*% D %*% iU P M = P for (i in 1:5000) { M = P %*% M } M

It prints:

> t = 0.1 > result$values [1] -1.734723e-18 -1.333333e-02 -1.333333e-02 -1.333333e-02 > k = exp(t*result$values) > D = diag(k) > P = U %*% D %*% iU > P [,1] [,2] [,3] [,4] [1,] 0.9990006664 0.0003331112 0.0003331112 0.0003331112 [2,] 0.0003331112 0.9990006664 0.0003331112 0.0003331112 [3,] 0.0003331112 0.0003331112 0.9990006664 0.0003331112 [4,] 0.0003331112 0.0003331112 0.0003331112 0.9990006664 > > M = P > for (i in 1:5000) { M = P %*% M } > M [,1] [,2] [,3] [,4] [1,] 0.2509532 0.2496823 0.2496823 0.2496823 [2,] 0.2496823 0.2509532 0.2496823 0.2496823 [3,] 0.2496823 0.2496823 0.2509532 0.2496823 [4,] 0.2496823 0.2496823 0.2496823 0.2509532

Looks good. How about doing it in Python?

Output:

Q [[-0.01 0.00333333 0.00333333 0.00333333] [ 0.00333333 -0.01 0.00333333 0.00333333] [ 0.00333333 0.00333333 -0.01 0.00333333] [ 0.00333333 0.00333333 0.00333333 -0.01 ]] evals [ -1.33333333e-02 4.33680869e-19 -1.33333333e-02 -1.33333333e-02] L [[ 0.99866756 0. 0. 0. ] [ 0. 1. 0. 0. ] [ 0. 0. 0.99866756 0. ] [ 0. 0. 0. 0.99866756]] U [[-0.8660254 0.5 -0.13387887 -0.2397665 ] [ 0.28867513 0.5 0.80991131 0.14009673] [ 0.28867513 0.5 -0.11709746 -0.62763021] [ 0.28867513 0.5 -0.55893497 0.72729998]] iU [[-0.8660254 0.08069103 0.58314599 0.20218839] [ 0.5 0.5 0.5 0.5 ] [-0. 0.84935676 -0.36809649 -0.48126026] [-0. 0.27697196 -0.85808026 0.58110829]] P [[ 9.99000666e-01 3.33111210e-04 3.33111210e-04 3.33111210e-04] [ 3.33111210e-04 9.99000666e-01 3.33111210e-04 3.33111210e-04] [ 3.33111210e-04 3.33111210e-04 9.99000666e-01 3.33111210e-04] [ 3.33111210e-04 3.33111210e-04 3.33111210e-04 9.99000666e-01]] M [[ 0.2509532 0.24968227 0.24968227 0.24968227] [ 0.24968227 0.2509532 0.24968227 0.24968227] [ 0.24968227 0.24968227 0.2509532 0.24968227] [ 0.24968227 0.24968227 0.24968227 0.2509532 ]]

Code listing:

from __future__ import division import math import numpy as np a = 0.01/3 Q = np.asarray( [ -3*a, a, a, a, a, -3*a, a, a, a, a, -3*a, a, a, a, a, -3*a ]) Q.shape = (4,4) print 'Q', '\n', Q evals,evecs = np.linalg.eig(Q) print 'evals', '\n', evals t = 0.1 L = np.diag([math.e**(k*t) for k in evals]) print 'L', '\n', L U = evecs print 'U', '\n', U iU = np.linalg.inv(U) print 'iU', '\n', iU P = np.dot(U,np.dot(L,iU)) print 'P', '\n', P M = P for i in range(5000): M = np.dot(P,M) print 'M', '\n', M

Models of sequence evolution 1

I've decided to work on Yang's Computational Molecular Evolution again. It's not an elementary book, but I like his style and clarity of explanation a lot. The goal is to understand models of sequence evolution, and more than that, to see how starting with the substitution-rate matrix Q, we can get to the transition probability matrix P. Once we have P, we can crank it as many times as we need to.

Q is a 4 x 4 matrix with the instantaneous rates for each nucleotide to change to another, or stay the same.

Everybody has a different order for the nucleotides! I usually use ACGT (because then I don't have to think twice about it), wikipedia uses AGCT some places and something else other places. Yang uses TCAG, and we'll stick with that since we're borrowing his whole exposition.

Usually, the columns of Q represent the instantaneous rate of going to a particular state, so row 1 is T => T (TT), TC TA and TG. There are exceptions (wikipedia some places, Dayhoff), but this is pretty standard.

I'm going to start with a more advanced model, called TN93 (after Tamura & Nei 1993 PMID 8336541). For TN93 we have 3 exchangeability parameters (α₁, α₂, β): these describe the rate at which either of the two transitions, or all transversions, occur.

α1 TC α2 AG β transversions

The model is time-reversible, so the exchangeability parameter for MN applies also to NM.

The equilibrium nucleotide frequencies are π_N. These are the frequencies to which a sequence would tend by many sequential applications of P. So here is Q. Notice the symmetry:

Q = - α1πC βπA βπG α1πT - βπA βπG βπT βπC - α2πG βπT βπC α2πA -

For example, the instantaneous rate of change from T to C is Q₁₂, which is α₁π_C. And that's counterintuitive: it involves π_C, even though that's where we're going.

There are two other rules. First, the nucleotide frequencies must sum to one: Σ πi = 1. Also the sum of rates in each row equals zero, so we calculate the diagonals accordingly.

Now, how do you get from Q to P? Here is the equation, which looks simple but definitely is not:

P(t) = eQt

But we can avoid using it for the moment. Yang gives an analytical solution (i.e. which is explicit in terms of the above parameters) for TN93. That's probably why he chose this example, because it's complex enough to be realistic but also has this kind of solution.

The eigenvalues (remember those?) are:

λ1 = 0 λ2 = -β λ3 = -(πRα2 + πYβ) λ4 = -(πYα1 + πRβ)

He defines e values:

e2 = exp(λ2t) = exp(-βt) e3 = exp(λ3t) = exp[-(πRα2 + πYβ)t] e4 = exp(λ4t) = exp[-(πYα1 + πRβ)t]

Then come the 16 terms for the P matrix..

row 1: πT + (πTπR/πY)e2 + (πC/πY)e4 πC + (πCπR/πY)e2 - (πC/πY)e4 πA(1 - e2) πG(1 - e2) row 2: πT + (πTπR/πY)e2 - (πT/πY)e4 πC + (πCπR/πY)e2 + (πT/πY)e4 πA(1 - e2) πG(1 - e2) row 3: πT(1 - e2) πC(1 - e2) πA + (πAπY/πR)e2 + (πG/πR)e3 πG + (πGπY/πR)e2 - (πG/πR)e3 row 4: πT(1 - e2) πC(1 - e2) πA + (πAπY/πR)e2 - (πA/πR)e3 πG + (πGπY/πR)e2 + (πA/πR)e3

Naturally, I wrote a simulation in Python. It chooses α1, α2, and β = 2, 1, 0.2. I have no idea yet if those are reasonable values.

But what you can see from the output is that, if we crank the P matrix a very long time, the nucleotide frequencies approach the values in the π vector, as they should.

P [[ 9.93036567e-01 5.96443233e-03 3.99600267e-04 5.99400400e-04] [ 3.97628822e-03 9.95024711e-01 3.99600267e-04 5.99400400e-04] [ 3.99600267e-04 5.99400400e-04 9.96011178e-01 2.98982117e-03] [ 3.99600267e-04 5.99400400e-04 1.99321411e-03 9.97007785e-01]] P * P [[ 9.86145739e-01 1.18582558e-02 7.98402131e-04 1.19760320e-03] [ 7.90550385e-03 9.90098491e-01 7.98402131e-04 1.19760320e-03] [ 7.98402131e-04 1.19760320e-03 9.92044626e-01 5.95936909e-03] [ 7.98402131e-04 1.19760320e-03 3.97291272e-03 9.94031082e-01]] 2500 rounds later: [[ 0.20134221 0.30201331 0.19865779 0.29798669] [ 0.20134221 0.30201331 0.19865779 0.29798669] [ 0.19865779 0.29798669 0.20134239 0.30201313] [ 0.19865779 0.29798669 0.20134209 0.30201344]]

from __future__ import division import math import numpy as np from numpy import dot # order TCAG pi = [0.2,0.3,0.2,0.3] T,C,A,G = pi Y = T + C R = A + G a1 = 2 # Y transitions a2 = 1 # R transitions b = 0.2 # transversions t = 0.01 e2 = math.e**(-b*t) e3 = math.e**(-(R*a2 + Y*b)*t) e4 = math.e**(-(Y*a1 + R*b)*t) TY,CY,AR,GR = T/Y,C/Y,A/R,G/R M = [ T + TY*R*e2 + CY*e4, C + CY*R*e2 - CY*e4, A*(1 - e2), G*(1 - e2), T + TY*R*e2 - TY*e4, C + CY*R*e2 + TY*e4, A*(1 - e2), G*(1 - e2), T*(1 - e2), C*(1 - e2), A + AR*Y*e2 + GR*e3, G + GR*Y*e2 - GR*e3, T*(1 - e2), C*(1 - e2), A + AR*Y*e2 - AR*e3, G + GR*Y*e2 + AR*e3 ] P = np.array(M) P.shape = (4,4) print 'P', '\n', P, '\n' M = dot(P,P) print 'P * P', '\n', M, '\n' print '2500 rounds later:' for i in range(2500): M = dot(P,M) print M