Lesson 19 of 20

여섯 시료를 처음부터 끝까지

【실습】 시료 준비에서 화산그림까지 한 줄로

실험 하나를 처음부터 끝까지 갑니다.

항목	정한 것	몇 강
시료	배양세포 · 대조군 3 + 처리군 3	8강
추출	실리카 컬럼 + on-column DNase	8강
품질	RIN 9.1~9.6 · Qubit 180~260 ng/µL	9강
라이브러리	polyA 선택 · 가닥특이적 · 500 ng 투입	10강
인덱스	고유형(UDI) 6쌍	11강
장비	NextSeq 1000 · P2 300 사이클 · 2×150	2·3강
깊이	시료당 읽기쌍 30백만 목표	2강
배치	여섯 시료를 한 런에 — 배치효과 없음	16강

High school여섯 시료를 처음부터 끝까지

단계	하는 일	강
1	시료 준비 · RNA 추출	8강
2	품질검사 (RIN 확인)	9강
3	라이브러리 제작	10강
4	인덱스 붙이고 풀링	11강
5	시퀀싱	3~6강
6	BCL → FASTQ	12강
7	품질검사와 트리밍	14강
8	정렬 · 카운트	15강
9	차등발현 · 화산그림	18강

먼저 런이 살았는지부터 봅니다.

Python Runs in your browser

# 1단계 — 런 요약을 받는다 (12강)
run = {
    "클러스터 밀도 K/mm2": 208,
    "%PF":                 89,
    "%Occupied":           93,
    "%>=Q30 (R1)":         92,
    "%>=Q30 (R2)":         88,
    "읽기쌍 (백만)":         196,
}
for k, v in run.items():
    print("%-22s %8s" % (k, v))
print()
print("시료 6개가 나눠 가지면 시료당 평균 %.1f백만 쌍" % (196 / 6))
print("목표 30백만에 조금 못 미치지만 쓸 만합니다.")

High school각 단계에서 기록할 것

단계	기록
시료	채취 시각 · 보존 방법 · 배치
RNA	농도 · A260/280 · A260/230 · RIN
라이브러리	키트 · 로트 · 인덱스 · 농도
시퀀싱	장비 · 실행 ID · 레인 · read 길이
분석	도구와 버전 · 참조 유전체 · 주석

이 기록이 20강의 방법 문단이 됩니다.

풀링이 잘 됐는지 여기서 확인합니다.

Python Runs in your browser

# 2단계 — 시료마다 몫이 얼마나 고른가 (11강)
import numpy as np

samples = ["Ctrl_1", "Ctrl_2", "Ctrl_3", "Trt_1", "Trt_2", "Trt_3"]
pairs_m = np.array([34.1, 29.8, 26.4, 38.2, 31.5, 36.0])

print("%-8s %12s %10s" % ("시료", "읽기쌍(백만)", "평균 대비"))
for s, p in zip(samples, pairs_m):
    print("%-8s %12.1f %9.0f%%" % (s, p, p / pairs_m.mean() * 100))
print()
print("가장 적은 시료 %.1f, 가장 많은 시료 %.1f — %.1f배 차이"
      % (pairs_m.min(), pairs_m.max(), pairs_m.max() / pairs_m.min()))
print()
print("1.4배 안쪽이면 정상 범위입니다. 정규화가 이 차이를 흡수합니다 (17강).")
print("3배를 넘으면 풀링을 다시 봐야 합니다.")

여기서 이상한 시료를 하나 찾아냅니다.

Python Runs in your browser

# 3단계 — 정렬 지표 (15·16강)
rows = [
    ("Ctrl_1", 94.2, 88.1, 1.8, 2.1),
    ("Ctrl_2", 93.8, 87.4, 2.0, 2.4),
    ("Ctrl_3", 94.5, 88.9, 1.6, 1.9),
    ("Trt_1",  93.1, 86.2, 2.2, 2.8),
    ("Trt_2",  94.0, 87.8, 1.9, 2.2),
    ("Trt_3",  88.4, 71.3, 9.6, 3.1),
]
print("%-8s %10s %10s %8s %8s" % ("시료", "정렬률", "엑손", "인트론", "rRNA"))
for s, aln, ex, intr, rrna in rows:
    flag = "  <-- ?" if ex < 80 else ""
    print("%-8s %9.1f%% %9.1f%% %7.1f%% %7.1f%%%s" % (s, aln, ex, intr, rrna, flag))
print()
print("Trt_3 의 엑손 비율이 혼자 낮고 인트론이 높습니다.")
print("gDNA 오염이거나 전구체 mRNA 가 많은 것입니다 (8·16강).")
print("정규화로는 못 고칩니다. 표시해 두고 PCA 에서 다시 봅니다.")

PC1 이 조건을 가르면 좋은 신호입니다. 배치를 가르면 문제예요.

Python Runs in your browser

# 4단계 — PCA 를 먼저 그린다 (16강)
import numpy as np
import matplotlib.pyplot as plt
rng = np.random.default_rng(19)

n_gene = 4000
samples = ["Ctrl_1", "Ctrl_2", "Ctrl_3", "Trt_1", "Trt_2", "Trt_3"]

base = rng.normal(6, 2.2, n_gene)
X = np.vstack([base + rng.normal(0, 0.35, n_gene) for _ in range(6)])
# 처리군에서 400개 유전자가 올라간다
up = rng.choice(n_gene, 400, replace=False)
X[3:, up] += 1.6
# Trt_3 은 gDNA 오염으로 전체가 조금 흔들린다
X[5] += rng.normal(0, 0.55, n_gene)

Xc = X - X.mean(axis=0)
U, S, Vt = np.linalg.svd(Xc, full_matrices=False)
pc = U[:, :2] * S[:2]
var = (S ** 2 / (S ** 2).sum())[:2] * 100

plt.figure(figsize=(5.6, 4.4))
for i, s in enumerate(samples):
    c = "#3b82f6" if s.startswith("Ctrl") else "#ef4444"
    plt.scatter(pc[i, 0], pc[i, 1], s=110, color=c)
    plt.annotate(s, (pc[i, 0], pc[i, 1]), fontsize=8,
                 xytext=(6, 4), textcoords="offset points")
plt.xlabel("PC1  (%.0f%%)" % var[0])
plt.ylabel("PC2  (%.0f%%)" % var[1])
plt.title("PCA — 조건별로 갈리는가")
plt.axhline(0, color="#cbd5e1", lw=.8); plt.axvline(0, color="#cbd5e1", lw=.8)
plt.grid(alpha=.25); plt.tight_layout(); plt.show()

High school중간중간 확인할 것

확인 시점	보는 것	이상하면
RNA 추출 후	RIN	다시 추출
라이브러리 후	크기 분포 · 농도	다시 제작
시퀀싱 중	클러스터 밀도 · Q30	중단 고려
FASTQ 후	read 수 · 품질	—
정렬 후	정렬률 · rRNA 비율	16강의 함정 확인
카운트 후	PCA로 배치 확인	보정 계획

마지막 줄이 중요합니다 — PCA에서 군이 아니라 배치로 갈리면 신호입니다.

카운트 표에서 유전자 목록이 나오는 자리.

Python Runs in your browser

# 5단계 — 차등발현 (17·18강)
import numpy as np
rng = np.random.default_rng(23)

n_gene = 12000
n = 3
disp = 0.08

mu = np.exp(rng.normal(4.2, 1.6, n_gene))
lfc = np.zeros(n_gene)
true_de = rng.choice(n_gene, 500, replace=False)
lfc[true_de] = rng.normal(0, 1.3, 500)

def nb(mean):
    p = 1 / (1 + mean * disp)
    return rng.negative_binomial(1 / disp, p)

ctrl = np.array([[nb(m) for m in mu] for _ in range(n)])
trt  = np.array([[nb(m * 2 ** l) for m, l in zip(mu, lfc)] for _ in range(n)])

# 크기인자 (중앙값비 — 17강)
allc = np.vstack([ctrl, trt]).astype(float)
ref = np.exp(np.mean(np.log(allc + 0.5), axis=0))
sf = np.median((allc + 0.5) / ref, axis=1)
norm = allc / sf[:, None]

base_mean = norm.mean(axis=0)
a, b = norm[:n], norm[n:]
obs_lfc = np.log2((b.mean(axis=0) + 1) / (a.mean(axis=0) + 1))
se = np.sqrt(a.var(axis=0, ddof=1) / n + b.var(axis=0, ddof=1) / n) / (base_mean + 1) / np.log(2)
stat = obs_lfc / np.maximum(se, 1e-3)

from scipy import stats as st
pval = 2 * (1 - st.norm.cdf(np.abs(stat)))

# BH 보정
order = np.argsort(pval)
m = len(pval)
bh = pval[order] * m / np.arange(1, m + 1)
bh = np.minimum.accumulate(bh[::-1])[::-1]
padj = np.empty(m); padj[order] = np.minimum(bh, 1)

keep = base_mean >= 10
sig = keep & (padj < 0.05) & (np.abs(obs_lfc) > 1)
print("검정한 유전자:            %s개" % format(int(keep.sum()), ","))
print("padj < 0.05:             %s개" % format(int((keep & (padj < 0.05)).sum()), ","))
print("padj < 0.05 이고 |LFC|>1: %s개" % format(int(sig.sum()), ","))
print("  올라간 것: %d개   내려간 것: %d개"
      % ((sig & (obs_lfc > 0)).sum(), (sig & (obs_lfc < 0)).sum()))
print()
print("문턱 두 개를 함께 걸어야 목록이 쓸 만한 크기가 됩니다.")

가로는 얼마나 변했나, 세로는 얼마나 믿을 만한가 — 둘 다 필요합니다.

Python Runs in your browser

# 6단계 — 화산그림
import numpy as np
import matplotlib.pyplot as plt
rng = np.random.default_rng(29)

n_gene = 12000
lfc = rng.normal(0, 0.35, n_gene)
de = rng.choice(n_gene, 500, replace=False)
lfc[de] = rng.normal(0, 1.6, 500)
base = rng.lognormal(4, 1.2, n_gene)

se = 0.9 / np.sqrt(base)
z = lfc / np.maximum(se, 0.08)
p = 2 * (1 - 0.5 * (1 + np.tanh(np.abs(z) * 0.79)))     # 정규 꼬리 근사
p = np.clip(p, 1e-30, 1)
nlp = -np.log10(p)

sig = (p < 0.001) & (np.abs(lfc) > 1)
plt.figure(figsize=(6.2, 4.6))
plt.scatter(lfc[~sig], nlp[~sig], s=4, alpha=.28, color="#94a3b8")
plt.scatter(lfc[sig & (lfc > 0)], nlp[sig & (lfc > 0)], s=7, alpha=.7, color="#ef4444")
plt.scatter(lfc[sig & (lfc < 0)], nlp[sig & (lfc < 0)], s=7, alpha=.7, color="#3b82f6")
plt.axvline(1, ls="--", lw=.9, color="#64748b")
plt.axvline(-1, ls="--", lw=.9, color="#64748b")
plt.axhline(3, ls="--", lw=.9, color="#64748b")
plt.xlabel("log2 fold change")
plt.ylabel("-log10 p")
plt.title("volcano — 오른쪽 위와 왼쪽 위가 목록이 된다")
plt.tight_layout(); plt.show()

High school화산그림까지

축	나타내는 것
가로	log2 fold change — 얼마나 달라졌나
세로	−log10(p 또는 q) — 얼마나 확실한가

왼쪽 위와 오른쪽 위가 크고 확실한 변화입니다.

여기까지 오면 유전자 목록이 나옵니다. 그다음은 경로 분석·기능 해석인데, 그건 다른 강의의 몫이에요.

🤔 Check yourself

Trt_3 의 엑손 비율만 유난히 낮고 인트론이 높았다. 가장 먼저 의심할 것은?

🤔 Check yourself

차등발현 목록을 자를 때 문턱을 두 개 거는 이유는?

Answer 반복이 많으면 1.05배 차이도 유의해집니다. 그래서 padj 와 |log2FC| 를 함께 겁니다.