6b285a53b036b309e3c7a9b61d3741731088a172 lrnassar Fri Jun 12 02:35:01 2026 -0700 varFreqs: switch affectedAF/backgroundAF from max-across-cohorts to pooled sum(AC)/sum(AN) so the rate matches the carrier count scale. Per-arm AN is derived as round(AC/AF) when both are reported. An optional "default_an" column was added to databases.tsv so AF-only cohorts (ABraOM, ALFA) can synthesize a denominator from their cohort size; without it those cohorts had been silently dropped from the pooled rate. New affectedAN and backgroundAN columns expose the pool denominator. The mouseOver now reads "Affected AC/AN: 33238 / 213153" so the ratio is visible. Per-arm cohorts that ship only AC and no default_an (MGRB, GREGoR AC_AFFECTED/UNAFFECTED/UNKNOWN, AllOfUs per-population) are still listed in affectedCohorts/backgroundSources but contribute 0 to the pool, preserving the invariant pool_AF <= 1. The build pipeline is unchanged: re-run vcfToBigBed.py --split-affected against the existing merged.annotated.vcf.gz. refs #36642 diff --git src/hg/makeDb/scripts/varFreqs/vcfToBigBed.py src/hg/makeDb/scripts/varFreqs/vcfToBigBed.py index 1cf9dbe9c73..35c96b421bb 100755 --- src/hg/makeDb/scripts/varFreqs/vcfToBigBed.py +++ src/hg/makeDb/scripts/varFreqs/vcfToBigBed.py @@ -139,35 +139,43 @@ else os.path.join(scripts_dir, db_file) with open(db_path) as f: for line in f: line = line.strip() if not line or line.startswith('#'): continue parts = line.split('\t') if len(parts) < 5: print(f"WARNING: skipping malformed line: {line}", file=sys.stderr) continue key, name, vcf, ac_field, af_field = ( parts[0], parts[1], parts[2], parts[3], parts[4]) is_disease = int(parts[5]) if len(parts) > 5 else 0 disease_role = parts[6].strip() if len(parts) > 6 else "" + default_an = 0 + if len(parts) > 7 and parts[7].strip(): + try: + default_an = int(parts[7].strip()) + except ValueError: + print(f"WARNING: bad default_an for {key}: {parts[7]}", + file=sys.stderr) databases[key] = { "name": name, "vcf": vcf, "ac_field": ac_field, "af_field": af_field, "is_disease": is_disease, "disease_role": disease_role, + "default_an": default_an, "pops": [], } pop_path = pop_file if os.path.isabs(pop_file) \ else os.path.join(scripts_dir, pop_file) with open(pop_path) as f: for line in f: line = line.strip() if not line or line.startswith('#'): continue parts = line.split('\t') if len(parts) < 5: continue db_key = parts[0] phenotype = parts[5].strip() if len(parts) > 5 else "" @@ -276,30 +284,55 @@ """Load pre-extracted TSV into dict keyed by pos:ref:alt.""" path = os.path.join(extract_dir, db_key, f"{chrom}.tsv") data = {} if not os.path.exists(path): return data with open(path) as f: for line in f: parts = line.rstrip('\n').split('\t') if len(parts) < 3: continue key = f"{parts[0]}:{parts[1]}:{parts[2]}" data[key] = parts[3:] return data +def _pool_arm(ac_val, af_val, default_an): + """Compute pooled (AC, AN) contribution for one cohort arm. + + Used by the affected and background pooled-AF calculations. Returns + (0, 0) when we can't determine AN, so the pool denominator never + includes a cohort's carriers without also including its allele number + -- the resulting pooled AF stays well-defined and bounded. + + Strategies, in order: + 1. Both AC and AF present with AF > 0: AN = round(AC / AF) (typical case). + 2. AF present but AC empty: synthesize AC = round(AF * default_an) + and use default_an as AN (e.g. ALFA, ABraOM, which ship only AF). + 3. AC present but AF empty/0: use default_an as AN (e.g. MGRB if it + had a configured default_an). + 4. None of the above: return (0, 0), arm does not contribute. + """ + if ac_val is not None and af_val is not None and af_val > 0: + return ac_val, max(1, round(ac_val / af_val)) + if af_val is not None and default_an > 0: + return max(0, round(af_val * default_an)), default_an + if ac_val is not None and default_an > 0: + return ac_val, default_an + return 0, 0 + + def process_chromosome(args): """Phase 2: Build the affected and background BEDs for one chromosome from the annotated VCF + pre-extracted per-cohort data. Two output rows are possible per variant, sharing one schema: - affected BED: variant seen in any affected/case arm of a disease cohort - background BED: variant seen in a population cohort or an unaffected/ control/unknown arm ("all other variants") A variant present in both groups is written to both (overlap is intended, so the case-vs-background comparison works). With split=False a single BED is written instead (one row per variant, score = max of the two summaries); used for tracks with no disease cohorts such as the genotyping-array combined track.""" chrom, annotated_vcf, databases, extract_dir, output_dir, split = args @@ -351,179 +384,190 @@ key = f"{pos}:{ref}:{alt}" consequence, gene, transcript, aa_change, dna_change = \ parse_bcsq(bcsq) r, g, b = get_color(bcsq) pos_int = int(pos) start = pos_int - 1 end = start + len(ref) ref_d = ref[:17] + "..." if len(ref) > 20 else ref alt_d = alt[:17] + "..." if len(alt) > 20 else alt name = f"{ref_d}>{alt_d}" var_type = get_vartype(ref, alt) - # Affected/case summary (affected arms + role=affected whole cohorts) - affected_af = 0.0 + # Pooled affected/case summary: sum AC, sum AN, AF = AC/AN. + # Switched from max-across-cohorts (which was dominated by tiny + # cohorts like GA4K when they reported high local AF) to a + # population-weighted ratio so the AF matches the AC scale. affected_ac = 0 + affected_an = 0 affected_cohorts = [] # Background summary = population cohorts + unaffected/control/unknown - # arms of disease cohorts ("all other variants"). - background_af = 0.0 + # arms of disease cohorts ("all other variants"), same pooling. background_ac = 0 + background_an = 0 background_sources = [] - db_ac_af = [] # per-database AC, AF - pop_ac_af = [] # per-population AC, AF (written AFTER all db fields) + db_ac_af = [] # per-database AC, AF (raw, for output columns) + pop_ac_af = [] # per-population AC, AF (raw, for output columns) for db_key, db_info in databases.items(): values = freq_data.get(db_key, {}).get(key, []) ac = values[0] if len(values) > 0 and values[0] not in \ (".", "") else "" af = values[1] if len(values) > 1 and values[1] not in \ (".", "") else "" is_disease_db = db_info.get("is_disease", 0) disease_role = db_info.get("disease_role", "") + default_an = db_info.get("default_an", 0) af_val = None if af: try: af_val = float(af) except ValueError: af_val = None ac_val = None if ac: try: ac_val = int(ac) except ValueError: ac_val = None - # Track this cohort's contribution to each group so the - # affectedCohorts / backgroundSources lists are accurate. + ac_add, an_add = _pool_arm(ac_val, af_val, default_an) + cohort_observes = (ac_val is not None) or (af_val is not None) + + # Track this cohort's appearance in each group's source list. + # A cohort that observes the variant but lacks a usable AN + # (e.g. MGRB ships AC only, GREGoR per-arm ships AC only) is + # still listed but contributes 0 to the pool. Future work: + # add default_an entries for these cohorts/arms. hits_affected = False hits_background = False if is_disease_db: - # Unified AC/AF slot: only meaningful when the whole cohort - # has a known role (e.g. GA4K = affected). Otherwise the - # per-arm populations below carry the phenotype signal. if disease_role == "affected": - if af_val is not None: - affected_af = max(affected_af, af_val) - hits_affected = True - if ac_val is not None: - affected_ac += ac_val + affected_ac += ac_add + affected_an += an_add + if cohort_observes: hits_affected = True elif disease_role == "unaffected": - if af_val is not None: - background_af = max(background_af, af_val) - hits_background = True - if ac_val is not None: - background_ac += ac_val + background_ac += ac_add + background_an += an_add + if cohort_observes: hits_background = True else: - # Population cohort feeds the background summary. - if af_val is not None: - background_af = max(background_af, af_val) - hits_background = True - if ac_val is not None: - background_ac += ac_val + background_ac += ac_add + background_an += an_add + if cohort_observes: hits_background = True db_ac_af.extend([ac, af]) for i, pop in enumerate(db_info["pops"]): idx = 2 + i * 2 pop_ac = values[idx] if len(values) > idx and \ values[idx] not in (".", "") else "" pop_af = values[idx + 1] if len(values) > idx + 1 and \ values[idx + 1] not in (".", "") else "" pop_ac_af.extend([pop_ac, pop_af]) pop_af_val = None if pop_af: try: pop_af_val = float(pop_af) except ValueError: pop_af_val = None pop_ac_val = None if pop_ac: try: pop_ac_val = int(pop_ac) except ValueError: pop_ac_val = None + # Per-arm default_an would let GREGoR per-arm rows pool + # cleanly; for now they fall through with default 0. + pop_default_an = pop.get("default_an", 0) + pop_ac_add, pop_an_add = _pool_arm( + pop_ac_val, pop_af_val, pop_default_an) + pop_observes = (pop_ac_val is not None) or \ + (pop_af_val is not None) pheno = pop.get("phenotype", "") if is_disease_db and pheno == "affected": - if pop_af_val is not None: - affected_af = max(affected_af, pop_af_val) - hits_affected = True - if pop_ac_val is not None: - affected_ac += pop_ac_val + affected_ac += pop_ac_add + affected_an += pop_an_add + if pop_observes: hits_affected = True elif is_disease_db and pheno in ("unaffected", "unknown"): # Unaffected relatives, controls, and unknown-phenotype - # individuals are all "not clearly affected". - if pop_af_val is not None: - background_af = max(background_af, pop_af_val) - hits_background = True - if pop_ac_val is not None: - background_ac += pop_ac_val + # individuals all feed the background. + background_ac += pop_ac_add + background_an += pop_an_add + if pop_observes: hits_background = True elif not is_disease_db: - # Ancestry population of a population cohort. - if pop_af_val is not None: - background_af = max(background_af, pop_af_val) - hits_background = True - if pop_ac_val is not None: + # Ancestry breakdown of a population cohort. The + # unified row above already pooled the cohort if it + # had AC+AF, so we deliberately don't double-count + # the per-pop AC/AN here. Per-pop AC and AF still + # write to their own bigBed columns above. + if pop_observes: hits_background = True if hits_affected: affected_cohorts.append(db_key) if hits_background: background_sources.append(db_key) + # Compute pooled allele frequencies. + affected_af = (affected_ac / affected_an) if affected_an > 0 else 0.0 + background_af = (background_ac / background_an) \ + if background_an > 0 else 0.0 + in_affected = 1 if (affected_ac > 0 or affected_af > 0) else 0 # Track length extremes for data-driven length filter ranges. ref_len = len(ref) alt_len = len(alt) var_len = alt_len - ref_len if ref_len > stats["max_ref_len"]: stats["max_ref_len"] = ref_len if alt_len > stats["max_alt_len"]: stats["max_alt_len"] = alt_len if var_len < stats["min_var_len"]: stats["min_var_len"] = var_len if var_len > stats["max_var_len"]: stats["max_var_len"] = var_len # Shared columns (score at index 4 is filled in per output below). base = [ chrom_name, str(start), str(end), name, "0", "+", str(start), str(end), f"{r},{g},{b}", ref, alt, str(ref_len), str(alt_len), str(var_len), var_type, normalize_consequence(consequence), gene, transcript, aa_change, dna_change, f"{affected_af:.6f}" if affected_af > 0 else "", str(affected_ac) if affected_ac > 0 else "", + str(affected_an) if affected_an > 0 else "", ",".join(affected_cohorts), f"{background_af:.6f}" if background_af > 0 else "", str(background_ac) if background_ac > 0 else "", + str(background_an) if background_an > 0 else "", ",".join(background_sources), str(in_affected), ] # Database AC/AF first, then population AC/AF — must match autoSql order base.extend(db_ac_af) base.extend(pop_ac_af) has_affected = affected_af > 0 or affected_ac > 0 has_background = background_af > 0 or background_ac > 0 if split: if has_affected: row = list(base) row[4] = str(min(1000, int(affected_af * 1000))) out_aff.write("\t".join(row) + "\n") @@ -574,36 +618,40 @@ f.write(' uint thickEnd; "Thick end"\n') f.write(' uint reserved; "Color by consequence"\n') # Variant info f.write(' lstring ref; "Reference allele"\n') f.write(' lstring alt; "Alternate allele"\n') f.write(' int refLen; "Reference length"\n') f.write(' int altLen; "Alternate length"\n') f.write(' int varLen; "Length change (alt-ref)"\n') f.write(' string varType; "Type (SNV/INS/DEL/MNV)"\n') # Consequence f.write(' string consequence; "Consequence"\n') f.write(' string gene; "Gene"\n') f.write(' string transcript; "Transcript"\n') f.write(' lstring aaChange; "AA change"\n') f.write(' lstring dnaChange; "DNA change"\n') - # Frequency summaries (shared by the affected and background tracks) - f.write(' string affectedAF; "Max allele frequency in affected/case individuals"\n') + # Frequency summaries (shared by the affected and background tracks). + # AF is pooled across contributing arms (sum AC / sum AN), not the + # max across arms, so the AF matches the AC and AN scale. + f.write(' string affectedAF; "Pooled allele frequency in affected/case individuals (sum AC / sum AN)"\n') f.write(' string affectedAC; "Summed allele count in affected/case individuals"\n') + f.write(' string affectedAN; "Summed allele number in affected/case individuals (pool denominator)"\n') f.write(' string affectedCohorts; "Disease cohorts contributing affected/case carriers"\n') - f.write(' string backgroundAF; "Max allele frequency in population cohorts + unaffected/control individuals"\n') + f.write(' string backgroundAF; "Pooled allele frequency in population cohorts + unaffected/control individuals (sum AC / sum AN)"\n') f.write(' string backgroundAC; "Summed allele count in population cohorts + unaffected/control individuals"\n') + f.write(' string backgroundAN; "Summed allele number in population cohorts + unaffected/control individuals (pool denominator)"\n') f.write(' string backgroundSources; "Cohorts contributing to the background (population + unaffected)"\n') f.write(' uint inAffected; "1 if seen in an affected/case arm, else 0"\n') # Per-database AC/AF for db_key, db_info in databases.items(): f.write(f' string {db_key}AC;' f' "{db_info["name"]} AC"\n') f.write(f' string {db_key}AF;' f' "{db_info["name"]} AF"\n') # Per-population AC/AF for db_key, db_info in databases.items(): for pop in db_info["pops"]: f.write(f' string {db_key}AC_{pop["key"]};' f' "{db_info["name"]} {pop["name"]} AC"\n') f.write(f' string {db_key}AF_{pop["key"]};' f' "{db_info["name"]} {pop["name"]} AF"\n') @@ -669,40 +717,44 @@ "altLen": ("Alternate Length", 1, len_stats["max_alt_len"]), "varLen": ("Length Change", len_stats["min_var_len"], len_stats["max_var_len"]), } for fld, (label, lo, hi) in len_ranges.items(): f.write(f" filterByRange.{fld} on\n") f.write(f" filterLabel.{fld} {label}\n") f.write(f" filter.{fld} {lo}:{hi}\n") f.write(f" filterLimits.{fld} {lo}:{hi}\n") # Affected and background frequency summaries (both tracks carry both, # so e.g. the Affected track can be filtered to variants rare in the # background). String range filters mirror the per-database AF/AC fields. f.write(" # Affected/case frequency summary\n") f.write(" filterByRange.affectedAF on\n") - f.write(" filterLabel.affectedAF Affected/case AF\n") + f.write(" filterLabel.affectedAF Affected/case AF (pooled)\n") f.write(" filterLimits.affectedAF 0:1\n") f.write(" filterByRange.affectedAC on\n") f.write(" filterLabel.affectedAC Affected/case AC\n") + f.write(" filterByRange.affectedAN on\n") + f.write(" filterLabel.affectedAN Affected/case AN (pool denominator)\n") f.write(" # Background (population + unaffected) frequency summary\n") f.write(" filterByRange.backgroundAF on\n") - f.write(" filterLabel.backgroundAF Background AF (population + unaffected)\n") + f.write(" filterLabel.backgroundAF Background AF (pooled)\n") f.write(" filterLimits.backgroundAF 0:1\n") f.write(" filterByRange.backgroundAC on\n") f.write(" filterLabel.backgroundAC Background AC (population + unaffected)\n") + f.write(" filterByRange.backgroundAN on\n") + f.write(" filterLabel.backgroundAN Background AN (pool denominator)\n") f.write(" # Affected/case membership flag\n") f.write(" filterByRange.inAffected on\n") f.write(" filterLabel.inAffected Seen in an affected/case arm (1=yes, 0=no)\n") f.write(" filter.inAffected 0:1\n") f.write(" filterLimits.inAffected 0:1\n") # Per-database AF and AC filters f.write(" # Per-database AF filters\n") for db_key, db_info in databases.items(): f.write(f" filterByRange.{db_key}AF on\n") f.write(f" filterLabel.{db_key}AF " f"{db_info['name']} AF\n") f.write(" # Per-database AC filters\n") for db_key, db_info in databases.items():