fix: multi-allelic statistics and capture filtering in variant-info - #40
Merged
Merged
Conversation
variant-info intersected the per-variant bitmaps with the sample filter only, so a sample with a call outside its own capture BED was listed as a carrier while query() excluded it from AC, AN and the genotype tallies. Use the same capture- and ploidy-aware eligible set as query(), so both commands agree. On the CSVS hs37d5 database this accounts for every reported mismatch, e.g. at 1:6508625 variant-info listed 120 het / 8 hom / 7 alt against N_HET=95 / N_HOM_ALT=7 / N_FAIL=0 from query.
…elic sites N_HOM_REF was the residual n_eligible - N_HET - N_HOM_ALT - N_FAIL - N_NO_COVERAGE. The eligible set depends only on the position, so at a site with several ALT alleles every sample carrying another allele was counted as hom-ref for this one. At chr13:44995309 G>A,T this reported N_HOM_REF=165 for G>A while 36 of those samples carry G>T. The tallies for one allele now leave out eligible samples that carry only another allele at the position. At a multi-allelic site the five categories therefore add up to less than n_eligible, by exactly that number of samples; biallelic sites are unchanged. AC, AN and AF are not affected. The documented invariant is updated accordingly. Coverage evidence (--min-pass, --min-observed, --min-quality-evidence) is now judged per position: a call for any allele shows the position was sequenced, and such a sample is never reported as N_NO_COVERAGE. The per-allele computation, previously repeated in query, region, batch, dump and annotate, now lives in QueryEngine._variant_stats. annotate also excludes carriers of stored alleles when the requested allele is absent from the database. The test oracle had the same residual and is corrected; a multi-allelic cohort is added to the oracle tests.
Most positions hold one allele, where the pooled evidence equals the row's own bitmaps. Pass None instead of building it, and index stored alleles by position in annotate rather than scanning them per record.
… duplicates normalize_vcf.sh ran bcftools norm without -m, so multi-allelic records were kept whole and their indels were not trimmed per allele: the 2:136592357 CAAAAAAA>C,CAAAAAAAAAAAA record stored its insertion as CAAAAAAA>CAAAAAAAAAAAA instead of C>CAAAAA, splitting one variant into two representations across samples. It also passed -d twice; the second value, "both", treats any two indels at one position as duplicates, so once the record is split one of the two alleles was discarded. Use -m -both -d exact. bcftools refuses to combine -m and -d before 1.20, so the script now checks the version and exits with a clear message on older releases. Closes #39
Codecov Report❌ Patch coverage is
Additional details and impacted files@@ Coverage Diff @@
## master #40 +/- ##
==========================================
+ Coverage 88.69% 89.40% +0.71%
==========================================
Files 22 22
Lines 3077 3077
Branches 482 488 +6
==========================================
+ Hits 2729 2751 +22
+ Misses 229 207 -22
Partials 119 119 ☔ View full report in Codecov by Harness. 🚀 New features to boost your workflow:
|
query prints nothing for a technology that does not cover the position, so the capture-index check now looks for the 'No variants found' message. State that --min-covered is still evaluated per allele at build time, while --min-quality-evidence counts carriers of any allele.
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Summary
Three independent problems around multi-allelic sites and capture regions, reported on the CSVS hs37d5 database.
1.
variant-infolisted carriers outside their capture BEDvariant_info()intersected the per-variant bitmaps with the sample filter only, so a sample with a call outside its own BED was listed as a carrier whilequery()excluded it from AC, AN and the genotype tallies. It now uses the same capture- and ploidy-aware eligible set asquery().On
afquery_db_hs37d5this accounts for every mismatch reported:variant-info)queryBehaviour change: the guide previously described the extra off-target carriers as expected; it is rewritten, and the debugging table no longer relies on that difference.
2.
N_HOM_REFcounted carriers of other alleles at multi-allelic sitesN_HOM_REFwas the residualn_eligible - N_HET - N_HOM_ALT - N_FAIL - N_NO_COVERAGE. The eligible set depends only on the position, so at a site with several ALT alleles every carrier of another allele was counted as hom-ref for this one. The tallies for an allele now leave those samples out.n_eligible. At a multi-allelic site they fall short by exactly the number of eligible samples carrying only another allele. AC, AN and AF are not affected. Documentation updated.--min-pass,--min-observedand--min-quality-evidencenow count calls for any allele at the position, and a sample with such a call is never reported asN_NO_COVERAGE.query, region, batch,dumpandannotate, now lives inQueryEngine._variant_stats.annotatealso excludes carriers of stored alleles when the requested allele is absent from the database.--min-coveredgate (filtered_bitmap) is still evaluated per allele. Carriers of other alleles are removed from it at query time, but at a multi-allelic site a rare allele's row can still report true hom-ref samples asN_NO_COVERAGEwhile the common allele's row counts them as hom-ref. Making it per position requires rebuilding the database.At 13:44995309
N_HOM_REFgoes from 165 (G>A) / 421 (G>T) to 141 for both, the number of G/G samples: 457 eligible minus the 316 carrying A, T or both (12 samples carry both). The commit message for this change says 36 samples were misclassified for G>A; the correct figure is 24, because the 12 samples carrying both alleles were already counted as carriers of A.3.
normalize_vcf.shdid not split multi-allelic recordsbcftools normran without-m, and with-dpassed twice; the effective-d bothtreats two different indels at one position as duplicates. Now-m -both -d exact. bcftools refuses that combination before 1.20 (checked with 1.18 and 1.19), so the script checks the version. VCFs normalized with the old script need re-normalizing and the database rebuilding; documented in the preprocessing guide.Closes #39
Validation
master.tests/test_multiallelic.py: point, region, batch (one allele requested), batch-multi, dump, annotate, variant-info and--min-passat a multi-allelic site. New oracle cohort with a1/2carrier, partial capture and a FILTER failure.masteragainst this branch onafquery_db_hs37d5, whole chr22 (585,928 rows, 116,802 at multi-allelic positions):N_HOM_REFchanges, on 114,232 multi-allelic rows; the invariant holds on every row;--min-pass 1: no biallelic row changes; AC/AN unchanged;N_NO_COVERAGEdrops on 2,292 multi-allelic rows.master, 18.1 s here (allele pooling at multi-allelic positions).normalize_vcf.shrun with bcftools 1.20 and 1.24 on the 2:136592357 record: two biallelic records,CAAAAAAA>CandC>CAAAAA, none dropped. Clear error with 1.18 and 1.19.mkdocs build --strictpasses.No database rebuild is needed for fixes 1 and 2; they apply at query time.