Favor Pileup over Full_alignment when SNP has high QUAL and Full alignment incorrectly assigns as Refcall

[wshropshire]

Hello,

I have a similar issue with #318 (although with bacteria/haploid data) where I can see in the VCFs as well as the pileup that under peculiar circumstances, I am seeing near high snp densities where there are pileup variant calls that are correctly assigning an alternate allele but in the full alignment at the same position it's incorrectly typed as as a refcall. Using the same run_clair.sh (v1.0.10) parameters shown below, I get for a serial isolate with the correct call. This is with bacteria data btw, and I found that in some instances when running with --haploid_precise, the GT was being called 0/1 and removed when the evidence for ALT had >.9 AF and high coverage (100X), thus I've dropped it for now. Here's a picture of the region in particular with the two serial isolates aligned to a reference from a different sequence type and the VCF showing the discordant pileup v. full alignment VCs:

CR-0005 Full Alignment VCF:

CR-0005 Pileup VCF:

CR-0063 Full Alignment VCF:

CR-0063 Pileup VCF:

Here's an IGV snapshot:

And command here:
clair3_cmd = [ "run_clair3.sh", f"--bam_fn={os.path.abspath(bam_output)}", f"--ref_fn={os.path.abspath(ref)}", f"--threads={threads}", "--platform=ont", f"--parallel={parallel_path}", f"--model_path={os.path.abspath(model_path)}", f"--output={clair3_output_dir}", f"--sample_name={sample}", "--include_all_ctgs", "--no_phasing_for_fa", "--enable_long_indel" ]

I'm using the r1041_e82_400bps_sup_v500 model. I can provide data if need be. Love your tool though! This and more consistent masking of recombinant regions and should have a decent core genome alignment pipeline.

Best,

Will

[zhengzhenxian]

Hi Will,

Thank you for reporting this.

It appears that the pileup output provided more reliable results for your bacterial data. The full-alignment model may have failed for some cases due to the high density of SNP variants and high-coverage sample. Given the relatively high-quality scores for the failed cases in pileup, you might consider feeding only a lower proportion of pileup calls into the full-alignment model by using the option --var_pct_full=0.3 (the default is 0.7).

[wshropshire]

Okay, just FYI in your argparse help, you mention the default is 0.3 for --var_pct_full. Also, if I increase --snp_min_af from 0.08 to 0.9, would that reduce the likelihood of 0/1 genotypes being called and thus allow me to use the --haploid_precise option?

[aquaskyline]

@wshropshire When the coverage gets higher, it is likely that signals that object a candidate as a true variant get more. Full-alignment was trained to be more sensitive to these signals than the pileup. One might argue that full-alignment should also take the additional supporting signals into account to counterbalance the objecting signals, but a tuning that favors an example might lead to a bunch of failures on other examples. While the input to the Clair3 can be of different species and different coverages, in fact there is no one single perfect solution that fits all usage scenarios.

Thus we wrote Clair3 to be as modulized as possible. Would you think of writing a small downstream processing script that takes both pileup and full-alignment output as inputs? In the cases you have shown, it seems that you would consider those variants called in pileup, albeit failed the full-alignment, but with good enough VAF and coverage support, as a true variant. You can also give those variants supported by both pileup and full-alignment a higher variant quality, and those supported by only pileup a lower quality, but still consider them as true. What I mean is, a customized downstream processing script gives you the best flexibility, while you can maintain the best recall performance at Clair3's raw outputs.

[wshropshire]

Hello Ruibang,

I'll take your suggestion and try to create a new merge module that takes in a number of conditionals to determine when to favor one call vs another. I frankly don't know how to create de novo QUAL scores so I'm not even going to touch that :).

How difficult would this be to train the full model with bacterial data? ONT now has bacterial/methylation aware models for dorado basecalling that have significantly improved it's basecalling just in the past year. Could I or clair3 developers think about creating a similar bacteria model?

[aquaskyline]

Ya, we kept an eye on any new data for model training. Technically, training a new or refining a model is not too challenging. But the quality of truths need to be really good to lead to a better performance after refinement. Are there any bacteria genomes that are with very good quality truth variants?

[wshropshire]

Hey, I have submitted quite a few bacterial genomes to NCBI that I have personally curated myself. I think given that it's been only in the past year I have considered ONT-only assemblies, they all should be hybrid ONT + ILMN assemblies. Regardless, I don't submit assemblies I have worked on unless I've done the proper curation and I can guarantee with some confidence that the assembly is 'error free'. Ryan Wick is another person who would be good to contact about potential 'ground truth' assemblies to train models.

As far as merging the pileup/full alignment models, I've created a merge script that parses through both respective VCFs using pysam and then uses a series of 'ifelse' statements to determine which variant to include in the merged file. It seems to be working fairly well, however, I've kind of hit a dead-end as far as improving results goes. The obstacle now I have narrowed it down to two predominant errors:

(1) false negative VCs with RefCalls occurring in both models at positions where the alignment indicates good support of the ALT allele (AF > .9) and the AF for the REF Allele is low. Here is an example:
Incorrect FN in full alignment:

Incorrect FN in pileup:

Screenshot of IGV pileup:

Note this occurs near a large deletion

(2) Incorrect false positive VC "PASS" calls in one model, but not the other:
Incorrect call in full alignment:

Correct call in pileup (Despite incorrect AF)

Here's another screenshot of an example near a MNP:

And here's a screenshot where both types of errors are occuring across two different samples:

These errors are occurring still predominantly in high SNP density regions and/or INDELs. Since the errors ARE in high density regions, I have been able to mask these regions with Gubbins with partial success, however, that success is still fairly limited.

So I think I've come to the conclusion that the precision/recall of this VF for bacteria really can only improve with better training data specific for bacteria and having a clair3-based model that works with bacteria. I predominantly work with Enterobacterales species, but I also have some other Gram Positive bacteria that I have a graduate student that did both ONT+ILMN assemblies with that you could potentially use. I didn't curate myself, but I can go back and check those as well :)

Let me know if this would be any help. Clair3 has a lot of promise for bacteria and it works well with populations of closely related sequence types, but it still seems to have mostly a specificity issue with more divergent genomes.

[aquaskyline]

Agreed with you that having bacteria data included during model training might lead to a better performance. Dropped you an email. I will need a few promising bacteria genomes that have good quality truths and might also contain the high small variant density regions like the cases you've shown.

[wshropshire]

Hey Ruibang, I didn't get your email yet, but you can reach me at wcshropshire@mdanderson.org.

I pinged Ryan Wick on Slack to see if he could help with providing a diverse collection of bacteria from different taxa since mine are almost exclusively belonging to one order of bacteria. Would be still happy to contribute as well :)

[aquaskyline]

Just in case, I sent via both my work email '@cs.hku.hk' and gmail.