Species abundance and rarefaction in metagenomics data

[OliverDeusch]

Hi,

I am doing a rarefaction analysis on a metagenomics dataset and ran into an issue. I will briefly describe the dataset before I ask the actual question.

For each sample I have approx. 50 million Illumina PE reads. After error correction I performed taxonomic and functional binning using MetaCV. I then summarised the results (applying a conservative cutoff) and obtained read counts for each species and finally sorted the species by abundance.

Looking at the summary about 2700 species were identified. The database contained ~3300 species. This feels like a lot of false positive species were called. This is supported by sorting the species by abundance. The 1200 most abundant species have 99.9 percent of the reads. The remaining 1500 species only 0.1 percent of the reads.

Clearly a species with 2 reads (out of millions) must be a false positive and should not be counted when estimating the number of species or doing a rarefaction analysis. But what cutoff does make sense? 10, 100, 1000 or even more reads?

Any input would be greatly appreciated.

Cheers,

Oliver

[JackieBadger]

You are right!
Far too many metagenomic studies include errors by classifying any sequence that has just a few reads as a "species". This is garbage. The base mismatch error rate of Illumina couple with ultra deep sequencing creates lots of erroneous variants.
Using an arbitrary global cut-off always seems dodgy to me. Plot a histogram of all your species depths... I'm guessing there will be a sharp drop at some point. Use this as your boundary between real things and artifacts. You can then demonstrate that the artifacts are descended from the real sequences using a network diagram. Id guess that the low copy number variants are very closely related to true sequences (base mismatch errors-which can occur systematically btw), or chimeras.