Dental project (6)

I want to show some more results from this project, namely, the UniFrac analysis. What I did for the paper was to cluster very closely related sequences (alignment > 450 and 0 or 1 mismatches), then upload them to RDP, which aligns the sequences as they are uploaded. The phylogenetic tree needs to be rooted, and I decided to use Thermotoga SL7 for this (Genbank AJ401017).

Rather than deal with the clustered OTUs for this post, I just uploaded all 1120 sequences, and carried out the analysis. The first time through (today) I forgot to include the outgroup! So that gives us a chance to see how much difference it makes.

Working in directory: temp

Check that seqs.fna from dental project dir has 1120 seqs

Rename to dental_1120.fna

Upload to RDP

Download as dental_1120_rdp.fna

Use R/ape to make a tree

setwd('Desktop/temp/rdp') library(ape) dna = read.dna('dental_1120_rdp.fna',format='fasta') tr = nj(dist.dna(dna)) plot(tr) write.tree(tr,'tree.txt')

Write a simple script to make the environment file

python write_env.py > environ.txt

It looks like this:

DAA_44 D_DAA DAA_43 D_DAA DQ_209 D_DQ ..

Upload tree and environment file to UniFrac

PCA (unweighted)

View as data table

Download data to pca_web.txt

Run the script below to plot the data using matplotlib.

Here it is:




From the UniFrac FAQ:

My tree was not rooted, but I was able to upload my file and perform an analysis. Are the results valid?

There is no way to tell based on a Newick string alone whether a tree is rooted or not. If an unrooted tree is input, UniFrac will usually assign an arbitrary root and allow you to perform the analysis on that tree. How the tree is rooted can affect the results of both UniFrac tests and the P test. You should redo the analysis with a tree that is rooted with an appropriate outgroup.

It turns out to be easy enough..go back to RDP and browse to find SL7 and then add it to the sequence cart. Repeat the download to dental_1120+_rdp.fna. Load the last 5 sequences into clustalx.app and look to make sure that SL7 is really properly aligned.

setwd('Desktop/temp/rdp') library(ape) dna = read.dna('dental_1120+_rdp.fna',format='fasta') tr = nj(dist.dna(dna)) > tr Phylogenetic tree with 1121 tips and 1119 internal nodes. Tip labels: DAA_44, DAA_43, DQ_209, DAA_45, DAA_40, DC_81, ... Unrooted; includes branch lengths. > grep('SL', tr$tip.label) > 1121

Root the tree appropriately and write it to disk

tr2=root(tr,1121) plot(tr2) write.tree(tr2,'rooted_tree.txt')

Go back to UniFrac



Repeat the PCA. You can look at the data in a spreadsheet app:



Now I plot it in matplotlib. The first image is what I plotted today for the rooted tree. The second is from the paper. Looks pretty good to me. Also, note some minor differences from the previous graphic where the tree we used was unrooted (and UniFrac rooted it for us however it does when it's not properly rooted).





plot_web.py

import sys import matplotlib.pyplot as plt from fileUtilities import load_data d = 0.5 fn = 'pca_web.txt' data = load_data(fn) data = data.strip().split('\n\n')[0] data = data.strip().split('\n')[1:] L = list() for e in data: name, x, y = e.split()[:3] x,y = float(x), float(y) x *= -1 name = name[2:] if name[1] in 'BCM': c = 'blue' else: c = 'red' plt.scatter(x,y,s=100,color=c) if name == 'DG': y += 0.03 plt.text(x+0.03,y-0.02,va='center', s=name[1:],color=c,fontsize=16) plt.plot((-d,d),(0,0),':',zorder=0) plt.plot((0,0),(-d,d),':',zorder=0) ax = plt.axes() ax.set_xlim(-d,d) ax.set_ylim(-d,d) plt.savefig('pca_web.png')

Dental project (5)

This post is one of a series (see dental project here or in the sidebar).

Last time I said I would show you how I make heatmaps these days. I've approached it several different ways over the past few years (R, Cocoa, matplotlib), but I think now that matplotlib is best, at least for me. Ultimately what I want is flexibility, and if you're a Python coder and you have matplotlib installed (as we've also discussed many times), then you'll have that. But I don't want to get into the technical details---and actually the script is a bit long, so I just put it (Heatmapper.py and its helper Preprocessor.py) into the zipped project files on Dropbox (here). The output from two different modes is at the bottom of the post. You just need a file data.csv in the same directory. It looks a little fuzzy and not as clean as I would like, but that's because there are so many samples, and partly because of the italic font. If you do savefig to a pdf file, and then blow it up, it looks great.

In this post I want to talk in a general way about the project and what I think it means. It began about four years ago, when we became aware that some folks in Dentistry at our school (WVU) were involved in a huge study of people from Appalachia (it's called COHRA). Poor oral health is a particular problem in West Virginia, and this study had collected thousands of samples along with patient histories and lots of clinical data. My belief is that the important thing about these samples is that the patients are young yet have serious periodontal issues. In any event, we convinced the people who actually run the project (based elsewhere) to let us have (a small part of) 8 samples out of all their thousands sitting in the freezers down the hall.

We did PCR with "universal" primers for the bacterial 16S rRNA gene, and cloned and sequenced the numbers you see in the table. It's not a big study (we don't have much money anyway), but we saw something which I think is truly significant. In high disease individuals, a broad group of microbes from the Clostridiales including an unusual clade called the Veillonellaceae are increased in abundance, whereas the sequences from control individuals in this clade were all very closely related to Veillonella parvula.

One reason this observation may be important is that the so-called "red complex", which is thought to be associated with serious periodontal disease, can only be recovered in about half the individuals with this diagnosis (not even considering abundance).

That story is in the modified version of the map above, where I drew a red box around the region of interest for the three controls, or "low disease" samples. Time went on, and another set of samples was added to the study from a different group, and we were able to get the work published. So that's why the study looks so old-fashioned, in an era of millions of reads, we've got about a thousand.

My role in all this was to actually do the analysis. I remember "we" wrote a grant (actually, someone else did!) and listed me as a technical expert in bioinformatics. Of course, the reviews were scathing. Dr. E doesn't have a degree in bioinfomatics. How could he know anything?

Well, I've learned a few things over the years. Rule one is, never make your own database: let someone else do it. That's why HOMD (and Greengenes and RDP) are so great. I particularly like the tools at the RDP site. It is very nice software.

And rule two is, if you live long enough, you will see work that took you months or years to accomplish be achieved using new tools in mere seconds or hours. Sequencing is a great example of this. When I was young I spent most of three months getting 500 bp; when I was a bit older I invested six months for 3.5 kb; still later it was a year for 20 kb.

This project is another example. I spent a year and more writing some 50 or so Python scripts (and rewriting them), and now QIIME does the whole thing in mere seconds.

Well, not quite the whole thing. I have a bit more to do with this project. I want to show you the UniFrac analysis of beta diversity, and show how to make what I think is a nicer plot of the PCoA results. Also, I want to show some phylogenetic trees detailing the increased diversity (species richness, really) in the Veillonellaceae that I mentioned.

And I should say: it's been fun. Even if I don't have that degree, or any papers with Rob Knight, I think I've learned something about Bioinformatics in the last 5 years.

Dental project (4)

This post is one of a series (see dental project here or in the sidebar).

After getting a set of sequences and removing chimeras, the next step is almost anticlimactic. We just copy a modified version of the shell script (from here, without the cd calls) or paste in the commands working from the dental directory (either individually, or all at once):

#!/bin/bash pick_otus.py -i seqs.fna -m uclust -s 0.97 -o otus pick_rep_set.py -f seqs.fna -i otus/seqs_otus.txt -m most_abundant -o otus/reps.txt align_seqs.py -i otus/reps.txt -m pynast -t ~/data/core.txt -o aln assign_taxonomy.py -i otus/reps.txt -m rdp -o tax filter_alignment.py -i aln/reps_aligned.txt -m ~/data/mask.txt -o aln2 make_phylogeny.py -i aln2/reps_aligned_pfiltered.fasta -o figs/tree.tre make_otu_table.py -i otus/seqs_otus.txt -t tax/reps_tax_assignments.txt -o figs/otu_table.txt summarize_taxa.py -i figs/otu_table.txt -o figs/otu_table_Level3.txt -L 3 plot_taxa_summary.py -i figs/otu_table_Level3.txt -l Phylum -o figs -k white make_otu_heatmap_html.py -i figs/otu_table.txt -o figs

It's all over in a few seconds.

The heatmap QIIME produced is at the top of the post. It is truly a remarkable html page, with a graphic where you can reorder the columns or rows by drag and drop, and redo the map at different threshholds for the OTUs, etc. I've never seen anything quite like it. But (and this is just me), it's not pretty enough.

So what I'd like to do from here is to show you how I currently make heatmaps with matplotlib, and we'll get into that next time.

First, I have to extract the data from QIIME. The script is complicated a bit by an additional job: I'm going to organize the rows and columns. (QIIME can do this too---see the tutorial).

The columns will be in the order they appear in sample_names.txt and the rows as they appear in genera_and_colors.txt. These files are in the same directory. The second one starts like this:

# Bacteria black Bacteria # Bacteroidetes green Bacteroidetes Bacteroidales Rikenella Prevotella Porphyromonas Saprospiraceae Sphingobacteriales

I just do this:

> python grab_qiime_data.py > data.csv

Here's the first part of data.csv:

,DB,DC,DM,DF,DL,DG,DI,DA,DT,DQ,DV,DAA,DZ Bacteria,,,,2,,,,,,,,,12 Bacteroidetes,,,,,,2,,,,,,, Bacteroidales,,,,,,1,,,,,,, Prevotella,,,,,3,2,,1,,,,, Porphyromonas,,,,,,1,,1,,,,, Sphingobacteriales,,,,,,,,,,1,,, Capnocytophaga,9,,6,2,3,6,1,19,31,35,2,8,7 Fusobacterium,1,,,4,,2,1,2,,2,1,3,1

The leading comma on line 1 is so the column headers line up properly in a spreadsheet. Speaking of spreadsheets, here is a screenshot after dropping the data.csv file onto Numbers (you could use Excel, of course):



That was painless! Zipped project files in Dropbox (here).

Dental project (3)

wikimedia

Before starting on analysis of the 1124 sequences from last time (here), we need to check for chimeras.

And at this point, I have a confession to make. It turns out there are 3 and perhaps 4 chimeras in the set of sequences from Genbank. I discovered this unwelcome fact a few weeks ago when playing with the QIIME toolkit. Since one of the pieces of software they recommend is ChimeraSlayer, I tried it out on these sequences.

Make a directory temp with a copy of seqs.fna. The sequences first need to be converted to NAST format, then we can run ChimeraSlayer.pl.

prog1=~/Software/microbiomeutil_2010-11-02/NAST-iEr/run_NAST-iEr.pl $prog1 --query_FASTA seqs.fna > seqs.nast prog2=~/Software/microbiomeutil_2010-11-02/ChimeraSlayer/ChimeraSlayer.pl $prog2 --query_NAST seqs.nast

It takes the better part of an hour on my slowest machine (a 5 year old iMac).

seqs.nast.CPS.CPC.wTaxons has flagged four sequences:

DA228 INTRA-GENUS DQ822 INTRA-PHYLUM DV55 INTRA-FAMILY DAA89 INTRA-FAMILY

I grab those four by hand into a new file suspects.fna (there is probably a better way) and do:

$prog1 --query_FASTA suspects.fna > suspects.nast $prog2 --query_NAST suspects.nast --printCSalignments option

The output shows there is definitely a problem. In suspects.nast.CPS.CPC.wTaxons we have:

ChimeraSlayer DQ_822 S000427388 S000260335 1.0566 98.74 100 0.7978 74.56 0 YES NAST:1861-1863 ECO:324-325 Streptococcus Streptococcus cristatus (T); NCTC12479; AB008313 Streptococcus cristatus Lachnospiraceae Incertae Sedis Clostridium aerotolerans (T); DSM 5434; X76163 Clostridium aerotolerans INTRA-PHYLUM Per_id parents: 73.80 Per_id(Q,A): 93.45 --------------------------------------------------- A: S000427388 99.64 78.63 ~~~~~~~~~~~~~~~~~~~~~~~~\ /~~~~~~~~~~~~~~~~~~~~~~~~ Q: DQ_822 DivR: 1.057 BS: 100.00 | Per_id(QLA,QRB): 98.74 | | (L-AB: 72.50) | (R-AB: 76.92) WinL:0-279 | WinR:280-396 | Per_id(QLB,QRA): 74.56 | DivR: 0.798 BS: 0.00 | ~~~~~~~~~~~~~~~~~~~~~~~~/ \~~~~~~~~~~~~~~~~~~~~~~~~~ Q: DQ_822 72.86 96.58 ---------------------------------------------------- B: S000260335 Per_id(Q,B): 79.85 DeltaL: 26.79 DeltaR: -17.95 ! CTAACGAGGAGGCGCTTGGTTAAGGGCTAGCTAAATTGCATGATGGTCAATGGGAAATCC A: S000427388 CTAACGAGGAGGCGCTTGGTTAAGGGCTAGCTAAATTGCATGATAGTCAATGGGAAATCC Q: DQ_822 TCCGACTAAGATTCGGAATCGGGCACAAGATCTCGACAGGCAGCACAGTGGAACTCCGGT B: S000260335 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!! ACCCTTGTCGTGACATC A: S000427388 ACCCTTGTCGTGACATC Q: DQ_822 GGTGCCTGGACCCAGCT B: S000260335 !!!!!!!!!!!!!!!!! ** Breakpoint ** !!!!!!!!!!!!! !!!!! !!!!!!! GACCGACGGTCGAGTGTATCGGGGTAAA A: S000427388 AGAGAGGACCTCGGTTATATGACAGCGG Q: DQ_822 AGAGAGGACCTCGAGTATATCACAGCAG B: S000260335 !! ! !

The first match is great for a while, then terrible, and the second is the converse.

I need to look into whether I should update the Genbank records, but I guess probably the answer is yes.

Anyway, I should have discovered this easily. I wrote a Python tool that looks for chimeras by BLAST of the front and back "halves" of each sequence against our local "boutique" database. It prints the top five hits for each. Here is the output for three of the suspects:

>DQ822 BLAST front(len = 207): 320 207/207 100.00 Streptococcus_clone_BP2-57_AB121930.1 319 205/207 99.03 Streptococcus_clone_502H08_AM420202.1 323 203/207 98.07 Streptococcus_cristatus_AB008313.1 321 200/207 96.62 Streptococcus_clone_BW009_AY005042.1 334 199/209 95.22 Streptococcus_sanguinis_SK36_SK36 BLAST back(len = 187): 349 179/187 95.72 Uncultured_clone_4.59_DQ346409.1 79 179/187 95.72 Clostridiales_clone_301C11_AM420062.1 121 178/187 95.19 Eubacterium_clone_DO008_AF385508.1 117 178/187 95.19 Eubacterium_clone_BP2-88_AB121960.1 115 178/187 95.19 Eubacterium_clone_BL026B96_AY806377.1 >DV55 BLAST front(len = 210): 353 207/210 98.57 Uncultured_clone_E105_DQ326659.1 100 206/210 98.10 Dialister_sp._E2_20_AF481209.1 94 205/210 97.62 Dialister_invisus_AY162469.1 31 190/207 91.79 Allisonella_clone_BL34_DQ130020.1 93 192/210 91.43 Dialister_clone_MCE7_134_AF481210.1 BLAST back(len = 190): 373 190/190 100.00 Veillonella_parvula_X84005.1 372 190/190 100.00 Veillonella_clone_X042_AF287781.1 370 190/190 100.00 Veillonella_clone_BU083_AF366266.1 369 190/190 100.00 Veillonella_clone_AA050_AF287782.1 371 182/183 99.45 Veillonella_clone_R1_DQ123569.1 >DAA89 BLAST front(len = 210): 283 209/210 99.52 Selenomonas_clone_CI002_AF287798.1 286 202/210 96.19 Selenomonas_clone_EQ054_AF385495.1 288 201/210 95.71 Selenomonas_clone_FT050_AY349403.1 298 189/195 96.92 Selenomonas_noxia_AF287799.1 297 200/210 95.24 Selenomonas_infelix_AF287802.1 BLAST back(len = 190): 373 184/187 98.40 Veillonella_parvula_X84005.1 370 184/187 98.40 Veillonella_clone_BU083_AF366266.1 372 181/184 98.37 Veillonella_clone_X042_AF287781.1 369 181/184 98.37 Veillonella_clone_AA050_AF287782.1 371 182/187 97.33 Veillonella_clone_R1_DQ123569.1

Note on sequence titles: I just introduced the underscore recently (as in DA_228), so this output doesn't have them.

It's pretty obvious that these guys are problematic. What happened is that I integrated the tool into the toolchain, but I never wrote code to look through the output and flag potential problems. I always did it manually, and as additional sequence samples were added to the experiment, I forgot to carry out this step.

Moral of the story: if you want to be sure something gets done, every time, you need to automate it completely! Otherwise you might forget.

We'll remove these from our sequence file by hand. Now there are 1120.

DA_228 DQ_822 DV_55 DAA_89

Dental project (2)

This is a series, first post here. Before we do anything else, we need to clean up the titles on the sequences. They come from Genbank like this:

>gi|324104022|gb|HQ894465.1| Uncultured bacterium clone DA19 16S ribosomal RNA gene, partial sequence

We want this:

>DA_19

.

We could do something like a regular expression, but that's overkill. Note that the alpha part is variable length, so we have to be a little bit smart. But these are so regular, it's easy. Also, UniFrac wants an underscore, so we add that.

Just do this from the command line:

> python retitle.py > seqs.txt

If you count the sequences, you should have 1124.

retitle.py

import utils as ut digits = '0123456789' data = ut.load_data('results.txt') data = data.strip().split('\n\n') for item in data: title,seq = ut.clean_fasta(item) e = title.split()[4] d = ''.join([c for c in e if c in digits]) s = ''.join([c for c in e if not c in digits]) print '>' + s + '_' + d print seq print

Dental project (1)

I'd like to spend a few posts talking about a project we just completed on surveying the bacterial species present in patients with gingivitis compared to normal controls. The main reason is to exercise QIIME with a different data set, but I'd also like to say a bit about the project.

I checked this morning at PubMed and found the paper has come out:

Olson 2011 PMID 21362199

Click on the link to download, 7.6 MB, that's just about the largest file size for a paper that I ever saw. Opening it up, I see why. It's still in manuscript form, and some of the figures are quite big. To play with this, we'll need to get the sequences. Luckily they were just posted by Genbank the other day.

I wrote a script to grab the sequences in chunks of 40, with a timer to sleep for 10 seconds between requests. The first sign of trouble was here:

HQ895465 HQ895504 1 URL Error HQ895505 HQ895544 1 ..

but eventually, we did another request for this batch which looked like it worked:

HQ895465 HQ895504 2

but actually, the file contains this near the end:

PubseqAccess cmd(id_gi_by_word HQ895505) failed with Cannot get server name from load balancer 'PUBSEQ_OS_PUBLIC' : errmsg='Service not found' HQ895505

and then more of the same. Looking at the sequences, it seems they cut us off with 1000 sequences.. stopping with HQ895464.1

I thought this should be OK. It's very early in the morning, with more than 3 seconds between requests, but apparently we ran up against some kind of limit.

I give the server some time to calm down, (and change the name of the file we've written to), edit the list and try again:

> python fetchSeqs.py HQ895465 HQ895504 1 HQ895505 HQ895544 1 HQ895545 HQ895584 1 HQ895585 HQ895588 1

then combine by hand.. Next time we'll take a look at them.

fetchSeqs.py

import urllib2, sys, time from utils import load_data ncbi = 'http://eutils.ncbi.nlm.nih.gov' eutils = ncbi + '/entrez/eutils/' efetch = 'efetch.fcgi?' def chunks(L,SZ): rL = list() while L: rL.append(L[:SZ]) L = L[SZ:] return rL def fetch(L): s = eutils + efetch s += 'id=' + L s += '&db=nucleotide&rettype=fasta&retmode=text' try: FH = urllib2.urlopen(s) data = FH.read() except: raise ValueError('URL Error') if 'NCBI C++' in data: raise ValueError('Empty') elif not data: raise ValueError('Empty') return data def run(FH): first = 894465 #first = 895465 last = 895588 L = ['HQ' + str(n) for n in range(first,last+1)] rL = list() L = chunks(L,SZ=40) L = [(e,1) for e in L] while L: sL,n = L.pop(0) print sL[0], sL[-1], n if n > 3: continue try: s = fetch(','.join(sL)) FH.write(s) except ValueError as e: print e L.append((sL,n+1)) time.sleep(10) if __name__ == '__main__': FH = open('results.txt','w') L = run(FH) FH.close()

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).