How do I measure the proportion of different bacteria in a sample from a high-throughput sequencer?

Question

[this question is based on a question that was asked on Reddit]

I have sequenced some [mostly cell-free] DNA from a sputum sample on a nanopore sequencing machine, and would like to know what is in it (i.e. the proportion of different bacteria). The run produced 2 million reads with an average read length of about 200bp.

What's the best way to go about doing a metagenomic classification of sequenced reads?

Answer 1

If you've got a fast solid state drive with >200G of space, then kraken2 + bracken works really well on both long and short reads (e.g. Nanopore, Illumina). It's also specifically designed for metagenomic analysis, unlike other general aligners like STAR, bwa, bowtie2, etc., so you get tree-based classification which can deal with partial / incomplete matches.

First, use kraken2 to map reads to taxa in memory-mapped mode (to reduce system memory consumption):

$ kraken2 --threads 10 --db /mnt/ultra_fast/kraken/metagenome --quick --memory-mapping reads/*.fastq.gz --report kraken2.report.txt > kraken2.out.txt

Then, use bracken to adjust the counts based on mapping / library probabilities:

$ ~/install/bracken/bracken -d /mnt/ultra_fast/kraken/metagenome -i kraken2.report.txt -o bracken.report.txt

I use pavian to look at the results, because it creates quite understandable Sankey plots:

But inspecting the _bracken_species.txt file is a pretty good text-based way to do it. If there's substantial amounts of a particular bacteria in the sample, it should come out in the results:

100.00  232939  0       R       1       root
99.99   232912  0       R1      131567    cellular organisms
68.52   159618  0       D       2           Bacteria
67.67   157635  0       P       1224          Proteobacteria
67.47   157152  0       C       1236            Gammaproteobacteria
61.82   144001  0       O       91347             Enterobacterales
61.79   143941  0       F       543                 Enterobacteriaceae
32.28   75183   0       G       561                   Escherichia
32.26   75147   75147   S       562                     Escherichia coli
0.02    35      35      S       1499973                 Escherichia marmotae
0.00    1       1       S       564                     Escherichia fergusonii

kraken2/bracken indexes can be downloaded from here (the build process takes ages; it makes more sense to use the pre-built indexes):

https://benlangmead.github.io/aws-indexes/k2

As far as I'm aware, kraken2 can't do strain-level assignment with the most commonly-used databases (although that may depend on what information that is stored in the database). For strain-level assignment, it seems like MetaMaps may work:

https://doi.org/10.1038/s41467-019-10934-2