2
$\begingroup$

[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?

$\endgroup$

1 Answer 1

4
$\begingroup$

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:

Sankey plot showing metagenomic sequencing results. The three most abundant organisms are Escherichia coli, Shigella flexneri, and Salmonella enterica

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

$\endgroup$
5
  • $\begingroup$ Very nice, not possible to classify at strain resolution I guess? $\endgroup$ Feb 10, 2023 at 17:35
  • $\begingroup$ Not as far as I'm aware; you could try metamaps, but I have no experience with that. $\endgroup$
    – gringer
    Feb 10, 2023 at 19:46
  • 2
    $\begingroup$ If you have bracken data from multiple samples and you’d like to compare them with relative abundance stacked bar plots, I made a shiny app for that: github.com/acvill/bracken_plot $\endgroup$
    – acvill
    Feb 11, 2023 at 13:39
  • 2
    $\begingroup$ Also, I believe taxonomic level S1 represents sub-species (strain) classification. Here are some example bracken output files representing different taxonomic levels, includingS1. $\endgroup$
    – acvill
    Feb 11, 2023 at 13:45
  • $\begingroup$ Cool question and answer ... thanks. Also thanks to @acvil $\endgroup$
    – M__
    Feb 13, 2023 at 2:26

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.