Lesson 18 of 20

차등발현 — 반복수가 p값보다 중요하다

왜 음이항이고, 왜 p값을 그대로 쓰면 안 되나

이제 「대조군과 처리군에서 다르게 발현된 유전자」를 찾습니다. 유전자마다 t검정을 하면 될 것 같지만, 두 가지가 막습니다.

① 반복이 너무 적다

시료가 3개씩
분산을 3개로 추정할 수 없다
t검정이 무의미해진다

② 검정을 2만 번 한다

p < 0.05 를 2만 번 하면
차이가 전혀 없어도
1,000개가 유의하게 나온다

High school반복수가 p값보다 중요하다

차등발현 분석에서 가장 중요한 것은 통계 방법이 아니라 생물학적 반복수입니다.

반복수	할 수 있는 것
1	분산을 추정 못 함 — 통계 불가
2	최소한 · 매우 불안정
3	흔한 최소
5~6 이상	안정적

2강에서 "깊이를 어디에 쓸 것인가"를 물었죠 — 답은 깊이보다 반복수입니다.

이 셀 하나가 padj 를 써야 하는 이유의 전부입니다.

Python Runs in your browser

# 다중검정 — 아무 차이가 없는데 몇 개가 유의하게 나오나
import numpy as np
rng = np.random.default_rng(21)

n_gene = 20000
# 두 조건 사이에 진짜 차이가 전혀 없다고 하자 -> p값은 균등분포
p = rng.uniform(0, 1, n_gene)

print("유전자 %s개, 진짜 차이는 하나도 없다" % format(n_gene, ","))
print()
print("p < 0.05 인 유전자: %s개" % format((p < 0.05).sum(), ","))
print("p < 0.01 인 유전자: %s개" % format((p < 0.01).sum(), ","))
print()

# Benjamini-Hochberg 보정
order = np.argsort(p)
ranked = p[order]
m = len(p)
bh = ranked * m / np.arange(1, m + 1)
bh = np.minimum.accumulate(bh[::-1])[::-1]
print("BH 보정(FDR) 뒤 0.05 미만: %d개" % (bh < 0.05).sum())
print()
print("보정 없이 발표하면 1,000개짜리 '발견' 목록이 나옵니다.")
print("전부 우연입니다.")

High school왜 깊이보다 반복수인가

같은 예산을	얻는 것
3시료 × 6,000만 read	깊이는 충분 · 분산 추정이 불안정
6시료 × 3,000만 read	깊이도 충분 · 분산 추정이 훨씬 좋음

발현 정량은 2,000만~3,000만 read면 포화됩니다 — 더 읽어도 얻는 것이 적어요.

흩어짐 하나가 들어가면 값이 이만큼 흔들립니다.

Python Runs in your browser

# 푸아송으로는 실제 자료를 못 담는다
import numpy as np
rng = np.random.default_rng(17)

mu = 500.0
n = 6

pois = rng.poisson(mu, n)

# 음이항 — 생물학적 개체차를 담는다 (흩어짐 0.1)
disp = 0.1
p_nb = 1 / (1 + mu * disp)
nb = rng.negative_binomial(1 / disp, p_nb, n)

print("평균 %d 인 유전자를 시료 %d개에서 관측한다면" % (mu, n))
print()
print("푸아송 가정:", list(pois))
print("   분산 %8.0f   (평균과 비슷하다)" % pois.var(ddof=1))
print()
print("실제 (음이항):", list(nb))
print("   분산 %8.0f   (훨씬 크다)" % nb.var(ddof=1))
print()
print("푸아송을 쓰면 분산을 낮게 잡아 p값이 지나치게 작아집니다.")
print("차이가 없는 유전자를 유의하다고 발표하게 됩니다.")

High school왜 음이항 분포인가

카운트 자료는 포아송 분포를 따를 것 같습니다. 그런데 실제로는 더 흩어져 있어요.

분포	분산	맞나
포아송	= 평균	기술적 반복에만
음이항	> 평균	생물학적 반복에

개체마다 진짜 발현이 다릅니다 — 그 차이가 추가 분산을 만듭니다.

반복수 — 여기가 진짜 문제다

반복수	잡을 수 있는 것	쓸 만한가
2	분산을 추정할 수 없다	❌ 하지 말 것
3	4배 이상 차이 나는 유전자	최소한
5~6	2배 차이	✅ 권장
10 이상	1.5배 차이	✅ 미묘한 효과를 볼 때

깊이를 두 배로 늘리는 것보다 시료를 하나 더 넣는 것 이 낫다

이것이 이 강의 전체의 결론 중 하나입니다. 2강에서 「깊이를 어디에 쓸 것인가」를 물었는데, 답이 여기 있어요.

High school분산을 어떻게 추정하나

3개 표본으로 유전자마다 분산을 추정하면 아주 불안정합니다. 그래서 빌려 옵니다.

방법	내용
정보 공유	비슷한 발현 수준의 유전자들에서 분산 경향을 학습
수축 (shrinkage)	개별 추정을 그 경향 쪽으로 당김

유전자 2만 개가 서로 도와 분산을 추정합니다 — 그래서 적은 반복수로도 됩니다.

예산을 어디에 쓸지 정하는 셈.

Python Runs in your browser

# 반복수와 깊이 — 어느 쪽에 돈을 쓸까
import numpy as np
rng = np.random.default_rng(31)

def power(n_rep, depth_m, lfc, trials=400):
    # 깊이는 기술적 잡음만 줄인다. 개체차는 그대로 남는다.
    bio_sd  = 0.35                        # 생물학적 개체차 (log2)
    tech_sd = 1.0 / np.sqrt(depth_m)      # 깊이가 늘면 준다
    sd = np.sqrt(bio_sd ** 2 + tech_sd ** 2)
    hits = 0
    for _ in range(trials):
        a = rng.normal(0.0, sd, n_rep)
        b = rng.normal(lfc, sd, n_rep)
        se = np.sqrt(a.var(ddof=1) / n_rep + b.var(ddof=1) / n_rep)
        if se > 0 and abs(b.mean() - a.mean()) / se > 2.5:
            hits += 1
    return hits / trials

print("2배 차이(log2FC = 1)를 잡을 확률")
print()
print("%10s %12s %12s" % ("", "깊이 20M", "깊이 40M"))
for n in [3, 4, 6, 8]:
    print("%8d개 %11.0f%% %11.0f%%"
          % (n, power(n, 20, 1.0) * 100, power(n, 40, 1.0) * 100))
print()
print("깊이를 두 배로 해도 거의 안 오릅니다.")
print("반복을 3에서 6으로 늘리는 쪽이 훨씬 큽니다.")
print("생물학적 개체차는 깊이로 줄일 수 없기 때문입니다.")

결과 표에서 보는 것

열	무엇	보는 법
`baseMean`	모든 시료의 평균 카운트	10 아래면 믿지 않는다
`log2FoldChange`	2를 밑으로 한 배수	1 이면 두 배
`lfcSE`	그 배수의 오차	크면 반복이 모자란다
`pvalue`	보정 전	쓰지 않는다
`padj`	BH 보정 뒤 (FDR)	이것으로 자른다 — 보통 0.05

자르는 관습은 padj < 0.05 그리고 |log2FC| > 1 입니다. 둘 다 걸어야 해요 — 통계적으로 유의해도 1.05배 차이는 뜻이 없습니다.

High school왜 p값을 그대로 쓰면 안 되나

유전자 2만 개를 각각 검정합니다. 차이가 없어도 1,000개가 우연히 p < 0.05가 돼요.

보정	결과
없음	가짜가 1,000개
본페로니	거의 아무것도 안 남음
BH (FDR)	고른 것 중 가짜 비율을 통제

q < 0.05 = "이 목록에 가짜가 5% 정도 섞여 있다"

그리고 fold change 문턱도 같이 겁니다. 통계적으로 유의해도 1.05배 차이는 생물학적 의미가 적어요.

🤔 Check yourself

유전자 2만 개를 p < 0.05 로 검정하면, 차이가 전혀 없어도 몇 개가 유의하게 나오나?

Answer 20,000 × 0.05 = 1,000개. 그래서 BH 보정을 거친 padj 로 자릅니다.

🤔 Check yourself

예산이 조금 더 생겼다. 무엇에 쓰는 것이 나은가?