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