Microbial diversity analysis using whole-genome metagenomic data

Question

[Antonio Camargo]

I have data, obtained from a single metagenomic DNA sample, that consists of two MiSeq FASTQ files (R1 and R2) that I merged using PEAR.

Now I want to estimate the abundances of the bacteria taxa to generate a figure like this one:

Figure from: Panosyan, Hovik, and Nils‐Kåre Birkeland. "Microbial diversity in an Armenian geothermal spring assessed by molecular and culture‐based methods." Journal of basic microbiology 54.11 (2014): 1240-1250.

The problem is that there wasn't a step of amplification of the 16S region as the goal of the sequencing was to discover new genes. I've already isolated 16S reads from my sample using SortMeRNA, but it seems like softwares that do OTU picking, taxonomic assignment and diversity analyses (such as mothur and QIIME) require that all the reads come from the same region of the 16S gene.

Is there a way of using these 16S reads that I've filtered using SortMeRNA in a diversity analysis using mothur/QIIME?

[Antonio Camargo]

Cross-posted on StackExchange

Answer 1

[Bioinformatics_NewComer]

If you've access to computing cluster then metagenomic tools like CLARK-S would be good to try. They give you abundances and allow you to perform other analyses.

[Antonio Camargo]

Thank you. For what I've read, tools like CLARK, KRAKEN and Kaiju are the answer for my problem.

Answer 2

[Brian Bushnell]

You could try assembling your 16S reads with an assembler that deals well with branches (possibly SPades), then aligning the resulting assemblies to other 16S sequences and trimming off the bases that go off the end (and are thus not 16S bases). I'm not sure how well that would work; depends on the data.

But I think your best bet would be to use the shotgun data as shotgun data, instead of trying to shoehorn it into 16S-based tools. We commonly assemble the whole metagenome, map the reads to the assembly to calculate coverage, and then align the contigs to existing databases like RefSeq to find out what they are. Once you know that contig_123 maps to E.coli, and has coverage of 43x, you can say you probably have 43x coverage of E.coli in your data. Whether this approach works depends on whether you have enough data to assemble; if only, say, 10% of your reads map to the assembly, then it's pretty much a failure and you'll need a different method.

One thing to try in that case is to compare reads directly to RefSeq to find what organisms they came from. You can get a list of organisms observed in your data with BBMap like this:

sendsketch.sh in=data.fq refseq records=400

Once you know which organisms are present, you can download their genomes and map reads to them for quantification purposes. Mapping to all of refseq directly normally takes too much time or memory to be practical.

You might also check out KRAKEN which looks like it is designed for this purpose. I have not tried it, though.

[seta]

Hi all,

Long delay, but may be helpful. I was wondering how I can study the microbiome (say microbial flora) from whole genome sequencing data of healthy people belong to a given population? Please kindly let me know your suggested procedure, any pipelines?

Thank you

[Bioinformatics_NewComer]

Please create a new thread. Thanks.