ad8e7d43125d15461861a22ccb4d99e2b46c88c0
lrnassar
  Fri Aug 21 10:46:04 2026 -0700
Fix TP53 gnomAD PASS filtering and FLOSSIES carrier count, per code review. refs #37399 refs #38139

Two independent correctness fixes from CR #38139:

- tp53AFfrequencies.py: skip non-PASS gnomAD records (mostly AC0, i.e. observed
in nobody after QC). They were entering the present-set and wrongly counting as
present in gnomAD, blocking the absent -> PM2_Supporting rule. 68 missense
variants now correctly get PM2, changing 6 provisional classes at the 5/6 and
-2/-1 boundaries. gnomAD frequency use is PASS-only regardless.

- tp53Flossies.py: count carriers as het + hom by summing the per-population het
and hom counts, instead of allele_count - hom_count. In this FLOSSIES export
allele_count is already het + hom (not the usual het + 2*hom), so the old
formula dropped homozygous carriers (e.g. P72R showed 4031 carriers, not 8372).
No BS2 tier changes in the current export; fixes the mouseover count and a
latent tier bug for low-count variants with homozygotes.

diff --git src/hg/makeDb/scripts/tp53/tp53AFfrequencies.py src/hg/makeDb/scripts/tp53/tp53AFfrequencies.py
index 96b8d49dbc3..b67572ee7ec 100644
--- src/hg/makeDb/scripts/tp53/tp53AFfrequencies.py
+++ src/hg/makeDb/scripts/tp53/tp53AFfrequencies.py
@@ -172,67 +172,76 @@
 
 def classify_and_build_rows(tx, chrom):
     """Read the source gnomAD v4.1 exomes VCF on hg38 via tabix and emit a list
     of classified rows keyed by an immutable hg38 identifier. The hg38 id is used
     as the 'name' field so the hg19 build can look up the same row after liftOver
     and rewrite the display text to reflect hg19 coords. Reading the source VCF
     (not the /gbdb bigBed) gives the per-ancestry faf95 that BA1/BS1 require."""
     region = "{}:{}-{}".format(chrom, tx['txStart'] + 1, tx['txEnd'])
     out = subprocess.run([TABIX, GNOMAD_VCF, region],
                          capture_output=True, text=True, check=True).stdout
     vcf_lines = [ln for ln in out.splitlines() if ln and not ln.startswith('#')]
     print("  {} variants in TP53 region (hg38)".format(len(vcf_lines)))
 
     classified = []   # list of dicts with all fields; hg38 coords fixed
     all_records = []  # every gnomAD variant, for the PM2 "present in gnomAD" set
-    stats = dict(total=len(vcf_lines), BA1=0, BS1=0, PM2=0, skipped=0, multi=0)
+    stats = dict(total=len(vcf_lines), BA1=0, BS1=0, PM2=0, skipped=0, multi=0,
+                 nonpass=0)
     for ln in vcf_lines:
         f = ln.split('\t')
         pos = int(f[1])          # 1-based VCF POS
         ref = f[3]
         alt = f[4]
         if ',' in alt:           # multi-allelic (none expected in v4.1 sites VCF)
             stats['multi'] += 1
             continue
+        # Skip non-PASS records (mostly AC0 = observed in nobody after QC, plus
+        # AS_VQSR / InbreedingCoeff filter failures). gnomAD frequency use is
+        # PASS-only; more importantly a non-PASS variant must not enter the
+        # present-set, or it would wrongly count as "present in gnomAD" and block
+        # the absent -> PM2_Supporting rule in the Provisional track.
+        if f[6] not in ('PASS', '.'):
+            stats['nonpass'] += 1
+            continue
         info = parse_info(f[7])
         c_start = pos - 1
         c_end = c_start + len(ref)
         all_records.append({'chrom': chrom, 'hg38_start': c_start,
                             'hg38_end': c_end, 'ref': ref, 'alt': alt})
         af_global = safe_float(info.get('AF'))
         af_grpmax = safe_float(info.get('AF_grpmax'))
         grpmax_pop = GRPMAX_POP_NAMES.get(info.get('grpmax'), info.get('grpmax'))
         faf_ba1bs1, faf_group = nonfounder_max_faf(info)
         hg38_name = "{}-{}-{}-{}".format(chrom, pos, ref, alt)
 
         code = classify(af_global, af_grpmax, faf_ba1bs1)
         if code is None:
             stats['skipped'] += 1
             continue
         stats[code if code in ('BA1', 'BS1') else 'PM2'] += 1
         classified.append({
             'hg38_name': hg38_name,
             'hg38_start': c_start,
             'hg38_end': c_end,
             'chrom': chrom,
             'ref': ref, 'alt': alt,
             'af_global': af_global, 'faf': faf_ba1bs1, 'faf_group': faf_group,
             'af_grpmax': af_grpmax, 'grpmax_pop': grpmax_pop,
             'code': code,
         })
     print("  classified: BA1={BA1} BS1={BS1} PM2={PM2} skipped={skipped} "
-          "multiallelic_skipped={multi}".format(**stats))
+          "multiallelic_skipped={multi} nonpass_skipped={nonpass}".format(**stats))
     return classified, all_records
 
 
 def write_present_set(all_records, db, outdir):
     """Write the set of every gnomAD variant key ('chrom-pos1-ref-alt') present
     at the TP53 locus, in the coordinates of the requested assembly. The
     Provisional track uses this to tell 'absent from gnomAD' (PM2_Supporting
     applies) apart from 'present but not rare enough to be coded' (no PM2). For
     hg19 the hg38 coords are lifted so the keys match the hg19 Provisional build."""
     path = os.path.join(outdir, "TP53AF_present_{}.txt".format(db))
     keys = []
     if db == 'hg38':
         for r in all_records:
             keys.append("{}-{}-{}-{}".format(
                 r['chrom'], r['hg38_start'] + 1, r['ref'], r['alt']))