Our group was hoping to use MEGAN as part of a long-read sample characterisation pipeline. In testing we have noticed some unexpected taxon classifications in cases where there are ambiguous protein matches with low read support: rather than going up higher in the taxonomy to find the lowest common ancestor that would cover the ambiguity, a specific assignment is made to one of the ambiguous taxa that does not seem to be a correct choice.
We are using MEGAN7 and the nr_full dataset, but can reproduce this issue with several custom databases of our own and/or using MEGAN6.
To save pasting a large BLASTX table from nr_full, here is a simple example using a DIAMOND database consisting of the entire RefSeq viral protein dataset. The sample was from a coronavirus OC43 infection and is publicly available. For the read in question, DIAMOND creates 3 protein matches:
| qseqid | sseqid | pident | length | mismatch | gapopen | qstart | qend | sstart | send | evalue | bitscore |
|---|---|---|---|---|---|---|---|---|---|---|---|
| ERR13982188.9357 | YP_009555238.1 | 90.4 | 591 | 32 | 11 | 1818 | 98 | 4383 | 4952 | 0 | 979 |
| ERR13982188.9357 | YP_009555260.1 | 90.4 | 591 | 32 | 11 | 1818 | 98 | 14 | 583 | 0 | 979 |
| ERR13982188.9357 | NP_150073.3 | 87.5 | 610 | 50 | 11 | 1871 | 98 | 4365 | 4952 | 0 | 973 |
| Taxon classifications for these DIAMOND hits are as follows: |
| MEGAN refID | MEGAN taxonID | MEGAN taxonName | taxonpicLineage |
|---|---|---|---|
| YP_009555238 | 31631 | Human coronavirus OC43 | Viruses;Riboviria;…;Betacoronavirus 1;Human coronavirus OC43 |
| YP_009555260 | 31631 | Human coronavirus OC43 | Viruses;Riboviria;…;Betacoronavirus 1;Human coronavirus OC43 |
| NP_150073 | 11128 | Bovine coronavirus | Viruses;Riboviria;…;Betacoronavirus 1;Bovine coronavirus |
Given the ambiguity between the 2 matched taxa and low support for both, we would expect this read to be classified as “Betacoronavirus 1” i.e. the LCA .
However, MEGAN7 daa-meganizer/daa2info classified the read as “Bovine coronavirus”.
Our question is whether this is indeed an unexpected result or are we misunderstanding the “longReads” LCA algorithm?
(For reference, our daa-meganizer options were: –longReads true –lcaAlgorithm longReads –lcaCoveragePercent 51.0 –minScore 50.0 –minReadLength 0.
We saw the same issue regardless of whether we used –ram alignedBases or –ram readCount.)
Overall, most reads do appear to be assigned to the expected taxa (from side-by-side comparison of MEGAN calls with the DIAMOND BLASTX output) but there always seem to be a handful of reads like the above whose classification does not appear consistent with the algorithm outlined in the manual.
We had hoped to use MEGAN to be capable of detecting cases where there is evidence from protein matches that we have a novel virus from the same family in our dataset: so specific assignments to lower ranks that seem zoonotic (for example) would function as a flag to us. However, due to the handful of misclassifications described above, we currently cannot, as in all cases we have examined the assignment should have been placed at a higher level of the taxonomy.
Many thanks in advance to the community for your help with this - any help/advice would be gratefully received!