metabuli classifiedRefiner issues

[MicroBTM]

In your documentation page, you have :

metabuli classifiedRefiner <i:read-by-read classification> <i:DBDIR> [options]

- read-by-read classification : The JobID_classifications.tsv file generated by the `classify` step.
- DBDIR : The same DBDIR used in the `classify` step.

, but I think you meant for the second input to be <i:TAXDUMP>, as in the directory with taxonomy dump files that was used to create the database.

When I attempt to use metabuli classifiedRefiner with the taxdump files for my custom database, it appears to finish, but it always thinks there's an existing file with the same name that it's trying to write to (when there isn't).

Metabuli Version (commit):                                                                      1.1.1
Remove unclassified reads                                                                       false
Exclude taxId as well as its children
Select taxId as well as its children
Select columns with number, (7:full lineage, generated if absent)
Make report of refined classification file                                                      false
Adjust classification to the specified rank                                                     species
0: without higher rank, 1: with higher rank, 2: separate file for higher rank classification    0
Threads                                                                                         24
Min. sequence similarity score                                                                  0

Loading nodes file ... Done, got 2729221 nodes
Loading merged file ... Done, added 96668 merged nodes.
Loading names file ... Done
Init computeSparseTable ...Done
Write refined classification result to:
results/btk_all_classifications_refined.tsv
results/btk_all_classifications_refined.tsv is already exists.
Time for processing: 0h 0m 3s 287ms

This blocks any refined output data from saving.

Note: I'm using a conda install on WSL-2

[jaebeom-kim]

Hi. Thank you for sharing issues. I think I solved the issue previously. Will be included in the next release.
Please try to use the latest pre-built binary.
#161

[MicroBTM]

Thanks for the speedy reply! I'm looking forward to testing the next release.

Another issue that I've encountered is an inability to detect E. coli phage phiX174 (GCF_000819615.1) with multiple databases that include it, like your pre-built refseq_prokaryote_virus database. Do you know why this phage might not be detectable? Even if you feed the phage's genomic FASTA file as input and treat it like long reads with --seq-mode 3, metabuli won't classify it.

[jaebeom-kim]

I queried the phage against the refseq_prokaryote_virus database.
After analyzing the internal k-mer match distribution, I found that the query matches not only to the virus itself but also to multiple bacterial species with very similar mapping qualities; therefore, their LCA is calculated as root (taxID 1).

Staphylococcus cohnii
Staphylococcus haemolyticus
Staphylococcus simulans
Sinsheimervirus phiX174
Cupriavidus taiwanensis
Cupriavidus neocaledonicus
Sphingorhabdus sp. Alg231-15
Erythrobacter sp. Alg231-14
Anderseniella sp. Alg231-50

This occurs because phage genomes are often integrated into bacterial genomes (prophages) or due to potential database contamination. It is an interesting observation and a problem that should be resolved, though it will take some time to address.

[MicroBTM]

That's such an odd/interesting hit list, because PhiX174 is strictly lytic without any integration machinery, and it's unable to infect gram-positives like Staphylococcus.

Last follow up question:
I’m curious whether you’ve explored sketch-based approaches like FracMinHash for classification (similar to what Sylph and Sourmash use). Do you think that approach could fit into Metabuli, or are there algorithmic reasons LCA works better for the design of the tool?

[jaebeom-kim]

Here is NCBI BLAST between the phage and Staphylococcus genus.
Maybe genome contamination.

[jaebeom-kim]

Thank you for interesting question.
Sketch-based approaches are suitable for containment or composition analysis, but not for classifying each short read.
For example, Sylph doesn't generate classification result for each read.
Sketching might be used for long-read classification.

Sketching-based methods tend to miss diverged species (which are not represented in database) and low-coverage species.
One of the goal of Metabuli is detecting diverged species or low-depth species.
For example, Metabuli can detect SARS-CoV-2 reads in COVID19 patient's RNA-seq data (very small viral reads compared to human reads) using SARS-CoV-1 genome in DB (SARS-CoV-2 is not represented in DB).

[MicroBTM]

Okay, yes, those are very clear cases of contamination based on the blast results. Thanks for sleuthing!

And thanks for the clarification on different approaches. I was asking because of your last point, specifically. Your nucleotide/protein k-mer strategy might extend the capabilities of sketching approaches into detection of diverged species. If so, it would be very cool to have both analysis options (for rapid or slow classification based on the user's needs) ... although, I'm sure that's much easier said than done lol

[MicroBTM]

Hi again. I've realized a more realistic/direct approach to fix the issue of failed virus classification for all cases (pro-phages, pro-viruses, or contamination in the reference sequences) could be to add an option which forces any read mapping to the "Viruses" node (taxid 10239) to forfeit all taxids which aren't connected back to that node.

[jaebeom-kim]

That's a wonderful idea. I will add it soon.

[jaebeom-kim]

I added --priority-taxid option in classify from v1.2.0. Please try it :)