jkbonfield
commented
6 months ago I'll also add that after several days of analysing this, I still believe fundamentally this isn't a useful strategy to take. Hence I'm happy to report my findings to you so you can act on them as you wish, but I am not willing to put time and effort into making PRs because it doesn't fit my own personal testing methodology.
In attempting to fix specific issues found here, I nearly always made the overall performance of my branch get worse. Focusing on a few specific examples is therefore detrimental.
It's useful to know where the tool is working and where it isn't, but what works and what doesn't changes every time we modify the software. I've found it much easier to take the following strategy.
-
Call a known truth set with current bcftools (eg HG002). Compare using bcftools isec to get correct, false positives and false negative sets. You can also split into alleles to handle the GT cases or compare GT separately to ensure the alleles are correctly recorded.
-
Do the same with a modified bcftools.
-
Then compare specific types of failure from 1 and 2 above. Eg do a bcftools isec on old-FP.vcf and new-FP.vcf to see which false positives we've fixed and which new ones we've now made. Similarly for false negatives.
-
For any comparison in 3 above, subdivide into type of variant. Eg are insertions more error prone than deletions? Are we improving overall but at the detrimental cost of GT 1/2 calls? Etc. Basically we can profile the changes and get an overall view of what's going on, with statistically significant numbers rather than 1-2 specific tests designed to demonstrate a specific issue. This avoids over-correcting to a small hand-crafted training set.
-
From any stage above, we can plot false positive vs false negative against a specific
QUAL>=$xfilter. This permits us to see how the trade-off works, as well as where the baseline starts (QUAL>=0) for maximum sensitivity. -
Finally given most of the published HG002 data is quite deep, use
samtools -sto produce shallower copies of the same data. This allows us to then verify the performance on both deep and shallow data, to see if we've introduced more biases there.
I've spent weeks, probably months overall, tinkering with bcftools and I've found the above to be vastly more useful that attempting to maintain a curated set of interesting cases. Indeed I can recreate any number of interesting cases simply using the above procedure, with the added benefit of (mostly) knowing what the right answer is.
(What I haven't done is analyse all the HG00* samples to produce a small population for joint analysis.)
pd3
Hello,
In going through all the tests here I've found numerous mistakes in this test corpus. Many "expected" values show one allele only of a two allele variant. This is mostly true for the CHM1/CHM13 calls, where we already have the SynDip truth set which is (usually) correct. It feels like the "truth" here was simply copied from whatever one specific version of bcftools once produced, rather than being curated. As such I think it's not particularly useful as any improvement in accuracy is reflected as a higher failure rate in this corpus. I did find one SynDip case where I disagreed with the published truth set however.
I spent several days manually going through this lot. Incase it's useful to you I may as well share the findings. Sorry it's quite long, but I've tried to be as detailed as possible in analysing each failed test. The ones that are absent from this list are ones where I agree mostly (bar differing AD numbers, which is a different problem with the tests).