Lesson 18 of 20

찾을 때까지 뒤지면 나온다 — 문턱과 유전자 목록

GAPDH 와 ACTB 가 '생존과 유의하게 관련' 되는 방법

앞의 강에서 자료를 읽는 법을 배웠습니다. 이제 가장 흔한 잘못 하나를 직접 저질러 보겠습니다. 저질러 봐야 안 저지를 수 있어요.

물음은 이것입니다 — TMB 가 높은 폐선암 환자는 오래 사는가. 면역치료 이야기에서 자주 나오는 주장이지요.

먼저 정직하게 해 봅니다. 중앙값으로 반을 갈라 로그순위 검정을 합니다.

High school찾을 때까지 뒤지면 나온다

이 강의의 경고입니다. 유전자 2만 개와 여러 문턱을 다 시도하면 무엇이든 유의해집니다.

시도할 수 있는 것	가짓수
유전자	2만
문턱 (중앙값·사분위·최적)	여럿
종점 (OS·DSS·DFI·PFI)	4
암종	33

조합이 수백만입니다 — 우연히 유의한 것이 반드시 나옵니다.

p = 0.53. 차이가 없습니다. 여기서 끝내는 것이 정직한 분석이에요.

Python Runs in your browser

import pandas as pd
import tcga

c = tcga.clinical()
s = tcga.luad_sample()
m = s.merge(c, on="sample", how="inner")
m = m[m["OS"].notna() & m["OS_time"].notna() & m["tmb"].notna()]

문턱 = m["tmb"].median()
높음 = m[m["tmb"] > 문턱]
낮음 = m[m["tmb"] <= 문턱]

chi2, p = tcga.logrank(높음["OS_time"], 높음["OS"], 낮음["OS_time"], 낮음["OS"])
print("TMB 중앙값 %.2f 로 분할" % 문턱)
print("   높음 %d명 · 중앙생존 %s일" % (len(높음), tcga.median_survival(높음["OS_time"], 높음["OS"])))
print("   낮음 %d명 · 중앙생존 %s일" % (len(낮음), tcga.median_survival(낮음["OS_time"], 낮음["OS"])))
print()
print("로그순위  chi2 = %.3f   p = %.4f" % (chi2, p))

그런데 논문을 써야 합니다. 그래서 이렇게 생각해 봅니다 — "중앙값이 최선의 문턱이 아닐 수도 있잖아? 여러 문턱을 다 해 보자."

17번을 훑었는데 가장 작은 p 도 0.157 입니다. 이 자료에서는 문턱을 아무리 옮겨도 안 나와요.

Python Runs in your browser

import pandas as pd
import tcga

c = tcga.clinical()
s = tcga.luad_sample()
m = s.merge(c, on="sample", how="inner")
m = m[m["OS"].notna() & m["OS_time"].notna() & m["tmb"].notna()]

print("%6s %8s %10s %9s" % ("백분위", "문턱", "표본(고/저)", "p"))
최소 = (1.0, None)
for q in range(10, 91, 5):
    문턱 = m["tmb"].quantile(q / 100)
    a, b = m[m["tmb"] > 문턱], m[m["tmb"] <= 문턱]
    if len(a) < 20 or len(b) < 20:
        continue
    p = tcga.logrank(a["OS_time"], a["OS"], b["OS_time"], b["OS"])[1]
    표시 = " ←" if p < 최소[0] else ""
    if p < 최소[0]:
        최소 = (p, q)
    print("%5d%% %8.2f %6d/%-4d %9.4f%s" % (q, 문턱, len(a), len(b), p, 표시))
print()
print("가장 작은 p = %.4f  (백분위 %d%%)" % 최소)

High schoolGAPDH도 유의해진다

GAPDH와 ACTB는 기준 유전자입니다. 생존과 관련 있을 리 없죠. 그런데 나옵니다.

어떻게	내용
여러 암종을 다 시도	어딘가에서 유의
문턱을 최적화	가장 잘 갈리는 지점을 찾음
종점을 바꿔 봄	—

"유의한 결과"를 만드는 것은 놀랍도록 쉽습니다.

문턱으로 안 되면 다음 유혹은 이것입니다 — "유전자를 전부 돌려 보자. 하나는 걸리겠지."

High school최적 문턱의 함정

가장 위험한 방법입니다.

단계	내용
1	모든 문턱에서 생존 차이를 계산
2	p값이 가장 작은 문턱을 고름
3	그 p값을 보고

수백 번 검정하고 가장 좋은 것만 보고하는 셈입니다 — 보정을 안 하면 사실상 무의미해요.

네 개가 걸렸습니다. 그런데 그 안에 ACTB 와 GAPDH 가 있어요.

Python Runs in your browser

import numpy as np
import pandas as pd
import tcga

c = tcga.clinical()
s = tcga.luad_sample()
e = tcga.expr()
m = s.merge(c, on="sample", how="inner")
m = m[m["OS"].notna() & m["OS_time"].notna()].set_index("sample")

공통 = [x for x in e.columns if x in m.index]
mm = m.loc[공통]

결과 = []
for 유전자 in e.index:
    v = np.log2(e.loc[유전자, 공통] + 1)
    있음 = v.notna().values                 # 발현이 빈 표본은 그 유전자에서만 뺀다
    v2, d = v[있음], mm[있음]
    문턱 = v2.median()
    hi, lo = d[(v2 > 문턱).values], d[(v2 <= 문턱).values]
    p = tcga.logrank(hi["OS_time"], hi["OS"], lo["OS_time"], lo["OS"])[1]
    결과.append((유전자, p))

결과.sort(key=lambda x: x[1])
print("유전자 %d개 · 표본 %d명 — 발현 중앙값으로 나눈 생존 검정" % (len(결과), len(공통)))
print()
for 유전자, p in 결과[:8]:
    print("   %-8s p = %.4f  %s" % (유전자, p, "*" if p < 0.05 else ""))
print("   ...")
print()
print("p < 0.05 인 유전자: %d개" % sum(1 for _, p in 결과 if p < 0.05))

보정하면 MKI67 과 ALK 만 겨우 남습니다. GAPDH 와 ACTB 는 떨어져요 — 살림유전자가 떨어지는 것이 정상입니다.

Python Runs in your browser

import numpy as np
import pandas as pd
import tcga

c = tcga.clinical(); s = tcga.luad_sample(); e = tcga.expr()
m = s.merge(c, on="sample", how="inner")
m = m[m["OS"].notna() & m["OS_time"].notna()].set_index("sample")
공통 = [x for x in e.columns if x in m.index]
mm = m.loc[공통]

결과 = []
for 유전자 in e.index:
    v = np.log2(e.loc[유전자, 공통] + 1)
    있음 = v.notna().values
    v2, d = v[있음], mm[있음]
    hi, lo = d[(v2 > v2.median()).values], d[(v2 <= v2.median()).values]
    결과.append((유전자, tcga.logrank(hi["OS_time"], hi["OS"], lo["OS_time"], lo["OS"])[1]))
결과.sort(key=lambda x: x[1])

n = len(결과)
문턱 = 0.05 / n
print("본페로니 문턱 = 0.05 / %d = %.5f" % (n, 문턱))
print()
print("%-8s %9s %9s %s" % ("유전자", "p", "보정p", "판정"))
for i, (유전자, p) in enumerate(결과[:6], 1):
    print("%-8s %9.4f %9.4f %s"
          % (유전자, p, min(1.0, p * n), "남음" if p < 문턱 else "탈락"))

High school정직하게 하는 법

방법	내용
문턱을 미리 정함	중앙값 등 — 결과를 보기 전에
다중검정 보정	FDR
독립 코호트에서 확인	가장 강력
시도한 것을 다 보고	몇 개를 봤는지
효과크기도 함께	p값만이 아니라

세 번째 줄이 결정적입니다 — TCGA에서 찾고 다른 자료에서 확인하세요.

🤔 Check yourself

문턱 17개를 훑어 가장 작은 p = 0.03 을 얻었습니다. 이 값은?

🤔 Check yourself

19개 유전자 검정에서 GAPDH 가 p = 0.022 로 나왔습니다. 가장 그럴듯한 해석은?