Incorprorating genomic 16S gene copy number in analysis

[Jorrit-Jan Hofstra]

Hi all,

I was wondering if anyone has a clear idea about the need for incorporating 16S gene copy number information when analizing samples. The number of genomic 16S copies varies greatly from bacteria to bacteria. Some species are reported to have >10 16S copies in their genome. Not taking this into account could lead to significant errors in reported relative abundances. Knowledge about exact 16S copy numbers may not be perfect so correcting for them may prove a challenge. Just ignoring differences in genomic 16S copy number differences does not seem to be the right approach to me, but as I'm relatively new to the field I would like to hear what others think about this. (reference: Steven W. Kempel and Martin Wu et al, PLOS Computational Biology October 2012)

Is there a way for QIIME to perform a correction for 16S copy number? If not, should such a step be introduced?

Best regards,
Jorrit

[Tony Walters]

Hello Jorrit,

Your concerns are legitimate, although not particularly easy to address. One limitation is the number of sequences that could be analyzed, as we would be limited to those with the whole genome sequenced (or otherwise had the SSU copy number determined). With these data on hand, we could take a closed-reference picked dataset (against the taxa with known SSU copy numbers), do a correction for counts based upon copy number, and convert this back to a biom format OTU table. I think this would work well for mock community data, but one would lose a lot of sequences from his or her data even with data that have a lot of genomes sequenced (e.g. human-associated microbes).

This functionality isn't currently in QIIME though, so it would have to be done by hand.

There are a number of other issues to consider if one wanted to improve the accuracy of the relative abundance of taxa: sampling bias towards certain taxa, DNA extraction bias, PCR bias (shorter amplicons favored over longer ones), primer bias (not just better matches to certain taxa with your primer set, but roughly equimolar degenerate primers with one sequence binding higher abundance taxa X with a perfect match and the other sequence binding lower abundance taxa Y with a perfect match would likely inflate the apparent abundance of taxa Y without strict limits on cycle number), bias towards shorter amplicons on sequencing platforms, and bias against certain taxa (e.g. those with homopolymers on 454) on the sequencing platform.

It may be better in the end to narrow down the particular taxa of interest in a given study and try to quantify the abundance with an alternative method (e.g. qPCR) rather than trying to correct for the various biases present in amplifying a community of microbes.

-Tony

--
 
---
You received this message because you are subscribed to the Google Groups "Qiime Forum" group.
To unsubscribe from this group and stop receiving emails from it, send an email to qiime-forum...@googlegroups.com.
For more options, visit https://groups.google.com/groups/opt_out.
 
 

[Jesse Zaneveld]

Hi Jorrit,

I have a bit of information to add to Tony's comments.   As you say, simply ignoring copy number is unsatisfying, but is also where the field has been until very recently.   If you would like to try correcting for 16S rRNA copy number, it is now possible to do so with a bit of extra work by using our PICRUSt software package (http://picrust.github.io/picrust/).    I'll talk a little more about how this works below.

HOWEVER, if you are new to analyzing 16S rRNA datasets, I would strongly recommend running through your data using standard workflows, and only then branching out to add in extra analyses like this one.  This has the advantage that you don't need to spend a long time worrying about whether the correction will affect your data before diving in.  Then you can just test it empirically by first running the analysis in the standard way, and following up with a second run that corrects for copy number, if desired.   I mention this only because I have seen some colleagues get bogged down in endless cycles of quality control before even looking at the biological results (that often may not change that much based on the fine details of the analysis).

That being said, Tony's e-mail accurately describes the current state-of-the-art. Although 16S rRNA copy number could, in theory, have significant effects on conclusions about beta diversity, it is not generally corrected for, mostly because of the issues (how do you do it given incomplete reference genome coverage?) that Tony described.  

However, recently it has become possible to get a reasonable estimate of 16S rRNA copy number even given incomplete availability of reference genomes.   PICRUSt takes a very similar approach to what Steve Kemble (along with Martin Wu, Jessica Green and Jonathan Eisen) lay out in the paper that you sent.   For us, predicting 16S rRNA copy number was an intermediate normalization step to predicting full metagenomes from 16S + sequenced genomes.   However, we have compared the accuracy of our 16S copy number predictions by leave-one-out cross-validation on all sequenced genomes, and we get accuracies that are very similar to the values published by Kemble et al.   Both PICRUSt and the Kemble et al approach are able to infer copy numbers by using ancestral state reconstruction to extend observed copy numbers to internal nodes.  So for any organism known only by 16S and not by complete genome, the 16S copy number prediction is the copy number reconstructed to be present at the last common ancestor of that tip and a sequenced genome.  Then the model of evolution is used (optionally) to build 95% confidence intervals around the prediction.  So of course the method is more accurate in groups that have more sequenced genomes and less accurate in groups that have fewer.   

So we now have some capability of normalizing for 16S rRNA copy number in practice.    To do so in practice you'll just need to pick OTUs against a special reference set (which has greengenes + IMG genomes) in QIIME, and PICRUSt can handle the rest using the same BIOM files as QIIME.   See the tutorial on the web page (above)

Since this capability is new, there has been very little rigorous research so far into how much of an effect this normalization really has, in practice, on conclusions about beta-diversity.

In the Kemble et al paper they showed that it might only matter in some datasets, and even then their analysis of clustering by beta diversity didn't test whether changes were significant, so we really don't know yet whether this normalization is important or can safely be ignored.  

Cheers,
Jesse

p.s. If you run into any issues running picrust we also have a user forum you can ask questions in: picrus...@googlegroups.com

[Jorrit-Jan Hofstra]

Hi Tony and Jesse,

Thank you very much! You are right. There are a number of issues to take into consideration if you want to improve the accuracy of relative abundances of the different taxa. Correcting for genomic 16S copy numbers seems easier said than done. I will definitely look into the suggested methods to attempt correcting for this to see if it affects conclusions on alpha and beta diversity outcomes.

Thanks!
All the best
Jorrit

[Doruk]

Hi Jesse, 

I am aware of the normalize_by_copy_number.py script in picrust, but the function accepts a .biom file for picked otus. Since I have already mapped my reads to the best mapping reference sequences (no otu picking), I have the counts how many times a reference 16S sequence is hit, (90,000 of them), and I would like to normalize them by their estimated copy numbers. This study http://rrndb.umms.med.umich.edu/estimate.php allows such a feature, but it forces the user to do the task in batches of small sequences. Given the number of sequences I need to use (90,000), it looks a bit unfeasible. What would be a good way to estimate the CNVs of individual sequences through picrust? Would picking OTUs (closed reference) for the mentioned 90,000 sequences (say I have a single sample, and 90,000 sequences inside, each occuring once), and then using the normalize_by_copy_number.py work? Since we know the abundance of each otus before normalization, I thought we can deduce the copy numbers by looking at the normalized otu abundances for each of the 90,000 sequences, so long as they happen to be in a OTU picked against greengenes. 

Another question, since those 90,000 sequences are from RDP, is there a way to reach those copy numbers by picking otus against RDP, to guarantee that I have every single 90,000 sequences assigned a copy number estimate?

Thanks, 

Doruk

[Michael Ricketts]

I am greatly interested in this topic as well... Was there ever an answer given to this??  I mean using PICRUSt's normalize OTU script on closed_reference biom tables?

Mike