From 8d91f50c1c861e899c347c781de64ee04a850dab Mon Sep 17 00:00:00 2001 From: kakyungkim Date: Mon, 31 Aug 2026 10:33:02 +0900 Subject: [PATCH] =?UTF-8?q?P5=20ATAC=E2=86=92=CE=B1=20=EB=AA=A8=ED=98=95?= =?UTF-8?q?=20=EA=B3=84=EC=97=B4=20=EA=B2=AC=EA=B3=A0=EC=84=B1=20=EC=A0=90?= =?UTF-8?q?=EA=B2=80(P5c)=20=E2=80=94=20=EB=B9=84=EC=84=A0=ED=98=95?= =?UTF-8?q?=EC=9D=B4=20=EC=84=A0=ED=98=95=EC=9D=84=20=EB=84=98=EC=A7=80=20?= =?UTF-8?q?=EB=AA=BB=ED=95=A8?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit P5b(atac_alpha_expression_confound.md)의 "선형 모형이라 신호를 놓친 것 아닌가" 반론에 답한다. 같은 ATAC peak 특징 6종, 같은 Leave-One-Lineage-Out 규약에서 모형 계열만 바꿔 비교했다. - 선형(RidgeCV) held-out ρ=+0.304 / 발현 통제 partial ρ=+0.109 (n=472) → P5b 기록(+0.309 / +0.112 / n=472) 재현 확인 - RandomForest +0.178 / +0.013, GradientBoosting +0.206 / +0.056 → 비선형은 열등하고 발현 통제 시 격차가 더 벌어짐 상태는 탐색적이며 원고 반영 전 팀 검토가 필요하다. 계보 홀드아웃 한 규약만 확인했고 하이퍼파라미터는 조정하지 않았다. --- .../results/atac_alpha_model_class_check.md | 53 ++++++++++++++++ .../scripts/p5c_alpha_model_class_check.py | 62 +++++++++++++++++++ 2 files changed, 115 insertions(+) create mode 100644 pipeline/hspc-velocity-benchmark/results/atac_alpha_model_class_check.md create mode 100644 pipeline/hspc-velocity-benchmark/scripts/p5c_alpha_model_class_check.py diff --git a/pipeline/hspc-velocity-benchmark/results/atac_alpha_model_class_check.md b/pipeline/hspc-velocity-benchmark/results/atac_alpha_model_class_check.md new file mode 100644 index 0000000..eeb9669 --- /dev/null +++ b/pipeline/hspc-velocity-benchmark/results/atac_alpha_model_class_check.md @@ -0,0 +1,53 @@ +# P5c — baseline ATAC→α 의 약한 신호가 선형 모형 탓인가 + +> 2026-08-30, kkkim. `atac_alpha_expression_confound.md`(P5b)의 후속 견고성 점검. +> **상태: 탐색적.** 커밋·원고 반영 전 팀 검토 필요. + +## 왜 + +P5b에서 baseline ATAC→α 는 raw held-out ρ=+0.309, 발현 통제 partial ρ=+0.112 (n=472)로 +발현 confound에 크게 잠식됐다. 리뷰어가 물을 수 있는 반론이 하나 남는다. +**"선형 모형이라 비선형 관계를 놓친 것 아닌가."** 모형 계열을 바꿔 확인한다. + +## 방법 + +- 특징: `atac_baseline_features.csv`의 진짜 ATAC peak 6종 + (prom_acc, enh_acc, enh_sum, n_prom, n_enh, prom_enh_ratio) +- 표적: `lag_model.csv`의 `fit_alpha`, 계보 라벨도 같은 파일 +- 분할: Leave-One-Lineage-Out (계보 6개). P5b와 같은 계보 홀드아웃 규약 +- 발현 통제: `coupling_per_gene.csv`의 abundance에 대해 예측·실측을 각각 순위 회귀한 잔차끼리 상관 +- 모형: RidgeCV(선형), RandomForest(500, min_samples_leaf=5), GradientBoosting(기본) + +## 재현 확인 + +선형 모형이 P5b 수치를 재현한다. **held-out ρ=+0.304**(기록 +0.309), +**partial ρ=+0.109**(기록 +0.112), **n=472**(기록 n=472). 설정이 일치한다. + +## 결과 + +| 모형 | held-out ρ | p | 발현 통제 partial ρ | +|---|---|---|---| +| linear (RidgeCV) | **+0.304** | 1.47e-11 | **+0.109** | +| RandomForest | +0.178 | 1.02e-04 | +0.013 | +| GradientBoosting | +0.206 | 6.64e-06 | +0.056 | + +개별 특징의 fit_alpha 상관: enh_n +0.352, enh_sum +0.318, enh_acc +0.228, +prom_enh_ratio −0.246, prom_acc −0.128, n_prom +0.046 + +## 해석 + +**비선형 모형은 선형을 넘지 못하고 오히려 떨어진다.** 발현을 통제하면 격차가 더 벌어져 +RandomForest의 partial은 +0.013으로 사실상 0이다. n=472에 특징 6개인 조건에서 유연한 모형이 +과적합해 일반화가 나빠지는 전형적인 양상이다. + +따라서 "선형 모형이라 신호를 놓쳤다"는 반론은 닫힌다. ATAC→α 신호가 발현 통제 후 약한 것은 +모형 계열의 한계가 아니라 신호 자체의 성질이다. P5b의 결론을 약화시키지 않고 오히려 보강한다. + +## 한계 + +- 계보 홀드아웃 한 규약만 봤다. 유전자 무작위 분할은 확인하지 않았다. +- 하이퍼파라미터를 조정하지 않았다(기본값 + 최소 규제). 튜닝하면 비선형이 선형에 근접할 수는 + 있으나 넘어설 근거는 이 표본 크기에서 기대하기 어렵다. +- 발현 통제는 순위 선형 잔차 방식이다. P5b가 쓴 통제 방식과 완전히 동일한지 대조하지 않았다. + +재현: `scripts/p5c_alpha_model_class_check.py` diff --git a/pipeline/hspc-velocity-benchmark/scripts/p5c_alpha_model_class_check.py b/pipeline/hspc-velocity-benchmark/scripts/p5c_alpha_model_class_check.py new file mode 100644 index 0000000..8e62310 --- /dev/null +++ b/pipeline/hspc-velocity-benchmark/scripts/p5c_alpha_model_class_check.py @@ -0,0 +1,62 @@ +"""P5c — baseline ATAC→α 의 약한 신호가 선형 모형 탓인지 확인한다. + +P5b(atac_alpha_expression_confound.md)에서 ATAC→α 는 raw held-out ρ=+0.309, +발현 통제 partial ρ=+0.112 로 발현 confound에 잠식됐다. "선형 모형이라 놓친 것 아닌가"라는 +반론에 답하기 위해 같은 특징·같은 계보 홀드아웃에서 모형 계열만 바꿔 비교한다. + +실행: python scripts/p5c_alpha_model_class_check.py (결과 dir 기준 상대경로) +""" +import numpy as np +import pandas as pd +from scipy.stats import spearmanr +from sklearn.ensemble import RandomForestRegressor, GradientBoostingRegressor +from sklearn.linear_model import RidgeCV +from sklearn.model_selection import LeaveOneGroupOut +from sklearn.preprocessing import StandardScaler + +R = "results" + +# 진짜 ATAC peak 특징을 쓴다. lag_model.csv의 base_acc/chrom_rng/acc_mean은 +# ATAC peak이 아니라 다른 양이므로 이 분석에 쓰지 않는다(부호가 달라진다). +feat = pd.read_csv(f"{R}/atac_baseline_features.csv").set_index("gene") +lagm = pd.read_csv(f"{R}/lag_model.csv", index_col=0) +abund = pd.read_csv(f"{R}/coupling_per_gene.csv").set_index("gene")["abundance"] + +FEATS = list(feat.columns) +d = (feat.join(lagm[["fit_alpha", "lineage"]], how="inner") + .join(abund, how="left") + .dropna(subset=FEATS + ["fit_alpha", "lineage"])) +print(f"병합 {d.shape} | 특징 {FEATS}") +print("개별 상관(vs fit_alpha):", + {c: round(spearmanr(d[c], d.fit_alpha)[0], 3) for c in FEATS}) + +X, y, groups = d[FEATS].values, d.fit_alpha.values, d.lineage.values +MODELS = { + "linear(RidgeCV)": lambda: RidgeCV(alphas=np.logspace(-3, 3, 25)), + "RandomForest": lambda: RandomForestRegressor( + n_estimators=500, min_samples_leaf=5, random_state=0, n_jobs=-1), + "GradBoost": lambda: GradientBoostingRegressor(random_state=0), +} + +ok = d.abundance.notna().values +rank_ab = np.argsort(np.argsort(d.abundance.values[ok])) + + +def residual_vs_abundance(v): + """발현 순위에 대한 선형 잔차. 발현 confound를 뺀 뒤 상관을 본다.""" + rv = np.argsort(np.argsort(v[ok])) + return rv - np.polyval(np.polyfit(rank_ab, rv, 1), rank_ab) + + +print(f"\n{'model':18s} {'held-out rho':>13s} {'p':>10s} {'partial(발현통제)':>18s}") +for name, make in MODELS.items(): + pred = np.full(len(y), np.nan) + for tr, te in LeaveOneGroupOut().split(X, y, groups): # 계보 하나를 빼고 학습 + sc = StandardScaler().fit(X[tr]) + pred[te] = make().fit(sc.transform(X[tr]), y[tr]).predict(sc.transform(X[te])) + rho, p = spearmanr(pred, y) + prho, _ = spearmanr(residual_vs_abundance(pred), residual_vs_abundance(y)) + print(f"{name:18s} {rho:+13.3f} {p:10.2e} {prho:+18.3f}") + +print("\nP5b 기록: held-out rho=+0.309, 발현통제 partial=+0.112, n=472") +print("선형이 이 값을 재현하면 설정이 일치한다는 뜻이다.")