Few reads are classified by Medi during kraken2 step

[richrr]

Hi

Thank you for this tool. I am excited to use this tool.

When I'm running human fecal shotgun metagenomics data through MEDI, I noticed that majority of the reads are being unclassified.

(medi) [rr@cn0116 tmpRR]$ kraken2 --db MEDI/medi/data/medi_db --confidence 0.3 --threads 24 --gzip-compressed --paired --output MF0035p.k2 --report MF0035p.tsv --memory-mapping test/work/e9/b9096a31dff7539598f895006c94b7/MF0035p_filtered_R1.fastq.gz test/work/e9/b9096a31dff7539598f895006c94b7/MF0035p_filtered_R2.fastq.gz
Loading database information... done.
21380521 sequences (5942.19 Mbp) processed in 60.246s (21293.4 Kseq/m, 5917.98 Mbp/m).
913 sequences classified (0.00%)
21379608 sequences unclassified (100.00%)

The following step returns an empty file and then the pipeline fails due to empty files being provided to the next commands.
architeuthis mapping architeuthis mapping filter MF0073p.k2 --data-dir MEDI/medi/data/medi_db/taxonomy --min-consistency 0.95 --max-entropy 0.1 --max-multiplicity 4 --out MF0073p_filtered.k2

Could you please help me figure out why the number of classified sequences are so low?

Is this an issue with my medi_db? It seems to have completed correctly during installation and the db building and hashing, but are there any diagnostics I can run to test it?

Do you have any sample paired-end fastq files and the resulting tables that you could provide me to test? Even if the results don't exactly match with my version of the database, it would help me troubleshoot if the pipeline is installed correctly and/or the expected results.

Thanks.

-Rich

[cdiener]

Hi Rich,

Yes that looks odd indeed. I would start by checking the size of the hash.k2d in the medi_db folder. That should be around 380GB or more.

drwxrwx--- 4 o_diener o_diener  499 Dec  2 15:50 .
drwxrwx--- 6 o_diener o_diener  191 Oct  7  2024 ..
-rwxrwxr-- 1 o_diener o_diener 4,1M Mar 25  2024 database100mers.kmer_distrib
-rwxrwxr-- 1 o_diener o_diener 1,1G Mar 25  2024 database100mers.kraken
-rwxrwxr-- 1 o_diener o_diener 3,8M Mar 25  2024 database150mers.kmer_distrib
-rwxrwxr-- 1 o_diener o_diener 994M Mar 25  2024 database150mers.kraken
-rwxrwxr-- 1 o_diener o_diener 221G Mar 24  2024 database.kraken
-rwxrwxr-- 1 o_diener o_diener 137M Nov 28 13:10 food_contents.csv.gz
-rwxrwxr-- 1 o_diener o_diener  81K Mar 26  2024 food_matches.csv
-rwxrwxr-- 1 o_diener o_diener 412G Mar 22  2024 hash.k2d
drwxrwx--- 9 o_diener o_diener  174 Mar 19  2024 library
-rwxrwxr-- 1 o_diener o_diener 117K Oct  8  2024 manifest.csv
-rwxrwxr-- 1 o_diener o_diener   64 Mar 22  2024 opts.k2d
-rwxrwxr-- 1 o_diener o_diener 733M Mar 19  2024 seqid2taxid.map
-rwxrwxr-- 1 o_diener o_diener 3,9M Mar 19  2024 taxo.k2d
drwxrwx--- 2 o_diener o_diener  527 Mar 19  2024 taxonomy
-rwxrwxr-- 1 o_diener o_diener  656 Mar 19  2024 unmapped.txt

If it is much smaller something might have gone wrong. The next thing to check would be the library folder there and verify that it includes all the decoys and the added folder.

drwxrwx--- 2 o_diener o_diener 48K Mar 19  2024 added
drwxrwx--- 2 o_diener o_diener 165 Mar 18  2024 archaea
drwxrwx--- 2 o_diener o_diener 165 Mar 19  2024 bacteria
drwxrwx--- 2 o_diener o_diener 129 Mar 18  2024 human
drwxrwx--- 2 o_diener o_diener 127 Mar 18  2024 plasmid
drwxrwx--- 2 o_diener o_diener 126 Mar 18  2024 UniVec_Core
drwxrwx--- 2 o_diener o_diener 165 Mar 18  2024 viral

Sorry for the issues, I'm already working on making the database build more resilient and improve the resume/caching support. If you e-mail me at mail [at] dienerlab.com with a Gmail address I can also share the built database with you so you can compare.

[richrr]

So my db definitely did not build correctly, it is much smaller. There might be issues within my library folder too. I've emailed you the list of the files inside the library folder.

[rodriguesrr@cn4261 medi]$ ll data/medi_db/
total 5.6G
drwxr-s--- 2 rodriguesrr MGC_CCRSF 4.0K Apr 10 15:12 ../
drwxr-s--- 2 rodriguesrr MGC_CCRSF 4.0K Apr 10 20:06 library/
drwxr-s--- 2 rodriguesrr MGC_CCRSF 4.0K Apr 10 20:28 taxonomy/
-rw-r----- 1 rodriguesrr MGC_CCRSF 9.5K Apr 10 20:31 unmapped.txt
-rw-r----- 1 rodriguesrr MGC_CCRSF 1.2G Apr 10 20:31 seqid2taxid.map
-rw-r----- 1 rodriguesrr MGC_CCRSF 4.2M Apr 10 20:32 taxo.k2d
-rw-r----- 1 rodriguesrr MGC_CCRSF 4.3G Apr 10 20:40 hash.k2d
-rw-r----- 1 rodriguesrr MGC_CCRSF   64 Apr 10 20:40 opts.k2d
-rw-r----- 1 rodriguesrr MGC_CCRSF  59K Apr 10 20:45 database.kraken
-rw-r----- 1 rodriguesrr MGC_CCRSF  36K Apr 10 20:45 database100mers.kraken
-rw-r----- 1 rodriguesrr MGC_CCRSF  121 Apr 10 20:45 database100mers.kmer_distrib
-rw-r----- 1 rodriguesrr MGC_CCRSF  36K Apr 10 20:45 database150mers.kraken
-rw-r----- 1 rodriguesrr MGC_CCRSF  121 Apr 10 20:45 database150mers.kmer_distrib
-rw-r----- 1 rodriguesrr MGC_CCRSF  77K Apr 10 20:45 food_matches.csv
-rw-r----- 1 rodriguesrr MGC_CCRSF 137M Apr 10 20:45 food_contents.csv.gz
drwxr-s--- 2 rodriguesrr MGC_CCRSF 4.0K Apr 10 20:45 ./
-rw-r----- 1 rodriguesrr MGC_CCRSF 106K Apr 10 20:45 manifest.csv

[rodriguesrr@cn4261 medi]$ ll data/medi_db/library/
total 0
drwxr-s--- 2 rodriguesrr MGC_CCRSF 4.0K Apr 10 15:12 UniVec_Core/
drwxr-s--- 2 rodriguesrr MGC_CCRSF 4.0K Apr 10 15:13 human/
drwxr-s--- 2 rodriguesrr MGC_CCRSF 4.0K Apr 10 15:14 plasmid/
drwxr-s--- 2 rodriguesrr MGC_CCRSF 4.0K Apr 10 15:14 archaea/
drwxr-s--- 2 rodriguesrr MGC_CCRSF 4.0K Apr 10 15:38 viral/
drwxr-s--- 2 rodriguesrr MGC_CCRSF 4.0K Apr 10 20:06 bacteria/
drwxr-s--- 2 rodriguesrr MGC_CCRSF 4.0K Apr 10 20:06 ./
drwxr-s--- 2 rodriguesrr MGC_CCRSF 4.0K Apr 10 20:28 added/
drwxr-s--- 2 rodriguesrr MGC_CCRSF 4.0K Apr 10 20:45 ../

Regarding, "making the database build more resilient and improve the resume/caching support", it seems that this fix/edit is still required for Downloading plasmid files from FTP. Not sure if this could be done within your installation files itself or would have to be initiated in the kraken2 code.

Thank you for your time and help.

Best,
Rich

[cdiener]

Hi Rich, yeah it failed in the Kraken2 default database download. Unfortunately, if that fails, the download has to be restarted from scratch (see the warning in the docs).

I also got your E-mail and will share the DB with you.

The good news is that the current development version is much more resilient and can resume interrupted downloads perfectly fine, but it is still not quite ready.

[nicodcroig]

Hello, I am also having trouble building the database, as I dont have too much resources. Is it possible that you send me the built database too?

Thank you in advance!

[cdiener]

Sure, just follow the instructions from the docs and send me a Gmail address to share it with.