5ad55adbb6a5cc72a393700130584aa87fef2c89 lrnassar Tue Jun 30 06:15:44 2026 -0700 varFreqs: add Top 3 source AFs to mouseOvers; audit excludes SGDP and SVatalog. refs #36642 Adds a Top 3 source AFs ranking to the varFreqsAffected and varFreqsBackground mouseOvers. Alongside the pooled allele frequency, the mouseOver now lists the three cohorts/arms with the highest per-source AF, formatted as "Source (AF), Source (AF), Source (AF)". Disease cohorts with phenotype splits carry the arm label (SPARK ASD, SCHEMA case, GREGoR unaffected); population cohorts use the bare key. Per-population sub-ancestries are deliberately excluded so a high sub-pop AF cannot crowd out actual project-level signals. vcfToBigBed.py adds a top_n_source_afs helper, collects per-arm AFs into affected_arm_afs / background_arm_afs, and emits two new fields topAffectedSources and topBackgroundSources. AS schema field count 163 -> 165. An AF-distribution sweep across all 28 source cohorts identified SGDP and SVatalog as encoding allele counts per genotyped individual (small N, AF defaults near 0.5), making their per-source AF unreliable for the ranking. Adds a skip_top_ranking column (col 9) to databases.tsv, set to 1 for SGDP and SVatalog, and gates the per-arm AF append in vcfToBigBed.py on this flag. Both cohorts still contribute to pooled backgroundAC/AN/AF and still appear in backgroundSources; they are only suppressed from the Top 3. Description pages varFreqsAffected.html and varFreqsBackground.html document the ranking; the latter also documents the SGDP/SVatalog exclusion. Build documentation in varFreqs.txt is updated. diff --git src/hg/makeDb/scripts/varFreqs/vcfToBigBed.py src/hg/makeDb/scripts/varFreqs/vcfToBigBed.py index 01bd3dd313f..7295749f6a6 100755 --- src/hg/makeDb/scripts/varFreqs/vcfToBigBed.py +++ src/hg/makeDb/scripts/varFreqs/vcfToBigBed.py @@ -146,36 +146,44 @@ 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) + skip_top_ranking = 0 + if len(parts) > 8 and parts[8].strip(): + try: + skip_top_ranking = int(parts[8].strip()) + except ValueError: + print(f"WARNING: bad skip_top_ranking for {key}: {parts[8]}", + 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, + "skip_top_ranking": skip_top_ranking, "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 "" @@ -284,30 +292,73 @@ """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 +# Short display label per disease-cohort arm key, used in the topAffectedSources +# and topBackgroundSources mouseOver fields so a reader of e.g. +# "SPARK ASD (0.144)" knows which arm of SPARK contributed that AF. Bare +# population-cohort entries use just the db_key. +ARM_DISPLAY = { + "AUT": "ASD", + "NON_AUT": "non-ASD", + "CASE": "case", + "CTRL": "ctrl", + "AFF": "affected", + "UNA": "unaffected", + "UNK": "unknown", +} + + +def source_label(db_key, pop_key=None): + """Display label for a (cohort, arm) contribution to top-N source AFs.""" + if pop_key is None: + return db_key + return f"{db_key} {ARM_DISPLAY.get(pop_key, pop_key)}" + + +def top_n_source_afs(arm_afs, n=3): + """Format the top-N (source, AF) entries by AF value. + + arm_afs is a list of (source_label, af_value) tuples. Sources appearing + more than once (shouldn't, but defensive) are deduped to the max AF. + Ties are broken alphabetically by source label. Returns "" when no + entry has AF > 0.""" + if not arm_afs: + return "" + by_source = {} + for src, af in arm_afs: + if af is None or af <= 0: + continue + if src not in by_source or af > by_source[src]: + by_source[src] = af + if not by_source: + return "" + top = sorted(by_source.items(), key=lambda x: (-x[1], x[0]))[:n] + return ", ".join(f"{src} ({af:.6f})" for src, af in top) + + 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. @@ -396,44 +447,58 @@ name = f"{ref_d}>{alt_d}" var_type = get_vartype(ref, alt) # 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"), same pooling. background_ac = 0 background_an = 0 background_sources = [] + # Per-arm (label, AF) contributions used to compute the top-N + # source AFs in the mouseOver. Disease-cohort arms get the arm + # label appended ("SPARK ASD"); population cohorts are bare. + # Per-population sub-ancestries of population cohorts are + # excluded so a high sub-pop AF in a reference panel does not + # crowd out actual project signals. + affected_arm_afs = [] + background_arm_afs = [] 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) + # When set, the cohort's per-source AF is unreliable for the + # Top-3 mouseOver ranking (e.g. SGDP and SVatalog encode AC/AN + # per genotyped individual, so AF defaults near 0.5). The + # cohort still feeds the pooled AC/AN/AF and appears in the + # Sources list, just not in the Top-3 ranking. + skip_top = db_info.get("skip_top_ranking", 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 ac_add, an_add = _pool_arm(ac_val, af_val, default_an) @@ -441,40 +506,49 @@ # 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: if disease_role == "affected": affected_ac += ac_add affected_an += an_add if cohort_observes: hits_affected = True + if af_val is not None and af_val > 0 and not skip_top: + affected_arm_afs.append( + (source_label(db_key), af_val)) elif disease_role == "unaffected": background_ac += ac_add background_an += an_add if cohort_observes: hits_background = True + if af_val is not None and af_val > 0 and not skip_top: + background_arm_afs.append( + (source_label(db_key), af_val)) else: background_ac += ac_add background_an += an_add if cohort_observes: hits_background = True + if af_val is not None and af_val > 0 and not skip_top: + background_arm_afs.append( + (source_label(db_key), af_val)) 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) @@ -488,87 +562,106 @@ 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": affected_ac += pop_ac_add affected_an += pop_an_add if pop_observes: hits_affected = True + if (pop_af_val is not None and pop_af_val > 0 + and not skip_top): + affected_arm_afs.append( + (source_label(db_key, pop["key"]), pop_af_val)) elif is_disease_db and pheno in ("unaffected", "unknown"): # Unaffected relatives, controls, and unknown-phenotype # individuals all feed the background. background_ac += pop_ac_add background_an += pop_an_add if pop_observes: hits_background = True + if (pop_af_val is not None and pop_af_val > 0 + and not skip_top): + background_arm_afs.append( + (source_label(db_key, pop["key"]), pop_af_val)) elif not is_disease_db: # 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. + # We also deliberately exclude these from the + # top-N source AFs: sub-populations like + # gnomAD-Finnish are not "sources", so they would + # crowd out actual project signals. 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 + # Top 3 source AFs for the mouseOver. Strings; empty if no source + # in the relevant group had AF > 0. + top_affected_sources = top_n_source_afs(affected_arm_afs, 3) + top_background_sources = top_n_source_afs(background_arm_afs, 3) + 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), + top_affected_sources, 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), + top_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") n_aff += 1 @@ -625,34 +718,36 @@ 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). # 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 topAffectedSources; "Top 3 affected/case cohorts by per-source AF"\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(' string topBackgroundSources; "Top 3 background cohorts by per-source AF"\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') f.write(')\n')