11강 / 20강

1만 7천 번 물으면 855번은 그냥 걸린다

Benjamini Y, Hochberg Y. J R Stat Soc B 57:289-300 (1995)

p < 0.05 는 '차이가 없는데도 이만큼 갈릴 확률이 5% 미만' 이라는 뜻입니다. 바꿔 말하면, 차이가 없는 유전자 100개를 검사하면 그중 5개는 그냥 걸립니다.

우리는 17,097개를 검사했습니다.

17,097 × 0.05 ≈ 855개

차이가 하나도 없더라도 855개쯤은 p < 0.05 로 나옵니다. 그냥 나옵니다.

차이를 0 으로 놓고 만든 자료인데도 800개 넘게 걸립니다. 전부 가짜입니다.

Python 브라우저에서 실행

# 정말 그런지 만들어서 확인합니다. 차이가 '없는' 자료로요.
import numpy as np
from scipy import stats

rng = np.random.default_rng(1995)
n_gene, n_msi, n_mss = 17097, 43, 565

# 한 개씩 1만 7천 번 돌리면 느립니다. 1,500개씩 묶어 한꺼번에 계산합니다.
hits = 0
step = 1500
for start in range(0, n_gene, step):
    n = min(step, n_gene - start)
    A = rng.normal(0, 1, (n, n_msi))     # 두 무리를 같은 분포에서 뽑습니다
    B = rng.normal(0, 1, (n, n_mss))     # 진짜 차이는 0 입니다
    p = stats.ttest_ind(A, B, axis=1, equal_var=False)[1]
    hits += int((p < 0.05).sum())

print("진짜 차이가 하나도 없는 유전자", n_gene, "개를 검사했더니")
print("p < 0.05 로 걸린 것:", hits, "개")
print("기대값:", round(0.05 * n_gene))

고등학생두 가지 안전 기준

	본페로니 (FWER)	BH (FDR)
통제하는 것	가짜가 하나라도 나올 확률	고른 것 중 가짜의 비율
기준	p < 0.05/17097	순위에 따라 문턱이 달라짐
결과	거의 아무것도 안 남음	수십~수백 개
쓰는 곳	임상 최종 판정	탐색적 발굴

우리는 표적 후보를 찾는 중입니다. 후보 중 5%가 가짜여도 나중에 실험으로 거르면 돼요.

거짓발견율 — 다르게 묻기

옛 방법(본페로니)은 기준을 통째로 조입니다 — 0.05 대신 0.05/17097 을 쓰는 식으로요. 확실하지만 너무 가혹해서 진짜도 다 놓칩니다.

벤저미니-호크버그(BH) 는 물음을 바꿉니다.

본페로니가 묻는 것

"가짜가 하나라도 섞일 확률을 5% 아래로"
→ 너무 엄격합니다

BH 가 묻는 것

"내가 고른 것들 가운데 가짜 비율을 10% 아래로"
→ 쓸 만합니다

이 비율을 거짓발견율(FDR, false discovery rate) 이라 하고, 유전자마다 계산한 값을 q 값이라고 부릅니다.

세 번째 줄을 보세요. 보정값은 0.091 인데 q 는 0.0588 입니다 — 뒤쪽(5번째)의 보정값이 더 작아서 끌어내린 것입니다.

Python 브라우저에서 실행

# BH 를 손으로 짜 봅니다. 정렬하고, 곱하고, 뒤에서부터 최소를 취합니다.
import numpy as np

p = np.array([0.001, 0.008, 0.039, 0.041, 0.042, 0.06, 0.5])
n = len(p)

order = np.argsort(p)
ranked = p[order]
q = ranked * n / np.arange(1, n + 1)        # 순위로 나눠 보정
q = np.minimum.accumulate(q[::-1])[::-1]    # 뒤에서부터 최소 — 순서가 뒤집히지 않게

print("p 값      순위   보정   q 값")
for i in range(n):
    print(str(ranked[i]).ljust(9), str(i + 1).ljust(6),
          str(round(ranked[i] * n / (i + 1), 4)).ljust(6), round(q[i], 4))

고등학생BH 절차 — 계단을 내려온다

단계	하는 일
1	p값을 작은 것부터 줄 세운다 (순위 i, 총 m개)
2	각 유전자에 문턱 (i/m)×0.05를 준다
3	p < 문턱인 가장 큰 순위 k를 찾는다
4	1위부터 k위까지 전부 유의로 선언

1등에게는 엄한 문턱, 아래쪽에는 느슨한 문턱 — 계단을 만듭니다.

Benjamini와 Hochberg(1995)가 제안한 방법입니다. 지금은 유전체 자료 분석의 표준이에요.

2,649개가 579개로 줄어듭니다. WRN 은 q = 1.03e-03 으로 넉넉히 살아남습니다.

Python 브라우저에서 실행

# 진짜 자료에 걸어 봅니다.
import wrn

g = wrn.genes()
g["myq"] = wrn.bh(g["p"].values)

print("p < 0.05  :", int((g["p"] < 0.05).sum()), "개")
print("q < 0.10  :", int((g["myq"] < 0.10).sum()), "개")
print("q < 0.05  :", int((g["myq"] < 0.05).sum()), "개")
print()
w = g[g["gene"] == "WRN"]
print("WRN  p =", format(float(w["p"].iloc[0]), ".3e"),
      "  q =", format(float(w["myq"].iloc[0]), ".3e"))

고등학생q값이 뜻하는 것

말	맞나
"이 유전자가 진짜일 확률 95%"	아니다 — q는 목록 전체의 성질
"이 목록에 가짜가 5% 정도"	맞다
"가짜가 정확히 몇 개다"	아니다 — 기댓값의 상한

q < 0.05로 100개를 고르면 그중 5개쯤이 가짜입니다 — 어느 것인지는 모릅니다.

그리고 FDR을 통과해도 12강의 체가 남아 있습니다. 통계적으로 유의한 것과 표적으로 쓸 만한 것은 다르니까요.

🤔 확인 문제

유전자 17,097개를 p < 0.05 로 검사할 때, 진짜 차이가 하나도 없어도 몇 개쯤 걸립니까?

🤔 확인 문제

q 값(FDR)이 본페로니 보정보다 실제 연구에서 널리 쓰이는 이유는?