QUAST v. 4.5 incorrectly assigning unaligned contigs

[BournSupremacy]

Hi there. I'm an MSc student using a pipeline developed by another group to identify contigs in human genomes that don't align to the reference genome. One step in the pipeline utilises QUAST v. 4.5 (the pipeline was made in 2019, I think) to identify the fully and partially unaligned reads using NUCmer, but I recently started doing some analysis of the final set of unaligned or "non-reference" sequences identified by QUAST and found that there was lots (and I mean like over 20 Mbp) of sequence that QUAST had identified as fully unaligned but that actually had quite strong alignment to the reference. Some as much as 90% identity and 100% coverage - which I know should have been labelled as aligned by QUAST. I've gone through the logs and found that NUCmer is correctly identifying the alignment, but somewhere along the way, the QUAST software appears to incorrectly analyse the NUCmer output and therefore incorrectly assign some aligned sequences as fully unaligned. I tested this again using QUAST v. 5.0.2 and Minimap2, and in this version the problem seems to be fixed and all unaligned contigs are correctly identified. This is where I saw that my amount of sequence decreased by 20 Mbp, meaning v. 4.5 had incorrectly assigned 20 Mbp of sequence as unaligned. I'm just wondering if this was a known issue in v. 4.5 or whether there is potentially something wrong with my setup of QUAST? It's not too important as I'm just going to adapt the pipeline to utilise the latest version of QUAST to avoid the error, but I will need to explain the change in my dissertation and I'd like to know if this problem had been seen and identified before. Thanks!

[alexeigurevich]

Hi!
We filter out alignments with identity below --min-identity (see here), the default cut-off is 95%. I assume that is why your 90% identity alignments were filtered out. Please check whether the pipeline explicitly sets this option or relies on the default value.

Still, if I remember correctly, the default cut-off is 95% since the very beginning (the option introduced in v.4.2; the default value was sometimes different, 90%, in the MetaQUAST mode, but it is not your case), so v.5.0.2 should also mark such contigs as "unaligned". If the real identity is close to 95%, than the switch from Nucmer to minimap2 may resulted in the dramatic change, e.g., the contig was 94% according to Nucmer and thus unaligned and 96% with minimap2 and thus aligned. There were some bugfixes regarding QUAST postprocessing of the raw alignments and calculation of their final identity value in v.5*; they are mostly minor and related to presence of long indels in the alignments.

In summary, v.5.* is a substantial update of v.4.5, so it could cause the difference. For nearly 100% IDY alignments there is no much difference, but for ~90% IDY the change of the core aligner inside QUAST could be really important.

[BournSupremacy]

Hey, thank you for the quick response. You are correct, the pipeline does not specify the --min-identity option and therefore relies on the default of 95%. However, if this is how QUAST is supposed to function then I think the definitions are maybe a bit misleading... Surely, even if the alignment identity is below the minimum value set, it is confusing to be classifying sequences as "fully unaligned" meaning they "have no alignment to the reference sequence" when that is not the case (for example, when they have an alignment with 94% identity and 100% coverage)? I suppose I'm just sad that the classification/labelling of the contigs is not as straightforward as it seems.

I also just went and looked at some differences in the Minimap2 and NUCmer alignment results, and it seems a majority of the differences come from the fact that Minimap2 is classifying a lot of the sequences identified as fully unaligned by NUCmer as "misassembled" and therefore excluding them from the "fully unaligned" dataset. Not sure what to think about that at the moment, but hopefully with some more analysis I will wrap my head around it. Ultimately, the change in QUAST has made an extremely massive difference to my results, and I guess I just need to decide whether I can justify using the older version when the semantics might be misleading to other people interpreting my results (as in, I'm looking for fully unaligned sequences, and I'm not sure I can classify a sequence with 94% identity alignment as fully unaligned).

Thank you anyway for the help!

[alexeigurevich]

Well, I agree, probably we should rename "unlaigned" to something like "not having alignments above min IDY threshold" but it seems too wordy for a .txt report..

[alexeigurevich]

You can also play around with various --min-identity values to find "real" (almost) fully unaligned sequences, e.g., having no matches or matches with identity below 80%.
Note that we use slightly different minimap2 presets for different --min-identity thresholds, since the aligner has some optimisations for closely related and distant genomes.