softwareDevelopment/GWAS
0
1import gradio as gr2import pandas as pd3import numpy as np4import matplotlib5matplotlib.use("Agg")6import matplotlib.pyplot as plt7import matplotlib.patches as mpatches8from matplotlib.colors import LinearSegmentedColormap9from matplotlib.lines import Line2D10from matplotlib.patches import Patch11import matplotlib.gridspec as gridspec12import matplotlib.patheffects as pe13from pycirclize import Circos14import warnings15warnings.filterwarnings("ignore")16 17from scipy import stats18from scipy.optimize import minimize_scalar19from scipy.cluster.hierarchy import dendrogram, linkage20from sklearn.decomposition import PCA21from sklearn.preprocessing import StandardScaler22from sklearn.linear_model import LinearRegression23from sklearn.manifold import MDS24from statsmodels.stats.multitest import multipletests25import io, re, os, zipfile, tempfile26from io import StringIO27from openpyxl import Workbook28from openpyxl.styles import Font, PatternFill, Alignment, Border, Side29from openpyxl.utils import get_column_letter30 31# ─────────────────────── DEFAULT COLOR PALETTES ──────────────────────────────────32BG_DARK = "#0D1117"33PANEL_DARK = "#161B22"34PANEL_MID = "#21262D"35ACCENT_TEAL = "#2DD4BF"36ACCENT_GOLD = "#F59E0B"37ACCENT_CORAL= "#F87171"38ACCENT_BLUE = "#60A5FA"39ACCENT_LIME = "#84CC16"40ACCENT_PURP = "#A78BFA"41ACCENT_PINK = "#F472B6"42ACCENT_CYAN = "#22D3EE"43ACCENT_ORG = "#FB923C"44TEXT_MAIN = "#E6EDF3"45TEXT_DIM = "#8B949E"46GRID_LINE = "#30363D"47SIG_COLOR = "#F59E0B"48 49# Light-theme equivalents -- used for ALL plots (PCA must NOT default to dark)50BG_LIGHT = "#FFFFFF"51PANEL_LIGHT = "#FAFBFF"52TEXT_DARK = "#2D3142"53GRID_LIGHT = "#E8ECF4"54ACCENT_TEAL_D = "#0F766E"55ACCENT_GOLD_D = "#B45309"56ACCENT_PINK_D = "#BE185D"57 58MODEL_COLORS = {59 "OLS": "#60A5FA",60 "MLM": "#2DD4BF",61 "EMMAX": "#F472B6",62 "FarmCPU": "#84CC16",63 "GEMMA": "#E879F9",64 "BLINK": "#FB923C",65 "mrMLM": "#A78BFA",66 "FASTmrMLM": "#F59E0B",67}68 69CHR_PALETTE = [70 "#4E79A7","#F28E2B","#E15759","#76B7B2","#59A14F",71 "#EDC948","#B07AA1","#FF9DA7","#9C755F","#BAB0AC",72 "#499894","#86BCB6","#D4A6C8","#FFBE7D","#8CD17D",73]74 75ALL_MODELS = ["OLS", "MLM", "EMMAX", "FarmCPU", "GEMMA", "BLINK", "mrMLM", "FASTmrMLM"]76 77def style_fig_white(fig, font_size=9, font_color="#2D3142", grid_color="#E8ECF4"):78 fig.patch.set_facecolor("#FFFFFF")79 for ax in fig.get_axes():80 ax.set_facecolor("#FAFBFF")81 ax.tick_params(colors=font_color, labelsize=font_size)82 ax.xaxis.label.set_color(font_color)83 ax.yaxis.label.set_color(font_color)84 ax.title.set_color(font_color)85 for spine in ax.spines.values():86 spine.set_edgecolor("#CCCCCC")87 ax.grid(color=grid_color, linewidth=0.5, linestyle="--", alpha=0.7)88 return fig89 90def apply_custom_style(fig, cfg):91 bg = cfg.get("bg_color", "#FFFFFF")92 panel = cfg.get("panel_color", "#FAFBFF")93 fc = cfg.get("font_color", "#2D3142")94 fs = cfg.get("font_size", 9)95 gc = cfg.get("grid_color", "#E8ECF4")96 fig.patch.set_facecolor(bg)97 for ax in fig.get_axes():98 ax.set_facecolor(panel)99 ax.tick_params(colors=fc, labelsize=fs)100 ax.xaxis.label.set_color(fc)101 ax.yaxis.label.set_color(fc)102 ax.title.set_color(fc)103 for spine in ax.spines.values():104 spine.set_edgecolor("#CCCCCC")105 ax.grid(color=gc, linewidth=0.5, linestyle="--", alpha=0.7)106 ax.xaxis.label.set_fontsize(fs)107 ax.yaxis.label.set_fontsize(fs)108 ax.title.set_fontsize(fs + 2)109 return fig110 111# ─────────────────────── DATA PARSING ────────────────────────────────────────────112def parse_data(file_obj=None, text_input=None):113 raw = None114 if file_obj is not None:115 try:116 raw = pd.read_csv(file_obj.name, sep=None, engine="python")117 except Exception:118 raw = pd.read_csv(file_obj.name, sep="\t")119 elif text_input and text_input.strip():120 try:121 raw = pd.read_csv(StringIO(text_input), sep=None, engine="python")122 except Exception:123 raw = pd.read_csv(StringIO(text_input), sep="\t")124 else:125 raise ValueError("No data provided.")126 127 chrom_info = {}128 if len(raw) > 1:129 second_row = raw.iloc[1].astype(str)130 chrom_vals = second_row.iloc[2:].astype(str)131 if len(chrom_vals) > 0 and chrom_vals.str.match(r'^\d+$').all():132 for i, col in enumerate(raw.columns[2:]):133 try:134 chrom_info[col] = int(second_row.iloc[2 + i])135 except:136 chrom_info[col] = 1137 raw = raw.drop(index=1).reset_index(drop=True)138 139 cols = list(raw.columns)140 id_col = cols[0]141 phen_col = cols[1]142 snp_cols = cols[2:]143 144 raw = raw[pd.to_numeric(raw[phen_col], errors="coerce").notna()].copy()145 raw[phen_col] = raw[phen_col].astype(float)146 geno = raw[snp_cols].apply(pd.to_numeric, errors="coerce")147 148 # Only these are unambiguous missing-value sentinels. -1 is deliberately149 # NOT in this list: many marker matrices (rrBLUP/GAPIT-style) use -1/0/1150 # as REAL genotype calls (homozygous ref / het / homozygous alt), not a151 # missing code. Blindly nuking -1 as missing in that case wipes out152 # ~30-55% of real genotype calls and corrupts MAF, kinship, and PCA153 # downstream (all of which assume a clean 0/1/2 dosage scale).154 geno = geno.where(~geno.isin([-9, -99, -999]), np.nan)155 156 # Detect the coding scheme from the observed value range and normalize157 # everything to a 0/1/2 additive dosage scale so downstream math is158 # always consistent regardless of which format the input used.159 obs_vals = geno.values[~pd.isna(geno.values)]160 if obs_vals.size > 0 and obs_vals.min() >= -1 and obs_vals.max() <= 1:161 # -1/0/1 coding detected -> rescale to 0/1/2162 geno = geno + 1163 else:164 # 0/1/2 coding: here -1 really is a missing-value sentinel165 geno = geno.where(geno != -1, np.nan)166 167 df = pd.DataFrame()168 df["ID"] = raw[id_col].values169 df["Phenotype"] = raw[phen_col].values170 for c in snp_cols:171 df[c] = geno[c].values172 173 chroms, positions = [], []174 # Fallback position counters MUST be per-chromosome, not a single running175 # index over the whole marker file. A global counter (the old176 # `len(positions) + 1`) makes later chromosomes carry huge positional177 # offsets (e.g. chr9 markers numbered ~3100-3400 instead of 1-300), which178 # plot_manhattan() masks by subtracting each chromosome's own min before179 # plotting, but plot_circos_density() does not -- producing a genuinely180 # empty arc at the start of every sector sized to the wrong (inflated)181 # chromosome length. Keeping the counter per-chromosome fixes it at the182 # source instead of relying on every downstream plot to compensate.183 chr_counters = {}184 185 def _next_idx(c):186 chr_counters[c] = chr_counters.get(c, 0) + 1187 return chr_counters[c]188 189 for m in snp_cols:190 if m in chrom_info:191 c = chrom_info[m]192 chroms.append(c)193 positions.append(_next_idx(c))194 else:195 m_clean = str(m).strip()196 mat = re.match(r"[Cc]hr(\d+)[_\-](\d+)", m_clean)197 if mat:198 c = int(mat.group(1))199 chroms.append(c)200 positions.append(int(mat.group(2)))201 else:202 num = re.findall(r"\d+", m_clean)203 c = int(num[0]) if num else 1204 chroms.append(c)205 positions.append(int(num[1]) if len(num) > 1 else _next_idx(c))206 207 return df, snp_cols, np.array(chroms), np.array(positions)208 209# ─────────────────────── UTILITY ─────────────────────────────────────────────────210def impute_genotypes(geno_matrix):211 mat = geno_matrix.copy().astype(float)212 for j in range(mat.shape[1]):213 col = mat[:, j]214 m = np.nanmean(col)215 col[np.isnan(col)] = m if not np.isnan(m) else 0216 mat[:, j] = col217 return mat218 219def compute_maf(G):220 freqs = np.nanmean(G, axis=0) / 2221 return np.where(freqs > 0.5, 1 - freqs, freqs)222 223def compute_call_rate(G):224 return 1 - np.isnan(G).mean(axis=0)225 226def build_kinship(G):227 p = np.nanmean(G, axis=0) / 2228 Z = G - 2 * p229 denom = 2 * np.sum(p * (1 - p))230 if denom == 0:231 denom = 1232 return Z @ Z.T / denom233 234def compute_pve_all(y, G, K=None, X0=None):235 """236 PVE (%) for every marker via genuine variance-component partitioning237 instead of a raw squared Pearson correlation.238 239 A plain r^2 between one marker and y only equals a valid PVE estimate240 when the model has no other covariates and no relatedness structure241 (it happens to coincide with SS_marker/SS_total for a bare242 single-predictor OLS fit). Once a kinship random effect or Q243 covariates are in the model, that raw correlation double-counts244 variance already absorbed by relatedness/structure and no longer245 matches what GAPIT/GEMMA-style PVE figures report.246 247 This instead computes, for each marker j:248 PVE_j = (RSS_null - RSS_full_j) / SS_total(y) * 100249 where RSS_null is the residual sum of squares of the model with just250 the intercept (+ X0 covariates, if given), RSS_full_j additionally251 includes marker j, and SS_total(y) is the total phenotypic sum of252 squares -- the same denominator convention GREML/GCTA-style "variance253 explained" figures use, so PVE values from different markers are254 directly comparable and summable, unlike per-marker r^2.255 256 If K (kinship) is supplied, y/X0/G are first weighted by the257 REML-estimated genetic/residual variance ratio (delta), exactly as258 mixed_model_gwas does, so a marker's PVE is *conditional on* the259 polygenic background already explained by relatedness rather than260 re-claiming variance the random effect already accounts for.261 """262 y = np.asarray(y, dtype=float)263 n, m = G.shape264 if X0 is None:265 X0 = np.ones((n, 1))266 267 SS_total = float(np.sum((y - y.mean()) ** 2))268 if SS_total <= 0:269 return np.zeros(m)270 271 if K is not None:272 try:273 evals, evecs = eigh_kinship(K)274 _, _, delta, _, _ = emma_reml(y, K, X0, evals=evals, evecs=evecs)275 w = 1.0 / (evals + delta)276 sw = np.sqrt(w)277 yt = (evecs.T @ y) * sw278 X0t = (evecs.T @ X0) * sw[:, None]279 Gt = (evecs.T @ G) * sw[:, None]280 except Exception:281 yt, X0t, Gt = y, X0, G282 else:283 yt, X0t, Gt = y, X0, G284 285 beta0, *_ = np.linalg.lstsq(X0t, yt, rcond=None)286 rss0 = float(np.sum((yt - X0t @ beta0) ** 2))287 288 pve = np.zeros(m)289 for j in range(m):290 xj = Gt[:, j]291 if np.std(xj) < 1e-8:292 continue293 Xf = np.column_stack([X0t, xj])294 try:295 beta, *_ = np.linalg.lstsq(Xf, yt, rcond=None)296 except np.linalg.LinAlgError:297 continue298 rss = float(np.sum((yt - Xf @ beta) ** 2))299 pve[j] = max(0.0, rss0 - rss) / SS_total * 100300 return pve301 302 303RENDER_DPI = 300 # publication-quality; was 150304EXPORT_DPI = 400 # all downloadable PNG/figure files must exceed 300 dpi305 306def fig_to_png_bytes(fig):307 buf = io.BytesIO()308 fig.savefig(buf, format="png", dpi=RENDER_DPI, bbox_inches="tight", facecolor=fig.get_facecolor())309 buf.seek(0)310 return buf.read()311 312def fig_to_pil(fig):313 buf = io.BytesIO()314 fig.savefig(buf, format="png", dpi=RENDER_DPI, bbox_inches="tight", facecolor=fig.get_facecolor())315 buf.seek(0)316 from PIL import Image317 return Image.open(buf).copy()318 319def get_top_mta(pvals, n=15):320 return np.argsort(pvals)[:min(n, len(pvals))]321 322def get_significant_idx(pvals, lod_threshold=None, fallback_top_n=10,323 pve_vals=None, pve_threshold=None):324 """325 Indices of markers that actually PASS the significance threshold --326 this is what should be highlighted as MTAs, not a blind top-N slice327 that highlights N markers regardless of whether any of them are real.328 329 lod_threshold is on the LOD scale (see neglog10p_to_lod) and is meant330 to be the permutation-derived empirical threshold -- the primary331 significance criterion for this pipeline, the same role a LOD332 threshold plays in linkage/QTL mapping. Returns indices sorted by333 increasing p-value (most significant first). If no permutation334 threshold is available (e.g. it failed to compute), falls back to a335 small top-N so plots don't render silently empty -- this is a safety336 net, not the intended mode of operation.337 338 pve_vals / pve_threshold add a SECOND, independent filter: if both are339 given, a marker must ALSO explain at least pve_threshold % of340 phenotypic variance to be kept as an MTA, on top of passing the LOD341 threshold. pve_threshold of None/0 disables this filter.342 """343 pvals = np.asarray(pvals)344 if lod_threshold is None:345 idx = get_top_mta(pvals, fallback_top_n)346 else:347 log_p = -np.log10(np.clip(pvals, 1e-300, 1))348 lod_vals = neglog10p_to_lod(log_p)349 idx = np.where(lod_vals >= lod_threshold)[0]350 idx = idx[np.argsort(pvals[idx])]351 352 if pve_vals is not None and pve_threshold:353 pve_vals = np.asarray(pve_vals)354 idx = idx[pve_vals[idx] >= pve_threshold]355 356 return idx357 358def compute_thresholds(n_snps, model_name=""):359 """360 Compute and explain Bonferroni and suggestive thresholds.361 Returns thresholds and explanation text.362 """363 bonferroni = 0.05 / n_snps364 suggestive = 1.0 / n_snps # 1 expected false positive genome-wide365 366 model_notes = {367 "OLS": "No structure/relatedness correction — baseline; expect inflation (λ > 1).",368 "MLM": "Q (kinship PCs, fixed) + K (polygenic random effect, REML) — Yu et al. 2006.",369 "EMMAX": "Pure kinship P3D mixed model; variance components estimated once genome-wide.",370 "FarmCPU": "Iterative FEM/REM with bin-based pseudo-QTN cofactors (Liu et al. 2016).",371 "GEMMA": "Exact per-marker REML LMM (Zhou & Stephens 2012); no post-hoc GC rescaling.",372 }373 note = model_notes.get(model_name, "")374 375 explanation = (376 f"──────────────────────────────────────────────\n"377 f" Threshold Calculations — {model_name}\n"378 f"──────────────────────────────────────────────\n"379 f" SNPs tested (m) : {n_snps:,}\n\n"380 f" 📌 Bonferroni threshold:\n"381 f" α / m = 0.05 / {n_snps:,} = {bonferroni:.3e}\n"382 f" -log₁₀(p) = {-np.log10(bonferroni):.2f}\n"383 f" Controls genome-wide Type-I error at 5%.\n"384 f" Assumes all tests are independent.\n\n"385 f" 📌 Suggestive threshold:\n"386 f" 1 / m = 1 / {n_snps:,} = {suggestive:.3e}\n"387 f" -log₁₀(p) = {-np.log10(suggestive):.2f}\n"388 f" Expects ~1 false positive per genome scan.\n"389 f" Lander & Kruglyak (1995) convention.\n\n"390 f" ℹ️ Model note: {note}\n"391 f"──────────────────────────────────────────────\n"392 )393 return bonferroni, suggestive, explanation394 395# ─────────────────────── EMMA / REML MIXED-MODEL ENGINE ──────────────────────────396# Shared machinery for MLM / EMMAX / GEMMA, following the eigendecomposition397# approach of Kang et al. 2008 (EMMA) and its extensions (EMMAX, FaST-LMM,398# GEMMA). The kinship matrix is decomposed ONCE; every subsequent variance-399# component search reuses the eigenvalues, which is what makes these methods400# tractable at genome scale instead of refitting an n x n GLS from scratch401# per marker.402 403def eigh_kinship(K):404 """Symmetric eigendecomposition of the kinship matrix. Eigenvalues are405 floored at a small positive number to keep downstream log/inversion406 steps numerically stable (K can be near-singular for small panels)."""407 n = K.shape[0]408 Ks = (K + K.T) / 2.0409 evals, evecs = np.linalg.eigh(Ks)410 evals = np.clip(evals, 1e-8, None)411 return evals, evecs412 413def _profile_loglik(log_delta, evals, yt, Xt, reml=True):414 """REML (or ML) profile log-likelihood as a function of delta =415 sigma_e^2 / sigma_g^2, evaluated in the rotated (eigenvector) basis416 where the covariance structure becomes diagonal: var(yt_i) = D_i + delta.417 """418 delta = np.exp(log_delta)419 w = 1.0 / (evals + delta)420 sw = np.sqrt(w)421 Xw = Xt * sw[:, None]422 yw = yt * sw423 beta, *_ = np.linalg.lstsq(Xw, yw, rcond=None)424 resid = yw - Xw @ beta425 n = len(yt)426 q = Xt.shape[1]427 df = max(n - q, 1) if reml else n428 rss = float(np.sum(resid ** 2))429 sigma2 = rss / df if df > 0 else rss / max(n, 1)430 if sigma2 <= 0:431 sigma2 = 1e-12432 ll = -0.5 * (df * np.log(2 * np.pi * sigma2) + np.sum(np.log(evals + delta)) + df)433 if reml:434 XtWX = Xw.T @ Xw435 sign, logdet = np.linalg.slogdet(XtWX)436 if sign > 0:437 ll -= 0.5 * logdet438 return ll439 440def emma_reml(y, K, X=None, evals=None, evecs=None):441 """442 Estimate the genetic/residual variance ratio delta via REML, using the443 eigendecomposition of K (Kang et al. 2008). Returns eigvals, eigvecs,444 delta_hat, sigma_g2_hat, sigma_e2_hat.445 """446 n = len(y)447 if X is None:448 X = np.ones((n, 1))449 if evals is None or evecs is None:450 evals, evecs = eigh_kinship(K)451 yt = evecs.T @ y452 Xt = evecs.T @ X453 454 res = minimize_scalar(lambda ld: -_profile_loglik(ld, evals, yt, Xt, reml=True),455 bounds=(-10, 10), method="bounded",456 options={"xatol": 1e-4})457 delta_hat = float(np.exp(res.x))458 459 w = 1.0 / (evals + delta_hat)460 sw = np.sqrt(w)461 Xw = Xt * sw[:, None]462 yw = yt * sw463 beta, *_ = np.linalg.lstsq(Xw, yw, rcond=None)464 resid = yw - Xw @ beta465 df = max(n - X.shape[1], 1)466 sigma_g2 = float(np.sum(resid ** 2) / df)467 sigma_e2 = sigma_g2 * delta_hat468 return evals, evecs, delta_hat, sigma_g2, sigma_e2469 470def _wald_test_rotated(y_col, X_cols, evals, delta):471 """GLS Wald test for the last column of X_cols (the marker), with y and472 X already rotated into the eigenvector basis. Weights 1/(D_i+delta) make473 this an ordinary weighted least squares problem — the whole reason the474 eigendecomposition trick works."""475 w = 1.0 / (evals + delta)476 sw = np.sqrt(w)477 Xw = X_cols * sw[:, None]478 yw = y_col * sw479 try:480 XtX_inv = np.linalg.inv(Xw.T @ Xw)481 except np.linalg.LinAlgError:482 return 1.0, 0.0, np.nan483 beta = XtX_inv @ (Xw.T @ yw)484 resid = yw - Xw @ beta485 dfree = max(len(y_col) - X_cols.shape[1], 1)486 sigma2 = float(np.sum(resid ** 2) / dfree)487 se = float(np.sqrt(max(sigma2 * XtX_inv[-1, -1], 0)))488 b = float(beta[-1])489 t_stat = b / (se + 1e-300)490 p = float(2 * stats.t.sf(abs(t_stat), df=dfree))491 return p, b, se492 493def mixed_model_gwas(y, G, K, covariates=None):494 """495 P3D mixed-model GWAS (Population Parameters Previously Determined):496 variance components are estimated ONCE from the null model (no marker497 effect), then reused for every marker's GLS test. This is what real498 MLM/EMMAX software does by default (TASSEL, GAPIT, EMMAX itself) since499 per-marker REML gives negligible accuracy gain for a large speed cost.500 """501 n, m = G.shape502 X0 = np.ones((n, 1))503 if covariates is not None and covariates.shape[1] > 0:504 X0 = np.column_stack([X0, covariates])505 506 evals, evecs, delta, sigma_g2, sigma_e2 = emma_reml(y, K, X0)507 yt = evecs.T @ y508 X0t = evecs.T @ X0509 Gt = evecs.T @ G # rotate all markers once — O(n^2 m), not O(n^3) per marker510 511 pvals = np.ones(m); betas = np.zeros(m); ses = np.zeros(m)512 for j in range(m):513 if G[:, j].std() < 1e-8:514 continue515 Xcols = np.column_stack([X0t, Gt[:, j]])516 p, b, se = _wald_test_rotated(yt, Xcols, evals, delta)517 pvals[j] = p; betas[j] = b; ses[j] = se518 return pvals, betas, ses, delta519 520# ─────────────────────── PERMUTATION (LOD) THRESHOLD ──────────────────────────────521def permutation_threshold_lod(y, G, K, extra_covariates=None, n_perm=100, alpha=0.05, seed=42):522 """523 Empirical permutation significance threshold (Churchill & Doerge 1994524 style -- the same logic used for LOD thresholds in linkage/QTL mapping),525 expressed on the LOD scale via neglog10p_to_lod so it's directly526 comparable to a QTL-mapping LOD threshold. This is meant to be THE527 primary significance criterion for this pipeline -- which MTAs get528 highlighted should be driven by whether they pass this threshold, not529 by an arbitrary top-N slice.530 531 Phenotype is shuffled n_perm times; for each permutation the maximum532 -log10(p) across the whole genome is recorded via the fast P3D mixed-533 model engine (variance components estimated once on the REAL data and534 held fixed across permutations -- standard practice for tractable535 mixed-model permutation, since re-running REML per permutation per536 marker would be computationally prohibitive). The threshold is the537 (1-alpha) quantile of that null max-statistic distribution, converted538 to LOD units.539 """540 rng = np.random.default_rng(seed)541 n, m = G.shape542 X0 = np.ones((n, 1))543 if extra_covariates is not None and extra_covariates.shape[1] > 0:544 X0 = np.column_stack([X0, extra_covariates])545 546 evals, evecs, delta, _, _ = emma_reml(y, K, X0)547 X0t = evecs.T @ X0548 Gt = evecs.T @ G549 valid_markers = np.where(G.std(axis=0) > 1e-8)[0]550 551 max_lp = np.zeros(n_perm)552 for b in range(n_perm):553 y_perm = rng.permutation(y)554 yt = evecs.T @ y_perm555 best = 0.0556 for j in valid_markers:557 Xcols = np.column_stack([X0t, Gt[:, j]])558 p, _, _ = _wald_test_rotated(yt, Xcols, evals, delta)559 lp = -np.log10(max(p, 1e-300))560 if lp > best:561 best = lp562 max_lp[b] = best563 564 threshold_neglog10p = float(np.quantile(max_lp, 1 - alpha))565 threshold_lod = float(neglog10p_to_lod(threshold_neglog10p))566 return threshold_lod, threshold_neglog10p, max_lp567 568 569def ols_gwas(y, G):570 """Simple linear regression per marker, no correction for structure or571 relatedness. The baseline every other model is judged against."""572 n, m = G.shape573 pvals = np.ones(m); betas = np.zeros(m); ses = np.zeros(m)574 for j in range(m):575 x = G[:, j] - G[:, j].mean()576 if x.std() < 1e-8: continue577 slope, _, _, p, se = stats.linregress(x, y)578 pvals[j] = p; betas[j] = slope; ses[j] = se579 return pvals, betas, ses580 581def get_kinship_pcs(K, n_pca, n_samples):582 """Top n_pca eigenvectors of the kinship matrix, used as fixed-effect583 'Q' structure covariates by MLM, mrMLM, FASTmrMLM, and (optionally)584 GEMMA. Shared here so the per-model PCA covariate count set in the UI585 (Cofactors tab) is applied identically across every model that uses it."""586 n_pca = max(min(int(n_pca), n_samples - 4), 0)587 if n_pca <= 0:588 return None589 evals_K, evecs_K = eigh_kinship(K)590 order = np.argsort(evals_K)[::-1]591 return evecs_K[:, order[:n_pca]]592 593def mlm_gwas(y, G, K, n_pca=3):594 """Q + K mixed linear model (Yu et al. 2006): top kinship eigenvectors595 as fixed 'Q' structure covariates, PLUS the polygenic random effect596 (via REML-estimated variance components) rather than eigenvectors597 substituting for it. This is the actual unified MLM, not a PCA-only598 correction. n_pca sets how many kinship PCs are used as covariates."""599 n = len(y)600 Q = get_kinship_pcs(K, n_pca, n)601 pvals, betas, ses, _ = mixed_model_gwas(y, G, K, covariates=Q)602 return pvals, betas, ses603 604def emmax_gwas(y, G, K):605 """EMMAX (Kang et al. 2008): pure kinship-based P3D mixed model, no606 separate Q structure covariates — relatedness alone absorbs population607 structure. Same REML/eigendecomposition engine as MLM, without the608 fixed PCs."""609 pvals, betas, ses, _ = mixed_model_gwas(y, G, K, covariates=None)610 return pvals, betas, ses611 612def _rem_prune_qtns(y, G, K, qtn_idx, alpha=0.01):613 """614 REM (Random Effect Model) step of FarmCPU: re-test each candidate615 pseudo-QTN's own effect in a kinship-based mixed model, with the other616 current pseudo-QTNs held in as fixed covariates, and drop any QTN that617 is no longer significant once the polygenic random effect is in the618 model. Real FarmCPU alternates exactly this -- a REM re-optimization of619 the pseudo-QTN set -- with the FEM genome scan; without it (fixed620 effects only, as in the previous version of this function) pseudo-QTNs621 that are only "significant" because of relatedness/background genetic622 correlation, rather than a real local effect, never get removed.623 """624 if not qtn_idx:625 return qtn_idx626 n = len(y)627 keep = []628 for idx in qtn_idx:629 others = [q for q in qtn_idx if q != idx]630 cov = None631 if others:632 cov = G[:, others] - G[:, others].mean(axis=0)633 try:634 p, _, _, _ = mixed_model_gwas(y, G[:, [idx]], K, covariates=cov)635 except Exception:636 keep.append(idx) # fit failed -- don't silently drop a candidate637 continue638 if p[0] <= alpha:639 keep.append(idx)640 return keep641 642def farmcpu_gwas(y, G, K, chroms=None, positions=None, max_iter=8, bin_size=1_000_000,643 rem_alpha=0.01):644 """645 FarmCPU (Liu et al. 2016): iterates (1) a bin-based candidate scan --646 a Fixed Effect Model (FEM) test of every marker with the current647 pseudo-QTNs as covariates, followed by one candidate kept per bin --648 and (2) a Random Effect Model (REM) step that re-tests each candidate649 pseudo-QTN in a kinship-based mixed model and drops any that aren't650 significant once relatedness is accounted for. This is the actual651 FEM/REM alternation FarmCPU is defined by, not FEM alone: the REM step652 is what keeps pseudo-QTNs that are only correlated with y through653 background relatedness (rather than a true local effect) from being654 locked in as covariates for the rest of the run.655 """656 n, m = G.shape657 if positions is None:658 positions = np.arange(m)659 if chroms is None:660 chroms = np.zeros(m, dtype=int)661 662 pvals = np.ones(m); betas = np.zeros(m); ses = np.zeros(m)663 pseudo_qtns, prev_qtns = [], None664 665 for _ in range(max_iter):666 X = np.ones((n, 1))667 if pseudo_qtns:668 cov = G[:, pseudo_qtns] - G[:, pseudo_qtns].mean(axis=0)669 X = np.column_stack([X, cov])670 671 pvals = np.ones(m); betas = np.zeros(m); ses = np.zeros(m)672 qtn_set = set(pseudo_qtns)673 for j in range(m):674 x = G[:, j]675 if x.std() < 1e-8:676 continue677 # A marker that is currently a pseudo-QTN is tested leave-one-out678 # (covariates = the OTHER pseudo-QTNs), not skipped -- otherwise679 # the true causal marker, which is exactly the one most likely680 # to have been selected as its own covariate, ends up at p=1.681 if j in qtn_set:682 others = [q for q in pseudo_qtns if q != j]683 if others:684 cov = G[:, others] - G[:, others].mean(axis=0)685 Xj = np.column_stack([np.ones((n, 1)), cov])686 else:687 Xj = np.ones((n, 1))688 else:689 Xj = X690 Xf = np.column_stack([Xj, x])691 try:692 beta = np.linalg.lstsq(Xf, y, rcond=None)[0]693 XtX_inv = np.linalg.pinv(Xf.T @ Xf)694 except np.linalg.LinAlgError:695 continue696 resid = y - Xf @ beta697 dfree = max(n - Xf.shape[1], 1)698 sigma2 = float(np.sum(resid ** 2) / dfree)699 se = float(np.sqrt(max(sigma2 * XtX_inv[-1, -1], 0)))700 b = float(beta[-1])701 t_stat = b / (se + 1e-300)702 pvals[j] = float(2 * stats.t.sf(abs(t_stat), df=dfree))703 betas[j] = b; ses[j] = se704 705 bonf = 0.05 / m706 cand = np.where(pvals < bonf)[0]707 if len(cand) == 0:708 cand = np.argsort(pvals)[:5]709 710 bins = {}711 for idx in cand:712 key = (int(chroms[idx]), int(positions[idx] // bin_size))713 if key not in bins or pvals[idx] < pvals[bins[key]]:714 bins[key] = idx715 new_qtns = sorted(bins.values(), key=lambda i: pvals[i])[:10]716 717 # REM step: prune bin-selected candidates that don't hold up once718 # tested in a kinship random-effect model (see _rem_prune_qtns).719 new_qtns = _rem_prune_qtns(y, G, K, new_qtns, alpha=rem_alpha)720 721 if prev_qtns is not None and set(new_qtns) == set(prev_qtns):722 pseudo_qtns = new_qtns723 break724 prev_qtns = new_qtns725 pseudo_qtns = new_qtns726 727 return pvals, betas, ses728 729def _fixed_effect_scan(y, G, cov_idx, test_leave_one_out=True):730 """One fixed-effect-model scan of every marker in G, with the markers731 listed in cov_idx included as covariates (removed from the covariate732 set for a marker that IS one of them, i.e. tested leave-one-out).733 Shared by BLINK / mrMLM / FASTmrMLM below so their iteration logic734 stays readable."""735 n, m = G.shape736 pvals = np.ones(m); betas = np.zeros(m); ses = np.zeros(m)737 cov_set = set(cov_idx)738 for j in range(m):739 x = G[:, j]740 if x.std() < 1e-8:741 continue742 if test_leave_one_out and j in cov_set:743 others = [q for q in cov_idx if q != j]744 else:745 others = list(cov_idx)746 if others:747 cov = G[:, others] - G[:, others].mean(axis=0)748 Xj = np.column_stack([np.ones((n, 1)), cov])749 else:750 Xj = np.ones((n, 1))751 Xf = np.column_stack([Xj, x])752 try:753 beta = np.linalg.lstsq(Xf, y, rcond=None)[0]754 XtX_inv = np.linalg.pinv(Xf.T @ Xf)755 except np.linalg.LinAlgError:756 continue757 resid = y - Xf @ beta758 dfree = max(n - Xf.shape[1], 1)759 sigma2 = float(np.sum(resid ** 2) / dfree)760 se = float(np.sqrt(max(sigma2 * XtX_inv[-1, -1], 0)))761 b = float(beta[-1])762 t_stat = b / (se + 1e-300)763 pvals[j] = float(2 * stats.t.sf(abs(t_stat), df=dfree))764 betas[j] = b; ses[j] = se765 return pvals, betas, ses766 767def _bic_fixed_effect_model(y, G, cand_idx):768 """Fit y ~ intercept + G[:, cand_idx] and return BIC, used by BLINK to769 pick the best-supported pseudo-QTN subset instead of a fixed top-N."""770 n = len(y)771 if len(cand_idx) == 0:772 X = np.ones((n, 1))773 else:774 X = np.column_stack([np.ones((n, 1)), G[:, cand_idx]])775 beta, *_ = np.linalg.lstsq(X, y, rcond=None)776 resid = y - X @ beta777 k = X.shape[1]778 rss = float(np.sum(resid ** 2)) + 1e-300779 return n * np.log(rss / n) + k * np.log(n)780 781def blink_gwas(y, G, K, chroms=None, positions=None, max_iter=8, ld_r2=0.7,782 ld_block_bp=1_000_000):783 """784 BLINK (Huang et al. 2019, Bayesian-information and Linkage-disequilibrium785 Iteratively Nested Keyway): like FarmCPU, alternates a fixed-effect-model786 scan with pseudo-QTN selection, but differs in two places that give787 BLINK its name/speed advantage over FarmCPU: (1) candidate pseudo-QTNs788 are pruned within local LD blocks (candidates on the same chromosome and789 within ld_block_bp of an already-kept, more-significant candidate are790 compared by pairwise r^2 and dropped if redundant; candidates outside791 that window are never compared and are always kept, since LD is a local792 phenomenon and two markers hundreds of kb apart on different LD blocks793 can be "collapsed" by chance correlation in a small panel even with no794 real linkage) instead of a fixed physical bin size, so two markers 50bp795 apart in high LD are collapsed the same way two markers 500kb apart in796 low LD are kept separate; (2) the number of pseudo-QTNs carried into the797 next iteration is chosen by minimising BIC over nested candidate798 subsets, not a fixed top-10 cutoff. No kinship matrix is used once QTNs799 are found — BLINK, like FarmCPU, assumes the pseudo-QTNs themselves800 absorb population structure/relatedness.801 """802 n, m = G.shape803 if positions is None:804 positions = np.arange(m)805 if chroms is None:806 chroms = np.zeros(m, dtype=int)807 808 pseudo_qtns, prev_qtns = [], None809 pvals = np.ones(m); betas = np.zeros(m); ses = np.zeros(m)810 811 for _ in range(max_iter):812 pvals, betas, ses = _fixed_effect_scan(y, G, pseudo_qtns)813 814 bonf = 0.05 / m815 cand = np.where(pvals < bonf)[0]816 if len(cand) == 0:817 cand = np.argsort(pvals)[:10]818 cand = cand[np.argsort(pvals[cand])]819 820 # LD-prune candidates: keep a candidate only if it isn't in high LD821 # with an already-kept, more-significant candidate FROM THE SAME LD822 # BLOCK (same chromosome, within ld_block_bp). Candidates outside823 # that window skip the correlation check entirely -- LD blocks are824 # local, so comparing markers genome-wide (the previous behaviour)825 # could collapse two unlinked-but-coincidentally-correlated distant826 # markers into one pseudo-QTN.827 kept = []828 for idx in cand:829 redundant = False830 for k_idx in kept:831 same_block = (chroms[idx] == chroms[k_idx] and832 abs(int(positions[idx]) - int(positions[k_idx])) <= ld_block_bp)833 if not same_block:834 continue835 r = np.corrcoef(G[:, idx], G[:, k_idx])[0, 1]836 if np.isfinite(r) and r ** 2 >= ld_r2:837 redundant = True838 break839 if not redundant:840 kept.append(idx)841 if len(kept) >= 15:842 break843 844 # BIC-based model selection: try nested prefixes of `kept`845 # (most-significant-first) and keep the prefix with lowest BIC846 best_bic, best_set = np.inf, []847 for size in range(0, len(kept) + 1):848 subset = kept[:size]849 bic = _bic_fixed_effect_model(y, G, subset)850 if bic < best_bic:851 best_bic, best_set = bic, subset852 new_qtns = best_set853 854 if prev_qtns is not None and set(new_qtns) == set(prev_qtns):855 pseudo_qtns = new_qtns856 break857 prev_qtns = new_qtns858 pseudo_qtns = new_qtns859 860 return pvals, betas, ses861 862def mrmlm_gwas(y, G, K, chroms=None, positions=None, screen_thresh=0.01,863 max_candidates=20, drop_thresh=0.05, n_pca=3):864 """865 mrMLM (multi-locus random-SNP-effect mixed linear model; Wang et al.866 2016): a two-step multi-locus approach. Step 1 runs a lenient867 single-locus P3D mixed-model scan (same engine as MLM/EMMAX) and keeps868 every marker below `screen_thresh` as a candidate QTN. Step 2 fits all869 candidates simultaneously as fixed effects (approximating mrMLM's870 random-SNP-effect / EM-Bayes shrinkage step) and backward-eliminates871 the least significant one at a time until every remaining candidate is872 significant at `drop_thresh` — mrMLM's real innovation is testing873 QTNs jointly rather than one-at-a-time, which is what's reproduced874 here. Markers never selected as candidates keep their Step-1 p-value875 so the full-length p-value vector still plots as a Manhattan track.876 """877 n, m = G.shape878 Q = get_kinship_pcs(K, n_pca, n)879 pvals, betas, ses, _ = mixed_model_gwas(y, G, K, covariates=Q)880 881 cand = np.where(pvals < screen_thresh)[0]882 if len(cand) == 0:883 return pvals, betas, ses884 cand = cand[np.argsort(pvals[cand])][:max_candidates].tolist()885 886 # backward elimination on the joint multi-locus fixed-effect model887 while len(cand) > 0:888 p_joint, b_joint, se_joint = _fixed_effect_scan(y, G, cand)889 cand_p = {c: p_joint[c] for c in cand}890 worst = max(cand_p, key=cand_p.get)891 if cand_p[worst] > drop_thresh and len(cand) > 1:892 cand.remove(worst)893 continue894 break895 896 p_final, b_final, se_final = _fixed_effect_scan(y, G, cand)897 for c in cand:898 pvals[c] = p_final[c]; betas[c] = b_final[c]; ses[c] = se_final[c]899 return pvals, betas, ses900 901def fastmrmlm_gwas(y, G, K, chroms=None, positions=None, screen_thresh=0.01,902 max_candidates=20, drop_thresh=0.05, n_pca=3):903 """904 FASTmrMLM (Tamba et al. 2017): the speed-optimised successor to mrMLM.905 The real method reuses a single kinship eigen-decomposition (computed906 once) to transform y and G, avoiding a per-marker REML fit; that's907 exactly what EMMAX/mixed_model_gwas already does here, so Step 1908 candidate screening below IS that fast rotated scan. Step 2 departs909 from mrMLM by fitting the joint multi-locus model as plain OLS910 (dropping the kinship term for the joint step, trading a small amount911 of type-I-error control for the speed gain that gives the method its912 name) with the same backward-elimination logic as mrMLM.913 """914 n, m = G.shape915 Q = get_kinship_pcs(K, n_pca, n)916 pvals, betas, ses, _ = mixed_model_gwas(y, G, K, covariates=Q)917 918 cand = np.where(pvals < screen_thresh)[0]919 if len(cand) == 0:920 return pvals, betas, ses921 cand = cand[np.argsort(pvals[cand])][:max_candidates].tolist()922 923 while len(cand) > 0:924 p_joint, b_joint, se_joint = _fixed_effect_scan(y, G, cand,925 test_leave_one_out=True)926 cand_p = {c: p_joint[c] for c in cand}927 worst = max(cand_p, key=cand_p.get)928 if cand_p[worst] > drop_thresh and len(cand) > 1:929 cand.remove(worst)930 continue931 break932 933 p_final, b_final, se_final = _fixed_effect_scan(y, G, cand,934 test_leave_one_out=True)935 for c in cand:936 pvals[c] = p_final[c]; betas[c] = b_final[c]; ses[c] = se_final[c]937 return pvals, betas, ses938 939def gemma_gwas(y, G, K, covariates=None):940 """941 Approximate GEMMA-style univariate LMM (Zhou & Stephens 2012): unlike942 the P3D models above, the genetic/residual variance ratio is943 re-estimated PER MARKER (the candidate SNP is included in the null944 model before REML), which is GEMMA's default exact behaviour. Slower945 than P3D but the most statistically rigorous of the five — no post-hoc946 p-value rescaling is applied, since a correctly specified LMM shouldn't947 need genomic-control patching.948 """949 n, m = G.shape950 X0 = np.ones((n, 1))951 if covariates is not None and covariates.shape[1] > 0:952 X0 = np.column_stack([X0, covariates])953 evals, evecs = eigh_kinship(K)954 yt = evecs.T @ y955 X0t = evecs.T @ X0956 Gt = evecs.T @ G957 958 pvals = np.ones(m); betas = np.zeros(m); ses = np.zeros(m)959 for j in range(m):960 if G[:, j].std() < 1e-8:961 continue962 Xt = np.column_stack([X0t, Gt[:, j]])963 res = minimize_scalar(lambda ld: -_profile_loglik(ld, evals, yt, Xt, reml=True),964 bounds=(-10, 10), method="bounded",965 options={"xatol": 1e-3})966 delta = float(np.exp(res.x))967 p, b, se = _wald_test_rotated(yt, Xt, evals, delta)968 pvals[j] = p; betas[j] = b; ses[j] = se969 return pvals, betas, ses970 971# ─────────────────────── CONFIG BUILDER ──────────────────────────────────────────972def build_cfg(w, h, fs, fc, gc="#E8ECF4", bg="#FFFFFF", panel="#FAFBFF",973 sig_col="#E15759", sug_col="#4E79A7", marker_size=12, alpha=0.75,974 font_color=None, dot_col=None, line_col=None, fill_col=None,975 axis_col=None):976 """977 Build a plot configuration dictionary.978 'fc' is the required positional font-color argument.979 'font_color' is an optional alias (takes precedence if provided).980 """981 effective_fc = font_color if font_color is not None else fc982 return {983 "width": w,984 "height": h,985 "font_size": fs,986 "font_color": effective_fc,987 "grid_color": gc,988 "bg_color": bg,989 "panel_color": panel,990 "sig_color": sig_col,991 "sug_color": sug_col,992 "marker_size": marker_size,993 "alpha": alpha,994 "dot_color": dot_col or "#4E79A7",995 "line_color": line_col or "#E15759",996 "fill_color": fill_col or "#A78BFA",997 "axis_color": axis_col or "#2D3142",998 }999 1000# ─────────────────────── MANHATTAN PLOT ──────────────────────────────────────────1001def plot_manhattan(pvals, chroms, positions, title="Manhattan Plot", model_name="",1002 cfg=None, top_n=15, lod_threshold=None,1003 pve_vals=None, pve_threshold=None, show_logp_axis=True):1004 """1005 LOD is the selection criterion for this pipeline, so the LEFT axis1006 always shows the LOD scale. The RIGHT axis is an optional secondary1007 -log10(p) scale (show_logp_axis toggles it on/off) -- both axes are1008 exact monotonic transforms of the same underlying data1009 (neglog10p_to_lod / lod_to_neglog10p), not independent scales.1010 """1011 if cfg is None:1012 cfg = build_cfg(16, 5, 9, "#2D3142")1013 1014 n_snps = len(pvals)1015 bonf_thresh = 0.05 / n_snps1016 sug_thresh = 1.0 / n_snps1017 1018 sig_line = -np.log10(bonf_thresh)1019 sug_line = -np.log10(sug_thresh)1020 1021 fig, ax = plt.subplots(figsize=(cfg["width"], cfg["height"]))1022 fig.patch.set_facecolor(cfg["bg_color"])1023 ax.set_facecolor(cfg["panel_color"])1024 1025 log_p = -np.log10(np.clip(pvals, 1e-300, 1))1026 1027 chrom_list = sorted(set(chroms))1028 chrom_colors = {c: CHR_PALETTE[i % len(CHR_PALETTE)] for i, c in enumerate(chrom_list)}1029 x_offset = 01030 tick_pos, tick_lab, x_coords = [], [], np.zeros(len(pvals))1031 gap = max(positions) * 0.018 + 11032 1033 for ch in chrom_list:1034 mask = chroms == ch1035 pos = positions[mask]1036 if len(pos) == 0: continue1037 order = np.argsort(pos)1038 x_vals = pos[order] - pos.min() + x_offset1039 x_coords[np.where(mask)[0][order]] = x_vals1040 tick_pos.append(x_vals.mean())1041 tick_lab.append(f"Chr{ch}")1042 ax.scatter(x_coords[mask], log_p[mask],1043 c=chrom_colors[ch], s=cfg["marker_size"],1044 alpha=cfg["alpha"], linewidths=0)1045 x_offset += (pos.max() - pos.min()) + gap1046 1047 ax.axhline(sug_line, color=cfg["sug_color"], ls=":", lw=1.5, alpha=0.9,1048 label=f"Suggestive (1/m = {sug_thresh:.1e}, -log₁₀={sug_line:.1f})")1049 ax.axhline(sig_line, color=cfg["sig_color"], ls="--", lw=2.0,1050 label=f"Bonferroni (0.05/m = {bonf_thresh:.1e}, -log₁₀={sig_line:.1f})")1051 if lod_threshold is not None:1052 lod_line_negp = lod_to_neglog10p(lod_threshold)1053 ax.axhline(lod_line_negp, color=ACCENT_LIME, ls="-", lw=2.0, alpha=0.95,1054 label=f"Permutation / LOD threshold (LOD={lod_threshold:.2f})")1055 1056 # PRIMARY significance criterion: markers passing the LOD/permutation1057 # threshold (optionally AND a minimum PVE%), not a blind top-N slice.1058 # If no threshold could be computed, falls back to a small top-N so the1059 # plot doesn't render empty.1060 sig_idx = get_significant_idx(pvals, lod_threshold, fallback_top_n=top_n,1061 pve_vals=pve_vals, pve_threshold=pve_threshold)1062 sig_mask = np.zeros(len(pvals), dtype=bool)1063 sig_mask[sig_idx] = True1064 1065 if sig_mask.any():1066 crit_bits = []1067 if lod_threshold is not None:1068 crit_bits.append("LOD threshold")1069 if pve_vals is not None and pve_threshold:1070 crit_bits.append(f"PVE ≥ {pve_threshold:.1f}%")1071 label = (f"{sig_mask.sum()} MTA(s) ≥ {' & '.join(crit_bits)}" if crit_bits1072 else f"Top {sig_mask.sum()} MTAs (no threshold available)")1073 ax.scatter(x_coords[sig_mask], log_p[sig_mask],1074 c=ACCENT_GOLD, s=max(cfg["marker_size"]*4, 60),1075 marker="D", zorder=10, edgecolors="#333333", linewidth=0.7,1076 label=label)1077 1078 ax.set_xticks(tick_pos)1079 ax.set_xticklabels(tick_lab, rotation=45,1080 fontsize=cfg["font_size"]-1, color=cfg["font_color"])1081 ax.set_xlabel("Chromosome", color=cfg["font_color"], fontsize=cfg["font_size"])1082 ax.set_title(f"{title}{' — ' + model_name if model_name else ''}",1083 color=cfg["font_color"], fontsize=cfg["font_size"]+2, pad=8)1084 ax.legend(fontsize=cfg["font_size"]-1, framealpha=0.9,1085 facecolor=cfg["bg_color"], edgecolor="#CCCCCC",1086 labelcolor=cfg["font_color"], loc='upper right')1087 ax.set_xlim(-x_offset * 0.01, x_offset * 1.01)1088 for spine in ax.spines.values():1089 spine.set_edgecolor("#CCCCCC")1090 ax.grid(color=cfg["grid_color"], linewidth=0.4, linestyle="--", alpha=0.6)1091 1092 # ── LEFT axis: the data is plotted on the -log10(p) scale, but every1093 # tick is *labeled* in LOD units via a direct formatter -- this always1094 # renders (unlike a stacked secondary axis, which can silently collapse1095 # to zero width under tight_layout). LOD and -log10(p) are exact1096 # monotonic transforms of each other, so this is just a relabeling.1097 ax.yaxis.set_major_formatter(1098 plt.FuncFormatter(lambda val, pos: f"{neglog10p_to_lod(val):.1f}")1099 )1100 ax.set_ylabel("LOD score", color=cfg["font_color"], fontsize=cfg["font_size"])1101 ax.tick_params(axis='y', colors=cfg["font_color"], labelsize=cfg["font_size"]-1)1102 1103 if show_logp_axis:1104 ax_logp = ax.secondary_yaxis('right')1105 ax_logp.set_ylabel("-log₁₀(p)", color=cfg["font_color"], fontsize=cfg["font_size"])1106 ax_logp.tick_params(colors=cfg["font_color"], labelsize=cfg["font_size"]-1)1107 ax_logp.spines["right"].set_edgecolor("#CCCCCC")1108 1109 fig.tight_layout()1110 return fig1111 1112 1113# ─────────────────────── QQ PLOT ─────────────────────────────────────────────────1114def _nature_qq_ticks(max_val):1115 """Pick clean, evenly-spaced tick values (Nature-style QQ plots use round1116 steps like 0,3,6,9,12 rather than whatever matplotlib defaults to)."""1117 if not np.isfinite(max_val) or max_val <= 0:1118 return [0, 1]1119 for step in (1, 2, 3, 5, 10, 15, 20, 25, 50, 100):1120 n_ticks = max_val / step1121 if n_ticks <= 5:1122 top = int(np.ceil(max_val / step)) * step1123 return list(np.arange(0, top + step, step))1124 top = int(np.ceil(max_val / 100)) * 1001125 return list(np.arange(0, top + 100, 100))1126 1127 1128def plot_qq_single(pvals, title="QQ Plot", model_name="", color=None, cfg=None,1129 panel_letter=None):1130 """Clean single-series QQ plot in a minimal Nature-journal style: plain1131 white background, only left/bottom axis lines, a thin black y=x1132 reference line, no gridlines, and a small unobtrusive lambda label."""1133 if cfg is None:1134 cfg = build_cfg(7, 7, 9, "#000000")1135 if color is None:1136 color = cfg.get("dot_color", "#3E7BD6")1137 1138 fig, ax = plt.subplots(figsize=(cfg["width"], cfg["height"]))1139 fig.patch.set_facecolor("#FFFFFF")1140 ax.set_facecolor("#FFFFFF")1141 1142 n = len(pvals)1143 p_sorted = np.sort(np.clip(pvals, 1e-300, 1))1144 obs = -np.log10(p_sorted)1145 exp = -np.log10((np.arange(1, n + 1) - 0.5) / n)1146 1147 ax.scatter(exp, obs, c=color, s=max(cfg["marker_size"] * 0.6, 8),1148 alpha=0.9, linewidths=0, zorder=5)1149 1150 axis_max = float(np.nanmax([exp.max() if n else 1, obs.max() if n else 1]))1151 ticks = _nature_qq_ticks(axis_max)1152 lim = ticks[-1]1153 ax.plot([0, lim], [0, lim], color="#000000", lw=1.3, zorder=4, solid_capstyle="round")1154 1155 lam = np.median(stats.chi2.ppf(1 - np.clip(pvals, 1e-10, 1), df=1)) / stats.chi2.ppf(0.5, df=1)1156 ax.text(0.97, 0.04, f"\u03bb = {lam:.3f}", transform=ax.transAxes,1157 fontsize=cfg["font_size"] - 2, color="#888888", ha="right", va="bottom")1158 1159 ax.set_xlim(0, lim)1160 ax.set_ylim(0, lim)1161 ax.set_xticks(ticks)1162 ax.set_yticks(ticks)1163 ax.set_xlabel("Expected \u2212log$_{10}$($P$ value)", fontsize=cfg["font_size"] + 1,1164 color="#000000")1165 ax.set_ylabel("Observed \u2212log$_{10}$($P$ value)", fontsize=cfg["font_size"] + 1,1166 color="#000000")1167 if title or model_name:1168 ax.set_title(f"{title}{' — ' + model_name if model_name else ''}",1169 color="#555555", fontsize=cfg["font_size"], pad=10)1170 if panel_letter:1171 ax.text(-0.16, 1.04, panel_letter, transform=ax.transAxes,1172 fontsize=cfg["font_size"] + 8, fontweight="bold", color="#000000",1173 ha="left", va="top")1174 1175 ax.tick_params(colors="#000000", labelsize=cfg["font_size"], length=4, width=1.1)1176 for side in ("top", "right"):1177 ax.spines[side].set_visible(False)1178 for side in ("left", "bottom"):1179 ax.spines[side].set_color("#000000")1180 ax.spines[side].set_linewidth(1.2)1181 ax.grid(False)1182 fig.tight_layout()1183 return fig1184 1185 1186def plot_qq_all_models(all_results, cfg=None):1187 """QQ plots for all models with regression lines."""1188 if cfg is None:1189 cfg = build_cfg(7, 6, 9, "#2D3142")1190 n_models = len(all_results)1191 ncols = min(3, n_models)1192 nrows = (n_models + ncols - 1) // ncols1193 fig, axes = plt.subplots(nrows, ncols,1194 figsize=(cfg["width"] * ncols, cfg["height"] * nrows))1195 fig.patch.set_facecolor(cfg["bg_color"])1196 1197 if n_models == 1:1198 axes = np.array([axes])1199 axes = np.array(axes).flatten()1200