#!/usr/bin/env python3 """Analysis — rs4988235 (lactase persistence) population differentiation. Adapted from reference_analyses/population_allele_freq.py. Compares rs4988235's allele-frequency differentiation across 1000G super-populations against a matched neutral locus on the same chromosome. Input: --kg CSV from fetch_1000g.py (LCT/MCM6 region containing rs4988235) --kg-neu CSV from fetch_1000g.py (neutral region upstream) --target-pos GRCh37 position of rs4988235 (default: 136608646) Output: figure.png — grouped bar chart (rs4988235 AF by pop) + box plot (Fst comparison) stats.json — chi-square + Fst for rs4988235, region-wide Fst distributions, and comparison with neutral locus """ from __future__ import annotations import argparse import json import sys import matplotlib matplotlib.use("Agg") import matplotlib.pyplot as plt import numpy as np import pandas as pd from scipy import stats as sp_stats SUPERPOP_2N = {"afr": 1322, "amr": 694, "eas": 1008, "eur": 1006, "sas": 978} POPS = ["afr", "amr", "eas", "eur", "sas"] TARGET_POS = 136608646 # rs4988235 GRCh37 def fst_pair(p1, p2, n1, n2): if p1 in (0, 1) and p1 == p2: return 0.0 num = (p1 - p2) ** 2 - p1 * (1 - p1) / (n1 - 1) - p2 * (1 - p2) / (n2 - 1) den = p1 * (1 - p2) + p2 * (1 - p1) return float(num / den) if den > 0 else 0.0 def fst_matrix_for_variant(af_dict): fst = {} for i, a in enumerate(POPS): for b in POPS[i + 1:]: fst[f"{a}-{b}"] = round(fst_pair(af_dict[a], af_dict[b], SUPERPOP_2N[a], SUPERPOP_2N[b]), 4) return fst def region_max_fst(df): """Compute max-pairwise Fst for every variant in a region.""" afmat = df[[f"{p}_af" for p in POPS]].to_numpy(dtype=float) fsts = [] for row in afmat: mx = 0.0 for i, a in enumerate(POPS): for b in POPS[i + 1:]: v = fst_pair(row[i], row[POPS.index(b)], SUPERPOP_2N[a], SUPERPOP_2N[b]) if v > mx: mx = v fsts.append(mx) return np.array(fsts) def main() -> int: ap = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter) ap.add_argument("--kg", required=True, help="LCT/MCM6 region CSV") ap.add_argument("--kg-neu", required=True, help="Neutral region CSV") ap.add_argument("--target-pos", type=int, default=TARGET_POS, help="GRCh37 position of rs4988235") ap.add_argument("--seed", type=int, default=1234) ap.add_argument("--outdir", default=".") args = ap.parse_args() # --- Load and validate target region --- df = pd.read_csv(args.kg, dtype=str).fillna("") for p in POPS: df[f"{p}_af"] = pd.to_numeric(df.get(f"{p}_af"), errors="coerce") df = df.dropna(subset=[f"{p}_af" for p in POPS]) if df.empty: _fail(args.outdir, "no variants with all-population AF in target region") return 1 # --- Identify rs4988235 --- target = df[df["pos"] == str(args.target_pos)] if target.empty: _fail(args.outdir, f"rs4988235 not found at pos {args.target_pos}") return 1 target = target.iloc[0] target_af = {p: float(target[f"{p}_af"]) for p in POPS} # Chi-square for rs4988235 alt_counts = [round(target_af[p] * SUPERPOP_2N[p]) for p in POPS] ref_counts = [SUPERPOP_2N[p] - a for a in alt_counts] table = np.array([alt_counts, ref_counts]) if table.min() < 0 or table.sum() == 0: _fail(args.outdir, "degenerate allele-count table for rs4988235") return 1 chi2, p_value, dof, _ = sp_stats.chi2_contingency(table) fst_target = fst_matrix_for_variant(target_af) max_fst_pair = max(fst_target, key=fst_target.get) # --- Region-wide Fst distribution (LCT region) --- region_fsts = region_max_fst(df) # --- Load neutral region --- df_neu = pd.read_csv(args.kg_neu, dtype=str).fillna("") for p in POPS: df_neu[f"{p}_af"] = pd.to_numeric(df_neu.get(f"{p}_af"), errors="coerce") df_neu = df_neu.dropna(subset=[f"{p}_af" for p in POPS]) neu_fsts = region_max_fst(df_neu) # --- Compute rank of rs4988235 within target region --- target_fst = fst_target[max_fst_pair] rank = int(np.sum(region_fsts >= target_fst)) # 1-based: how many variants have Fst >= rs4988235 percentile = float(rank / len(region_fsts) * 100) # Mann-Whitney U comparing rs4988235's max-Fst to region distributions # One-sided: is rs4988235 more differentiated than random variants? # Compare to neutral region Fst distribution mw_stat, mw_p = sp_stats.mannwhitneyu( np.array([target_fst]), neu_fsts, alternative='greater' ) # --- Build result --- result = { "label": "rs4988235 (lactase persistence) vs neutral locus", "seed": args.seed, "reference_template": "population_allele_freq", "target_variant": f"rs4988235 at 2:{args.target_pos}", "target_af_by_pop": target_af, "n_variants_target": int(len(df)), "n_variants_neutral": int(len(df_neu)), "target_chi2": float(chi2), "target_dof": int(dof), "target_p_value": float(p_value), "target_pairwise_fst": fst_target, "target_max_fst_pair": max_fst_pair, "target_max_fst": float(fst_target[max_fst_pair]), "target_region_median_fst": float(np.median(region_fsts)), "target_region_mean_fst": float(np.mean(region_fsts)), "target_rank_in_region": rank, "target_percentile_in_region": percentile, "neutral_region_median_fst": float(np.median(neu_fsts)), "neutral_region_mean_fst": float(np.mean(neu_fsts)), "mannwhitney_u": float(mw_stat), "mannwhitney_p_greater": float(mw_p), "significant_at_0.05": bool(p_value < 0.05), "outcome": "success" if p_value < 0.05 else "null_result", "test": "chi-square of allele-count homogeneity across 5 super-populations (rs4988235)", "headline_statistic": f"chi2 p = {p_value:.2e}, max Fst = {fst_target[max_fst_pair]:.3f} ({max_fst_pair}), region rank {rank}/{len(region_fsts)}", } with open(f"{args.outdir}/stats.json", "w") as fh: json.dump(result, fh, indent=2) _render(args.outdir, target_af, result, region_fsts, neu_fsts, fst_target, max_fst_pair) print(json.dumps(result, indent=2)) return 0 def _fail(outdir, reason): result = {"label": "rs4988235", "reference_template": "population_allele_freq", "outcome": "failed", "reason": reason} with open(f"{outdir}/stats.json", "w") as fh: json.dump(result, fh, indent=2) print(json.dumps(result, indent=2)) def _render(outdir, target_af, result, region_fsts, neu_fsts, fst_target, max_fst_pair): fig, axes = plt.subplots(1, 2, figsize=(12, 5)) # Panel A: rs4988235 AF by super-population ax = axes[0] colors = ["#8e44ad", "#e67e22", "#16a085", "#2980b9", "#c0392b"] af_vals = [target_af[p] for p in POPS] bars = ax.bar([p.upper() for p in POPS], af_vals, color=colors, alpha=0.85, edgecolor="black", linewidth=0.5) ax.set_ylabel("Alt allele frequency") ax.set_title(f"rs4988235 allele frequencies\n" f"chi2 p = {result['target_p_value']:.2e}", fontsize=11) ax.set_ylim(0, max(af_vals) * 1.2 + 0.05) ax.grid(axis="y", alpha=0.3) for bar, val in zip(bars, af_vals): ax.text(bar.get_x() + bar.get_width() / 2, bar.get_height() + 0.01, f"{val:.3f}", ha="center", va="bottom", fontsize=9) # Panel B: Fst comparison (box plot) ax = axes[1] bp = ax.boxplot([neu_fsts, region_fsts], tick_labels=["Neutral\nregion", "LCT/MCM6\nregion"], patch_artist=True, widths=0.5) bp["boxes"][0].set_facecolor("#bdc3c7") bp["boxes"][1].set_facecolor("#3498db") for box in bp["boxes"]: box.set_alpha(0.7) # Overlay rs4988235's Fst as a red diamond ax.scatter([2], [result["target_max_fst"]], marker="D", color="#e74c3c", s=100, zorder=5, label=f"rs4988235 (Fst={result['target_max_fst']:.3f})") ax.set_ylabel("Max pairwise Fst") ax.set_title("Population differentiation\n" f"Mann-Whitney p = {result['mannwhitney_p_greater']:.2e}", fontsize=11) ax.legend(fontsize=9, loc="upper left") ax.grid(axis="y", alpha=0.3) fig.suptitle("Lactase persistence variant (rs4988235) population differentiation\n" f"{result['target_variant']} — 1000 Genomes phase 3", fontsize=13, fontweight="bold", y=1.02) fig.tight_layout() fig.savefig(f"{outdir}/figure.png", dpi=150, bbox_inches="tight") plt.close(fig) if __name__ == "__main__": sys.exit(main())