요약
두 집단(또는 여러 조건) 사이에서 발현량이 유의하게 달라지는 유전자들을 찾는 분석. 즉, 두 집단에 어떤 유전자가 더 많이 발현되는지(up-regulated), 어떤 유전자가 덜 발현되는지(down-regulated)를 찾는 방법. t-test
사전지식
1. library size: 한 샘플에서 얻은 총 read의 수(총 카운트의 수)
2. CPM normaliation (Counts Per Million) = CPM = (gene count / sample total count) * 10^{6}. 10^{6}을 곱하는 이유는 단순히 보기좋게 스케일로 바꾸기 위한 관례
3. RNA-seq: 샘플안에 RNA을 시퀀싱해서, 각 유전자에서 나온 조각(read)이 몇개 관측되었는지 세어내는 방법.
- 세포/조직에서 RNA를 추출
- RNA를 더 다루기 쉬운 cDNA(complementary)로 바꿈 (안정적, PCR증폭쉬움, DNA시퀀싱 플로우적)
- 그 cDNA를 잘게 나눔
- 시퀀서가 그 조각들을 읽어서 수많은 short read를 생성
- 그 read들을 reference genome/transcriptome에 매핑
- 어떤 유전자에 속한 read가 몇 개인지 세어서 gene-by-sample count matrix를 만듦
4. FDR adjusted p-value: 여러 유전자를 한번에 검정할때 생기는 다중비교문제를 보정한 p-value. 예를 들어, p-value가 0.05면 차이가 우연히있다고 할만한 확률이 0.05나됨. 20,000개의 유전자를 0.05하면 1,000개나 되는데, 이정도 false-positive가 나올 수 있어서 FDR(False discovery rate)을 시행. 즉 p-value조정이며, 구체적인 방법은 "Benjamini-Hochberg(BH)"을 이용함. 이는 유의하다고 고른것들중에 false-positive rate을 너무 높지않게하자는거임. p-value을 작은 순서대로 정렬한다음에 p-value가 자신의 순위에 비해 충분히 작은지 보는 방식. 예를 들어, 1등은 p-value가매우작고, N등은 p-value가 좀 덜작아도됨. 총 검정 개수를 m, 정렬된 FDR 기준을 q라고하면, 각 순위 i에 대해서 $\frac{i}{m}q$의 기준을 제시
실전예시: bulk RNA-seq 데이터(airway dataset)
이 데이터셋은 4개의 human airway smooth muscle cell line(세포주)에서 각기 control vs dexamethasone-treated를 짝지은 총 8개 샘플로 구성된 고전적인 공개 예제입니다. 즉, 4개의 샘플에서 2개의 control 및 덱사메타손을 준 샘플입니다. 여러 DESeq2 튜토리얼과 Bioconductor 워크플로에서 반복해서 쓰이는 데이터입니다. 원 출처는 GEO GSE52778 / Himes et al. 2014.
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.stats import ttest_rel
from statsmodels.stats.multitest import multipletests
# -----------------------------
# 1) 공개 데이터 불러오기
# -----------------------------
counts_url = "https://github.com/GuangchuangYu/treedata-workshop/raw/refs/heads/master/_data/airway_scaledcounts.csv"
meta_url = "https://github.com/bioconnector/workshops/raw/refs/heads/master/data/airway_metadata.csv"
anno_url = "https://github.com/bioconnector/workshops/raw/refs/heads/master/data/annotables_grch38.csv"
counts = pd.read_csv(counts_url) # gene x sample
meta = pd.read_csv(meta_url) # sample metadata
anno = pd.read_csv(anno_url) # Ensembl -> symbol annotation
print("counts shape:", counts.shape)
print("meta shape:", meta.shape)
print(meta)
counts shape: (38694, 9)
meta shape: (8, 4)- ensgene은 유전자 ID 컬럼이며, SRR1039508 ~ SRR1039521 : 각 샘플의 count 값 컬럼입니다. 그리고, 각 필드값에는 count(정수)형이 있는데, 각 8개의 샘플에서 측정된 read count을 의미합니다.
# -----------------------------
# 2) count matrix 정리
# -----------------------------
sample_ids = meta["id"].tolist()
mat = counts.set_index("ensgene")[sample_ids]
# library size
lib_sizes = mat.sum(axis=0)
# CPM normalization
cpm = mat.div(lib_sizes, axis=1) * 1e6
# log transform
log_cpm = np.log2(cpm + 1)- CPM바꿔서 샘플별 발현량 평균치를 바꿔주고, log transform을 해줍니다.
# -----------------------------
# 3) 저발현 유전자 필터링
# 예: CPM >= 1 이 2개 이상 샘플에서 관찰
# -----------------------------
keep = (cpm >= 1).sum(axis=1) >= 2
log_cpm_f = log_cpm.loc[keep]
print("tested genes after filtering:", log_cpm_f.shape[0])- CPM 1 이상인경우의 유전자가 최소 2샘플이상 있는 기로 합니다. 이 경우는 굳이 너무 극소하거나 noisy alignment일수 있어 제외합니다.
# -----------------------------
# 4) fold change만 보는 간단한 탐색
# -----------------------------
meta_sorted = meta.sort_values(["celltype", "dex"])
control_ids = meta_sorted.loc[meta_sorted["dex"] == "control", "id"].tolist()
treated_ids = meta_sorted.loc[meta_sorted["dex"] == "treated", "id"].tolist()
fc_df = pd.DataFrame({
"ensgene": cpm.index,
"control_mean_CPM": cpm[control_ids].mean(axis=1).values,
"treated_mean_CPM": cpm[treated_ids].mean(axis=1).values,
})
fc_df = fc_df[(fc_df["control_mean_CPM"] > 1) | (fc_df["treated_mean_CPM"] > 1)].copy()
fc_df["log2FC_treated_vs_control"] = np.log2(
(fc_df["treated_mean_CPM"] + 1) / (fc_df["control_mean_CPM"] + 1)
)
fc_df = fc_df.merge(anno[["ensgene", "symbol"]], on="ensgene", how="left")
top_fc = fc_df.iloc[
fc_df["log2FC_treated_vs_control"].abs().sort_values(ascending=False).index
][["ensgene", "symbol", "control_mean_CPM", "treated_mean_CPM", "log2FC_treated_vs_control"]]
print("\nTop genes by absolute fold change:")
print(top_fc.head(10).to_string(index=False))- log변환 후에 발현량에 차이가 있는 유전자들만 모아서 정렬해봅니다.
# -----------------------------
# 5) paired t-test
# 같은 cell line의 control/treated를 쌍으로 비교
# -----------------------------
control = log_cpm_f[control_ids].to_numpy()
treated = log_cpm_f[treated_ids].to_numpy()
t_stat, pvals = ttest_rel(treated, control, axis=1, nan_policy="omit")
log2fc = treated.mean(axis=1) - control.mean(axis=1)
padj = multipletests(pvals, method="fdr_bh")[1]
res = pd.DataFrame({
"ensgene": log_cpm_f.index,
"mean_log2CPM_control": control.mean(axis=1),
"mean_log2CPM_treated": treated.mean(axis=1),
"log2FC_treated_vs_control": log2fc,
"t_stat": t_stat,
"pvalue": pvals,
"padj": padj,
})
res = res.merge(
anno[["ensgene", "symbol", "description"]],
on="ensgene",
how="left"
).sort_values(["padj", "pvalue"])
print("\nNumber of genes with FDR < 0.05:", (res["padj"] < 0.05).sum())
print("Number of genes with FDR < 0.05 and |log2FC| > 1:",
((res["padj"] < 0.05) & (res["log2FC_treated_vs_control"].abs() > 1)).sum())
top_sig = res[
(res["padj"] < 0.05) &
(res["log2FC_treated_vs_control"].abs() > 1)
][["ensgene", "symbol", "log2FC_treated_vs_control", "pvalue", "padj"]].head(15)
print("\nTop significant genes:")
print(top_sig.to_string(index=False))- paired t-test을 진행합니다. case-control이 동일샘플이기 때문에 paired t-test을 진행합니다.
# -----------------------------
# 6) CRISPLD2 확인
# -----------------------------
cris = res[res["symbol"] == "CRISPLD2"][
["ensgene", "symbol", "mean_log2CPM_control", "mean_log2CPM_treated",
"log2FC_treated_vs_control", "pvalue", "padj"]
]
print("\nCRISPLD2 result:")
print(cris.to_string(index=False))- CRISPLD2라는 대표적인 유전자 확인. 덱사메타손처리하고나면 over-expression되는 대표적인 glucocorticoid-responsive gene임.
# -----------------------------
# 7) Volcano plot
# -----------------------------
plot_df = res.copy()
plot_df["neglog10_padj"] = -np.log10(plot_df["padj"])
sig_mask = (plot_df["padj"] < 0.05) & (plot_df["log2FC_treated_vs_control"].abs() > 1)
plt.figure(figsize=(7, 5))
plt.scatter(
plot_df["log2FC_treated_vs_control"],
plot_df["neglog10_padj"],
s=8, alpha=0.5
)
plt.scatter(
plot_df.loc[sig_mask, "log2FC_treated_vs_control"],
plot_df.loc[sig_mask, "neglog10_padj"],
s=10, alpha=0.8
)
plt.xlabel("log2 fold change (treated vs control)")
plt.ylabel("-log10(FDR)")
plt.title("Airway dataset: DEG example")
plt.tight_layout()
plt.show()
- x축은 치료군/control로 나눈 리드수의 배수를 의미하고(=발현증가/감소), y축은 FDR-adjusted p-value을 의미함. 작을수록유의하니까(예, 0.000001) -log_{10}을쳐서 유의성이클수록 숫자가 커지도록 구함.
# -----------------------------
# 8) CRISPLD2 boxplot
# -----------------------------
gene_id = "ENSG00000103196" # CRISPLD2
gene_vals = log_cpm.loc[gene_id, sample_ids].reset_index()
gene_vals.columns = ["id", "log2CPM"]
gene_vals = gene_vals.merge(meta, on="id")
plt.figure(figsize=(5, 4))
groups = [
gene_vals.loc[gene_vals["dex"] == "control", "log2CPM"],
gene_vals.loc[gene_vals["dex"] == "treated", "log2CPM"]
]
plt.boxplot(groups, tick_labels=["control", "treated"])
x_map = {"control": 1, "treated": 2}
for _, sub in gene_vals.groupby("celltype"):
sub = sub.sort_values("dex")
xs, ys = [], []
for _, row in sub.iterrows():
x = x_map[row["dex"]] + (np.random.rand() - 0.5) * 0.08
y = row["log2CPM"]
plt.plot(x, y, "o")
xs.append(x)
ys.append(y)
plt.plot(xs, ys, "-", alpha=0.6)
plt.ylabel("log2(CPM + 1)")
plt.title("CRISPLD2 (ENSG00000103196)")
plt.tight_layout()
plt.show()
Q. count matrix을 왜 정리해야하나?
A. count matrix의 raw data을 그대로 비교하면 샘플간 공정한 비교가 안될 수 있음. 예들 들어, A샘플은 총 1,000만 reads가 시퀀싱되고, B샘플은 총 2,000만개의 read가 시퀀싱되어있으면 유전자의 카운트가 B가 상대적으로 되게 over-expression된것으로 볼 수 있음.
Q. 왜 log2으로 계산하는지?
RNA-seq값들은 분포가 엄청 치우쳐져있음. 어떤유전자는 0,1,2 수준 이지만, 어떤 유전자는 100 또는 10,000경우도있어서 스케일이 너무큼. 이 상태로 분석을 진행하면 1) 큰 값들의 유전자들이 분석 및 시각화에 엄청 의존적이고, 2) 분산이 평균에 따라서 커지고, 3) t-test가정과 안맞음. 그래서 log(CPM+1)으로 보정함.