Repository navigation
fix: consistent allele accounting at multi-allelic sites, failed calls out of AN - #41
Merged
Merged
Conversation
…c sites The build stores filtered_bitmap per ALT allele, counting only that allele's quality-passing carriers. At a multi-allelic site a panel with carriers of one allele but none of another had its non-carriers reported as N_NO_COVERAGE for the second allele and as N_HOM_REF for the first. At positions holding several alleles the --min-covered gate is now recomputed at query time from the pooled quality_pass bitmaps, so every allele of the site gets the same verdict. Single-allele positions keep using the stored bitmap, which is equivalent there.
…abase An input ALT allele with no stored row got N_HOM_REF = eligible minus carriers of the stored alleles, with N_NO_COVERAGE always 0. Under --min-pass, --min-observed or --min-quality-evidence the same samples were uncovered for the stored alleles of the record and hom-ref for the absent one. When other alleles are stored at the position, the absent allele is now computed like any other, with no carriers and the pooled site evidence. Positions with nothing stored keep every eligible sample as hom-ref.
A sample with a failed call (FILTER!=PASS) at a position was left out of AC but counted as reference in AN, so AF was computed as if every failed sample were hom-ref. At multi-allelic sites the alleles did not add up either: AN exceeded the reference alleles plus the AC of all ALT alleles by the alleles of the failed samples. A failed call for any allele at a position now makes the sample's genotype there unknown, for every allele of the site: it is counted in N_FAIL, which is the same on every line of the position, and left out of AC and AN. AN therefore equals the reference alleles plus the AC of all ALT alleles. When every eligible sample failed, AN is 0 and AF is undefined (NA in query output, an empty cell in dump, "." in annotate). annotate wrote AFQUERY_N_FAIL as the sum of the per-allele counts, so a sample masked as ./. at a failed multi-allelic site was counted once per allele. It now reports the per-site N_FAIL. variant-info lists the failed samples for every allele of the position, matching N_FAIL. The test fixture stored some samples in both het and fail bitmaps, which a real build never does; it and its source VCFs now agree. BREAKING CHANGE: AN no longer includes samples with a failed call, so AN and AF change at every position where N_FAIL > 0. N_FAIL counts failed samples per position rather than per allele.
Failed calls are now left out of AN. Some analyses still want the previous estimate, in which every failed sample is a hom-ref non-carrier and AF is a conservative lower bound. --fail-as-ref on query, annotate and dump (fail_as_ref in the Python API and SampleFilter) counts those samples as hom-ref in AN. AC, N_FAIL and the genotype tallies are unchanged. The AFQUERY_AN header of annotated VCFs states which convention was used.
Ploidy comes from the sample's sex, and each allele's AC is counted on its own. A call carrying more alleles than the sample has, such as 1/2 on chrX non-PAR in a male without fixploidy, then added one allele per ALT, so the AC of the site exceeded the alleles that sample contributes to AN. Ingestion now checks each sample's calls per position against its ploidy: two different ALT alleles on a haploid position, or more than two copies on a diploid one, make every call of the sample there a failed call, counted in N_FAIL and left out of AC and AN. create-db and update-db --add-samples log a warning per sample listing the positions. ingest_all and ingest_sample now take the genome build (and the sample's sex) to decide ploidy. Existing databases need a rebuild to apply it. The preprocessing guide claimed that diploid male chrX calls inflate AN; AN always used the sample's sex, and the guide now says so.
Text and TSV query output always include N_NO_COVERAGE, but most examples stopped at N_FAIL, and two TSV examples printed AF with five decimals instead of six. The tutorial examples were re-run against the demo dataset. The concepts page also said ingestion is PASS-only, while failed calls are stored and reported as N_FAIL.
An input allele at a position where no sample has a call still reported every eligible sample as hom-ref, even under --min-pass, --min-observed, --min-quality-evidence or a build-time --min-covered. A position with no evidence at all was thus trusted more than one with a single carrier. Such alleles are now computed like any other absent allele, from an empty site: under an active gate, the non-carriers of partially-covered techs are N_NO_COVERAGE. The absent-allele branch also reuses the site evidence already pooled for multi-allelic positions.
The check keyed alleles by their ALT string alone, so a haploid sample carrying a SNV G>A and a deletion GA>A at one position was not flagged, and summed copies across records, so a repeated identical record made a diploid 1/1 look like four copies. Alleles are now keyed by (REF, ALT), and a repeated record counts once.
AN is computed per variant in _variant_stats, so the AN that _compute_eligible still returned was only used to test for an empty eligible set, at the cost of extra bitmap work per position.
A failed call is not only FILTER!=PASS: genotypes inconsistent with the sample's ploidy are recorded as failed too. Refer to N_FAIL instead.
Codecov Report❌ Patch coverage is Additional details and impacted files@@ Coverage Diff @@
## master #41 +/- ##
==========================================
+ Coverage 89.40% 89.58% +0.18%
==========================================
Files 22 22
Lines 3077 3131 +54
Branches 488 501 +13
==========================================
+ Hits 2751 2805 +54
+ Misses 207 206 -1
- Partials 119 120 +1 ☔ View full report in Codecov by Harness. 🚀 New features to boost your workflow:
|
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
2 participants
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
At every position, AN now equals the reference alleles plus the AC of all its ALT alleles, and every allele of a multi-allelic site gets the same verdict.
N_FAILand left out of AC and AN for every allele there.N_FAILis now a per-site value. When every eligible sample failed, AN is 0 and AF is undefined:NAin query output, an empty cell indump,.inannotate.--fail-as-refonquery,annotateanddump(fail_as_refin the Python API) counts failed calls as hom-ref in AN, as before.AFQUERY_N_FAILcounted a sample once per ALT allele; it now counts each failed sample once per site.--min-coveredis judged per position at multi-allelic sites. Before, a panel could be "uncovered" for one allele and covered for another at the same position.annotateapplies coverage evidence to input alleles absent from the database: they share the evidence of the stored alleles at their position, and at a position with nothing stored, every active gate finds no evidence.1/2on chrX in a male without fixploidy, or on chrM) are recorded as failed calls at ingestion, with a warning listing the positions. Alleles are told apart by REF and ALT.Breaking changes
N_FAIL > 0. Use--fail-as-reffor the previous AN.N_FAILcounts failed samples per position rather than per allele.dump,AFandAF_<group>are empty when the matching AN is 0 (AF_<group>used to be0.0).annotatenow reports novel positions inside partially-covered regions asN_NO_COVERAGEfor those techs.Test plan
pytest: 640 passed. Includes new multi-allelic scenarios, an extended independent oracle (site-level FAIL, ploidy check,fail_as_ref) and ingest unit tests.mkdocs build --strictpasses.