2
$\begingroup$

I built a protozoa kraken 2 database ("db1"), from which I further built two other databases ("db2" and "db3") by adding extra genomes in the Leishmania sub-genus*.

When comparing results obtained when analyzing sequencing data using these databases, I noticed a strange increase in the number of reads assigned to the Leishmania chagasi species (which is a terminal node in the taxonomy for both databases) when going from db2 to db3.

According to my superficial understanding of kraken, I didn't expect this to happen, because db3 was built using a higher diversity of genomes outside the L. chagasi species than db2. If I expected any effect, this should have reduced the number of k-mers that are found specifically in L. chagasi, because some of these k-mers in db2 could have been present in the new genomes added to db3.

This increase in analysis reports parallels an increase in the number of minimizers as reported by kraken2-inspect. Below is a table comparing the kraken2-inspect reports for db2 and db3.

I only show the lines corresponding to Leishmania subgenus and its parents.

Compared to the kraken2-inspect output, I removed the taxid and re-ordered columns, for readability. nb_min and nb_min_lvl correspond to the number of minimizers reported on columns 2 and 3 of the kraken2-inspect default output (See https://github.com/DerrickWood/kraken2/blob/master/docs/MANUAL.markdown#inspecting-a-kraken-2-databases-contents).

The line for L. chagasi is in bold, and shows a substantial increase (more than 20-fold) from 63874 minimizers in db2 to 1389504 in db3:

taxname lvl % (db2) nb_min (db2) nb_min_lvl (db2) % (db3) nb_min (db3) nb_min_lvl (db3)
root R 100.0 228330589 0 100 254009353 0
cellular organisms R1 100.0 228330589 0 100 254009353 0
Eukaryota D 100.0 228330589 4685 100 254009353 6669
Discoba D1 23.4 53431766 0 31.14 79110530 0
Euglenozoa P 23.4 53431766 0 31.14 79110530 0
Kinetoplastea C 23.4 53431766 0 31.14 79110530 0
Metakinetoplastina C1 23.4 53431766 0 31.14 79110530 0
Trypanosomatida O 23.4 53431766 0 31.14 79110530 0
Trypanosomatidae F 23.4 53431766 1645 31.14 79110530 2029
Leishmaniinae F1 19.87 45358915 0 27.97 71037679 0
Leishmania G 19.87 45358915 259764 27.97 71037679 305088
Leishmania G1 14.67 33503508 4587885 23.3 59182271 12771940
Leishmania donovani species complex G2 3.61 8243592 6215217 5.01 12729716 1378715
Leishmania donovani S 0.63 1439307 1439307 3.71 9428437 8055955
Leishmania donovani Ld 2001 S1 0.48 1207432 1207432
Leishmania donovani Ld 39 S1 0.06 165050 165050
Leishmania chagasi S 0.03 63874 63874 0.55 1389504 1389504
Leishmania infantum S 0.23 525194 520206 0.21 533060 528882
Leishmania infantum JPCM5 S1 0.0 4988 4988 0 4178 4178
Leishmania mexicana species complex G2 4.39 10033887 4875993 4.18 10626363 5076482
Leishmania amazonensis S 1.22 2782684 2782684 1.23 3123079 3123079
Leishmania mexicana S 1.04 2375210 0 0.96 2426802 2191283
Leishmania mexicana MHOM/GT/2001/U1103 S1 1.04 2375210 2375210 0.09 235519 235519
Leishmania martiniquensis S 3.44 8749402 8749402
Leishmania major species complex G2 2.49 5691012 0 2.35 5980712 0
Leishmania major S 2.49 5691012 5634913 2.35 5980712 5350171
Leishmania major strain LV39c5 S1 0.17 419323 419323
Leishmania major strain SD 75.1 S1 0.08 192440 192440
Leishmania major strain Friedlin S1 0.02 56099 56099 0.01 18778 18778
Leishmania tropica species complex G2 1.89 4803669 0
Leishmania tropica S 1.89 4803669 4796954
Leishmania tropica L590 S1 0 6715 6715
Leishmania aethiopica species complex G2 2.17 4947132 0 1.39 3520469 0
Leishmania aethiopica S 2.17 4947132 0 1.39 3520469 3057668
Leishmania aethiopica L147 S1 2.17 4947132 4947132 0.18 462801 462801
Viannia G1 5.08 11595643 5588978 4.55 11550320 5553385
Leishmania braziliensis species complex G2 1.33 3037436 0 1.19 3032336 0
Leishmania braziliensis S 1.33 3037436 0 1.19 3032336 0
Leishmania braziliensis MHOM/BR/75/M2904 S1 1.33 3037436 3037436 1.19 3032336 3032336
Leishmania guyanensis species complex G2 1.3 2969229 0 1.17 2964599 0
Leishmania panamensis S 1.3 2969229 2969229 1.17 2964599 2964599

My question: How is such an increase possible, when the only genomes added are outside the clade?

I don't know if this can help understanding what happens, but L. chagasi could actually be considered a sub-clade of L. infantum (chagasi is the descendant of infantum brought to the New world by European conquistadores), but in the NCBI taxonomy used to build the databases, they are sister S-level taxa within the G2 "L. donovani species complex" (together with the L. donovani species, that has two S1-level strains added in db3 compared to db2).

* Here is how I built the databases

kraken2-build --download-taxonomy --db protozoa
kraken2-build --download-library protozoa --db protozoa
cp -r protozoa db2
cp -r protozoa db3
mv protozoa db1
kraken2-build --build --threads 16 --db db1
# Download genomes
# Add genomes to db2 using `kraken2-build --add-to-library <genome>/*.fna --db db2`
kraken2-build --build --threads 16 --db db2
# Add even more genomes do db3 using `kraken2-build --add-to-library <genome>/*.fna --db db3`
kraken2-build --build --threads 16 --db db3
Improve this question
$\endgroup$
2
  • $\begingroup$ I will look at this question over the next week. There's clearly a lot going on here. What I need to think about is "outside the clade" ... I'm not sure about that. $\endgroup$
    – M__ ♦
    Commented Jul 28 at 2:31
  • $\begingroup$ To elaborate: If I had added genomes belonging to L. chagasi (for instance, at the S1 level), it would not have surprised me to see a larger number of k-mers associated to L. chagasi at the species level. New genomes should give a better idea of the diversity observable within a given clade. However, with genomes only added elsewhere, I don't get it. Maybe something to do with the distinction between k-mers and minimizers? $\endgroup$
    – bli
    Commented Jul 29 at 13:20

1 Answer 1

Reset to default
1
$\begingroup$

In one of its possible output modes, kraken2 will give you read-level and kmer-level information about how it came to its conclusions about taxonomic assignment. You've provided a lot of detail in your question about your process of generating the data, and the taxon assignment table, but I expect your specific question will best be answered by digging deeper into your own data.

To explain a bit further, individual reads are assigned to a taxon based on the kmer assignments for that read (i.e. each kmer within a read essentially votes for the read taxon). The behaviour of assignment after changing the database is unpredictable at a read level because the kmer assignments will change. Even if diversity is increased, that increased diversity will alter the "votes" from the kmers, and could result in a particular taxon crossing a threshold (and increasing in apparent abundance) because other taxa now have fewer votes.

In any case, Kraken2 is known to assign taxa in a way that is not necessarily a good representation of the actual content. If you want to make statements about the proportion of different taxa in your samples, I'd recommend that you try out using Bracken on your Kraken2 results. It will correct some assignments based on sequence similarities in the database:

Bracken is a companion program to Kraken 1, KrakenUniq, or Kraken 2 While Kraken classifies reads to multiple levels in the taxonomic tree, Bracken allows estimation of abundance at a single level using those classifications (e.g. Bracken can estimate abundance of species within a sample).

https://github.com/jenniferlu717/Bracken?tab=readme-ov-file#step-1-generate-the-bracken-database-file-databasexmerskmer_distrib

Improve this answer
$\endgroup$
2
  • $\begingroup$ Thanks for the ref to bracken. You wrote "I expect your specific question will best be answered by digging deeper into your own data.". Note that the puzzling observation I made is actually not dependent on my data, since it seems to stem from some effect at database construction. $\endgroup$
    – bli
    Commented Jul 27 at 21:01
  • 1
    $\begingroup$ @bli I've added a bit more detail to explain my reasoning. The main thing is that database changes lead to unpredictable effects. You can't always assume that class assignment fractions will go down when database diversity is increased. That's why it's important to use tools like Bracken that take this kmer weirdness into account. $\endgroup$
    – gringer ♦
    Commented Jul 28 at 1:11

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.