pydeseq2로 RNA-seq 차등발현 — R 없이 파이썬만으로 DESeq2 (카운트→볼케이노)
TL;DR — pydeseq2는 RNA-seq 차등발현 분석의 표준인 DESeq2 (R/Bioconductor)를 파이썬으로 다시 구현한 라이브러리입니다.
pip install pydeseq2한 줄이면, R을 깔지 않고도 정수 카운트 행렬과 조건 메타데이터만으로 DESeq2 워크플로 — 크기 인자 정규화·음이항 분산 추정·GLM 적합·Wald 검정·BH 다중검정보정 — 를 그대로 돌릴 수 있습니다. 이 글은 pydeseq2에 동봉된 합성 데이터로DeseqDataSet→.deseq2()→DeseqStats→.summary()까지 한 줄씩 따라가고, 결과표(log2FoldChange·pvalue·padj)로 볼케이노 플롯을 그린 뒤, R DESeq2와 통계·생태계·속도를 비교합니다.흐름은 R DESeq2와 사실상 같습니다. 카운트와 메타데이터로 데이터셋 객체를 만들고(
DeseqDataSet),.deseq2()로 정규화·분산추정·모델적합을 한 번에 끝낸 뒤, 대비(contrast)를 지정한DeseqStats로 검정하고.summary()로 결과표를 얻습니다. 같은 음이항(negative binomial) GLM에 같은 shrinkage 아이디어라, 결과도 R과 거의 일치합니다.재현성을 위해 예제는 어디서 돌려도 결과가 같은 pydeseq2 내장 합성 데이터(100 샘플 × 10 유전자)를 씁니다. 초보가 자주 밟는 지뢰 — 카운트에 소수·음수가 섞인 오류, 샘플 순서 불일치,
design_factors/design인자 이름 변경 — 도 5절에서 실제 메시지와 함께 짚습니다.
🔗 관련 글 · R로 같은 분석을 하는 원조는 DESeq2로 RNA-seq 차등발현(DEG) 분석에서 다뤘고, 이 글은 그 파이썬판입니다 · 결과를 그리는 볼케이노 플롯 시각화와 matplotlib·seaborn 생물 데이터 시각화 · padj가 왜 필요한지는 p-value와 FDR, 다중검정보정 · 카운트를 비교 가능한 값으로 바꾸는 RNA-seq 정규화란?에서 이어집니다
이 글에서 만들 것
차등발현(differential expression) 분석은 RNA-seq에서 가장 자주 하는 일입니다. 두 조건(예: 처리 vs 대조)에서 유전자마다 발현이 통계적으로 유의하게 달라졌는지를 판정하고, 그 결과를 볼케이노 플롯으로 요약하죠. 이 표준을 만든 도구가 R/Bioconductor의 DESeq2 (Love, Huber & Anders 2014)인데, 파이썬으로 파이프라인을 짜는 사람에게는 R을 따로 띄우는 게 늘 걸림돌이었습니다.
pydeseq2 (Muzellec 외, 2023)는 그 DESeq2를 순수 파이썬으로 다시 구현한 라이브러리입니다. 같은 음이항 GLM·같은 분산 shrinkage·같은 Wald 검정을 쓰되, 입출력이 전부 pandas·numpy라 scanpy·scikit-learn 같은 파이썬 생태계와 매끄럽게 이어집니다. 원논문은 pydeseq2 결과가 R DESeq2와 log2 fold change·pvalue에서 거의 일치함을 보고했습니다.
구체적으로 이 글은 pydeseq2 내장 합성 데이터로 ① 카운트 행렬·메타데이터를 준비하고, ② DeseqDataSet을 만들어 .deseq2()로 크기 인자·분산·LFC를 한 번에 적합하고, ③ DeseqStats로 조건 대비를 Wald 검정해 .summary()로 결과표(log2FoldChange·padj)를 얻고, ④ 그 결과로 볼케이노 플롯을 그리고, ⑤ R DESeq2와 비교합니다. 명령은 pydeseq2 0.5.x·파이썬 3.10+ 기준입니다.
# 핵심 흐름 한눈에
from pydeseq2.dds import DeseqDataSet
from pydeseq2.ds import DeseqStats
dds = DeseqDataSet(counts=counts, metadata=meta, design="~condition") # ① 데이터셋 객체
dds.deseq2() # ② 정규화·분산·적합
ds = DeseqStats(dds, contrast=["condition", "B", "A"]) # ③ 대비 지정
ds.summary() # ④ Wald 검정 → 결과표
res = ds.results_df # log2FoldChange · pvalue · padj (⑤ 결과 DataFrame)
1. 사전 준비 — 설치와 예제 데이터
pydeseq2는 pip로 바로 설치됩니다. R이나 Bioconductor는 필요 없고, 의존성은 numpy·pandas·scipy·anndata 정도입니다.
# pip (권장) — 가상환경 안에서
pip install pydeseq2
# 시각화용 (볼케이노 플롯)
pip install matplotlib# 설치 확인 (버전 출력)
import pydeseq2
print(pydeseq2.__version__)0.5.4이제 실습 데이터를 준비합니다. 재현이 쉽도록, pydeseq2가 테스트용으로 동봉한 합성(synthetic) 데이터를 씁니다. load_example_data가 카운트와 메타데이터를 pandas DataFrame으로 바로 돌려줘, 파일을 내려받을 필요가 없습니다.
from pydeseq2.utils import load_example_data
# 카운트 행렬 — 행: 샘플, 열: 유전자 (정수 카운트)
counts = load_example_data(modality="raw_counts", dataset="synthetic")
# 메타데이터 — 행: 샘플, 열: 조건 변수 (condition A/B, group X/Y)
meta = load_example_data(modality="metadata", dataset="synthetic")
print(counts.shape) # (샘플 수, 유전자 수)
print(counts.iloc[:3, :5])(100, 10)
gene1 gene2 gene3 gene4 gene5
sample1 7 1 1 1 11
sample2 2 3 6 2 4
sample3 3 5 2 2 15핵심은 입력 카운트가 정수여야 한다는 점입니다. DESeq2 계열은 음이항 분포를 가정하므로, TPM·FPKM 같은 정규화·소수 값이 아니라 원시 리드 카운트(raw count)를 넣어야 합니다. 정규화는 도구가 내부에서 크기 인자(size factor)로 처리하니, 미리 정규화하지 마세요 — 이 구분은 RNA-seq 정규화 편에서 다뤘습니다.
# 메타데이터 구조 확인 — 각 샘플이 어느 조건인지
print(meta.head()) condition group
sample1 A X
sample2 A Y
sample3 A X
sample4 A Y
sample5 A X메타데이터의 condition 열이 우리가 비교할 조건(A vs B)입니다. pydeseq2는 이 열을 설계식(design formula) ~condition으로 받아, 조건에 따라 샘플을 나눠 발현 차이를 모델링합니다. 실제 데이터라면 이 자리에 counts.csv·metadata.csv를 pd.read_csv(..., index_col=0)로 읽어 넣으면 됩니다. 이때 두 DataFrame의 샘플 인덱스(행 이름)가 같은 순서로 맞아야 합니다.
2. DeseqDataSet — 데이터셋 객체와 모델 적합
DeseqDataSet 만들기
pydeseq2 워크플로의 출발점은 DeseqDataSet 객체입니다. 카운트·메타데이터·설계식을 받아, 이후 정규화와 분산 추정 결과를 모두 이 객체 안에 저장합니다. 내부적으로 anndata (scanpy가 쓰는 그 형식) 위에 얹혀 있어, dds.X(카운트)·dds.obs(샘플 정보)·dds.var(유전자 정보)로 접근할 수 있습니다.
from pydeseq2.dds import DeseqDataSet
from pydeseq2.default_inference import DefaultInference
# 병렬 추론기 — n_cpus로 코어 수 지정(생략 시 가용 코어 전부)
inference = DefaultInference(n_cpus=4)
dds = DeseqDataSet(
counts=counts, # 정수 카운트 (샘플 × 유전자)
metadata=meta, # 샘플별 조건
design="~condition", # 설계식: condition으로 비교
refit_cooks=True, # Cook's 거리 이상치 자동 재적합(권장)
inference=inference,
)
print(dds)AnnData object with n_obs × n_vars = 100 × 10
obs: 'condition', 'group'여기서 설계식 ~condition이 통계 모델의 뼈대입니다. R DESeq2의 design = ~condition과 같은 의미로, "발현을 condition으로 설명하겠다"는 선언이죠. 배치 효과 같은 공변량을 보정하려면 design="~group + condition"처럼 항을 추가합니다(관심 변수를 맨 뒤에 두는 게 관례). 참고로 pydeseq2 0.4.x까지는 design_factors=["condition"]처럼 리스트로 받았는데, 0.5부터 R과 같은 설계식 문자열 design="~condition"으로 바뀌었습니다.
.deseq2() — 정규화·분산추정·적합을 한 번에
# DESeq2 핵심 계산을 한 번에 수행
dds.deseq2().deseq2() 한 줄이 DESeq2의 3대 계산을 순서대로 처리합니다. ① 크기 인자(size factor) 추정 — 샘플별 시퀀싱 깊이 차이를 median-of-ratios로 보정하고, ② 유전자별 분산(dispersion) 추정 — 발현이 비슷한 유전자끼리 정보를 빌려 shrinkage로 안정화하고, ③ 음이항 GLM 적합 — 각 유전자에 log2 fold change를 추정합니다. 계산 결과는 dds 객체 안에 쌓입니다.
# 추정된 크기 인자 (샘플별 정규화 계수) 확인
print(dds.obs["size_factors"].head())sample1 1.038
sample2 0.921
sample3 1.114
sample4 0.972
sample5 1.005
Name: size_factors, dtype: float64크기 인자가 1보다 크면 그 샘플이 상대적으로 깊게 읽혔다는 뜻이라, 카운트를 그만큼 나눠 비교 가능한 수준으로 맞춥니다. 이 값들을 직접 볼 일은 드물지만, 특정 샘플의 크기 인자가 유독 튀면 라이브러리 품질을 의심할 단서가 됩니다. 유전자별 분산 추정값은 dds.var["dispersions"]에 들어 있습니다.
3. DeseqStats — Wald 검정과 결과표
.deseq2()가 모델을 적합했다면, 이제 "조건 A와 B 사이에 유의하게 다른 유전자는?"을 검정할 차례입니다. DeseqStats가 대비(contrast)를 받아 Wald 검정과 다중검정보정을 수행합니다.
contrast 지정과 summary()
from pydeseq2.ds import DeseqStats
# contrast = [변수, 검정 수준, 기준 수준] → B를 A에 견줘 비교
ds = DeseqStats(
dds,
contrast=["condition", "B", "A"], # log2FC > 0 이면 B에서 상향(up)
alpha=0.05, # 유의수준(padj 기준)
inference=inference,
)
ds.summary()Running Wald tests...
Log2 fold change & Wald test p-value: condition B vs A
baseMean log2FoldChange lfcSE stat pvalue padj
gene1 8.981682 0.258287 0.312043 0.827732 0.407843 0.582633
gene2 5.412374 -0.377981 0.375134 -1.007589 0.313653 0.522755
gene3 8.204553 0.129292 0.298734 0.432793 0.665128 0.738987
...contrast=["condition", "B", "A"]는 "B를 A와 비교하라"는 뜻입니다. 그래서 log2FoldChange가 양수면 B에서 발현이 올라간(상향(up)) 유전자, 음수면 내려간(하향(down)) 유전자입니다. 순서를 ["condition", "A", "B"]로 바꾸면 부호가 뒤집힙니다. alpha=0.05는 뒤에서 유의성을 나눌 padj 문턱이고요.
결과표는 ds.results_df에 pandas DataFrame으로 들어 있습니다. 컬럼은 R DESeq2와 이름·의미가 같습니다.
# 결과표를 padj 오름차순으로 정렬해서 보기
res = ds.results_df
print(res.sort_values("padj").head()) baseMean log2FoldChange lfcSE stat pvalue padj
gene8 12.4103 2.8571 0.4126 6.9245 4.37e-12 4.37e-11
gene4 9.8820 -1.9438 0.3852 -5.0462 4.51e-07 2.26e-06
gene2 5.4124 -0.3780 0.3751 -1.0076 0.313653 0.522755
gene1 8.9817 0.2583 0.3120 0.8277 0.407843 0.582633
gene6 7.1055 0.1902 0.3344 0.5688 0.569470 0.632745각 컬럼의 뜻은 이렇습니다. baseMean은 모든 샘플에서 정규화된 평균 발현량, log2FoldChange는 조건 간 발현비의 log2 값(1이면 2배·2면 4배), lfcSE는 그 표준오차, stat은 Wald 통계량, pvalue는 원시 pvalue, 그리고 padj는 다중검정보정(BH, Benjamini-Hochberg)을 거친 조정 pvalue입니다. 유전자 수천 개를 동시에 검정하면 우연히 낮은 pvalue가 쏟아지므로, 반드시 padj로 판단해야 합니다 — 왜 그런지는 p-value와 FDR 편에서 정리했습니다.
lfc_shrink — fold change 축소(선택)
발현이 낮거나 분산이 큰 유전자는 log2 fold change 추정이 요동칩니다. lfc_shrink는 이런 잡음을 줄여, 볼케이노·MA 플롯에서 신뢰할 만한 크기만 남깁니다(R DESeq2의 lfcShrink와 같은 역할).
# LFC 축소 — coeff는 설계식이 만든 계수 이름
ds.lfc_shrink(coeff="condition[T.B]")축소는 pvalue를 바꾸지 않고 log2FoldChange 값만 0쪽으로 당깁니다. 유의성 판정은 그대로 두면서 효과 크기 시각화를 안정화하는 셈이죠. 이제 이 결과표로 볼케이노 플롯을 그립니다.
4. 볼케이노 플롯 — 결과 시각화
볼케이노 플롯(volcano plot)은 차등발현 결과를 한 장으로 요약하는 표준 그림입니다. x축에 log2FoldChange(변화 방향·크기), y축에 −log10(padj)(유의성)를 놓아, 오른쪽·왼쪽 위 구석에 몰린 점이 "크게 변하고 통계적으로도 유의한" 관심 유전자입니다. matplotlib만으로 그릴 수 있습니다.
import numpy as np
import matplotlib.pyplot as plt
res = ds.results_df.dropna(subset=["padj"]) # padj 없는 행 제거
# 유의 기준: padj < 0.05 이고 |log2FC| > 1 (2배 이상 변화)
sig = (res["padj"] < 0.05) & (res["log2FoldChange"].abs() > 1)
up = sig & (res["log2FoldChange"] > 0) # 상향(up)
down = sig & (res["log2FoldChange"] < 0) # 하향(down)
x = res["log2FoldChange"]
y = -np.log10(res["padj"]) # 작을수록 위로
plt.figure(figsize=(6, 5))
plt.scatter(x[~sig], y[~sig], s=12, c="#B8B0A2", label="비유의(ns)") # 회색: 유의하지 않음
plt.scatter(x[up], y[up], s=16, c="#E76F51", label="상향(up)") # 코랄: 상향
plt.scatter(x[down], y[down], s=16, c="#2E7D6B", label="하향(down)") # 그린: 하향
plt.axhline(-np.log10(0.05), ls="--", c="#9AA0A6", lw=1) # padj 0.05 기준선
plt.axvline( 1, ls="--", c="#9AA0A6", lw=1) # log2FC +1
plt.axvline(-1, ls="--", c="#9AA0A6", lw=1) # log2FC -1
plt.xlabel("log2 fold change (B vs A)")
plt.ylabel(r"$-\log_{10}$ adjusted $p$")
plt.legend(frameon=False)
plt.tight_layout()
plt.savefig("volcano.png", dpi=150)# volcano.png 저장됨 — 유의 유전자가 좌우 상단에 강조됨이 코드는 어느 결과표에나 그대로 적용됩니다. 다만 내장 합성 데이터는 유전자가 10개뿐이라 볼케이노가 성글게 나옵니다. 유전자 수천 개짜리 실제 데이터(예: airway·GEO 카운트)를 넣으면, 아래처럼 전형적인 화산 모양이 나타납니다.
# (실데이터 예시) 유의 유전자 개수 세기
n_up = int(up.sum())
n_down = int(down.sum())
print(f"상향 {n_up}개 · 하향 {n_down}개 · 전체 검정 {len(res)}개")상향 214개 · 하향 187개 · 전체 검정 12648개여기서 상향 214·하향 187은 유전자 수천 개짜리 실데이터를 돌렸을 때의 대표 예시입니다. 볼케이노에서 좌우 상단으로 갈수록 "크게, 그리고 확실하게" 달라진 유전자라, 후속 검증(qPCR)이나 경로 분석(pathway analysis)의 1순위 후보가 됩니다. 점 색·기준선·라벨은 볼케이노 시각화 편과 matplotlib 시각화 편의 방식으로 더 다듬을 수 있습니다.
5. 자주 발생하는 에러 & 해결
- 카운트가 정수가 아님 — 가장 흔한 실수 —
DeseqDataSet에 TPM·FPKM처럼 소수·정규화된 값을 넣으면ValueError: counts must be non-negative integers류의 오류가 나거나, 결과가 통계적으로 무의미해집니다. DESeq2 계열은 원시 리드 카운트(raw count)를 전제로 음이항 분포를 적합하므로,featureCounts·htseq-count·salmon의 정수 카운트를 그대로 넣어야 합니다. 정규화는.deseq2()가 내부에서 크기 인자로 처리합니다. - 샘플 순서·이름 불일치 — 카운트의 행(샘플)과 메타데이터의 행이 서로 다른 순서거나 이름이 어긋나면, 조건이 엉뚱한 샘플에 붙어 결과가 통째로 틀립니다.
counts.index.equals(meta.index)로 두 인덱스가 같은지 먼저 확인하고, 어긋나면meta = meta.loc[counts.index]로 맞추세요. pydeseq2는 이 정렬을 자동으로 해 주지 않습니다. design_factorsvsdesign인자 이름 — 오래된 튜토리얼을 따라DeseqDataSet(..., design_factors="condition")으로 쓰면, pydeseq2 0.5+에서는TypeError: unexpected keyword argument 'design_factors'가 납니다. 0.5부터 R과 같은 설계식 문자열design="~condition"으로 바뀌었기 때문입니다. 반대로 0.4.x에design=을 쓰면 인식하지 못하니,pydeseq2.__version__을 먼저 확인하세요.- 저발현 유전자로 인한 padj NaN — 발현이 너무 낮은 유전자는 독립 필터링(independent filtering)·Cook's 이상치 처리 과정에서 padj가
NaN으로 남습니다. 이는 오류가 아니라 "검정에서 제외됨"이라는 신호입니다. 볼케이노를 그리기 전res.dropna(subset=["padj"])로 걸러야 −log10 계산에서 에러가 안 납니다. 사전에counts.loc[:, counts.sum(axis=0) >= 10]처럼 총 카운트가 낮은 유전자를 미리 쳐내면 검정력도 올라갑니다. - 결과가 비었거나 대비 오류 —
contrast=["condition", "B", "A"]에서B·A가 실제 메타데이터의condition값과 정확히 일치해야 합니다. 대소문자·공백만 달라도KeyError가 나거나 빈 결과가 나옵니다.meta["condition"].unique()로 실제 수준(level) 이름을 확인한 뒤 그대로 넣으세요.
6. 확장 — 더 해 볼 것
- 실제 공개 데이터로 돌리기 — 합성 데이터 대신 GEO·recount3의 실제 카운트(예: airway 데이터셋)를
pd.read_csv로 읽어 넣으면, 수천 유전자짜리 진짜 볼케이노가 나옵니다. 전처리는 저발현 유전자 필터링(sum >= 10)만 하면 충분합니다. - 다요인 설계·배치 보정 —
design="~batch + condition"처럼 공변량을 넣으면, 배치 효과를 보정한 채 조건 차이만 뽑을 수 있습니다. 관심 변수를 설계식 맨 뒤에 두는 게 관례입니다. - scanpy·단일세포와 잇기 — pydeseq2는 anndata 기반이라, scanpy로 다룬 단일세포(single-cell) 데이터를 유사벌크(pseudobulk)로 합친 뒤 그대로 넘겨 차등발현을 볼 수 있습니다. R을 거치지 않고 파이썬 한 흐름으로 끝나는 게 최대 장점이죠 — 단일세포 개념은 scRNA-seq 편에서 다뤘습니다.
- 결과 저장·경로 분석 —
res.to_csv("deg.csv")로 결과표를 저장하고, 유의 유전자 목록을 gseapy(GSEA 파이썬판)나 Enrichr로 넘기면 경로 분석까지 파이썬 안에서 이어집니다. - R DESeq2와 교차검증 — 같은 카운트를 R DESeq2와 pydeseq2 양쪽에 돌려 log2FoldChange·padj를 산점도로 비교하면, 두 구현이 거의 일직선으로 일치함을 눈으로 확인할 수 있습니다. 방법론이 궁금하면 R DESeq2 편과 나란히 놓고 보세요.
R DESeq2 vs pydeseq2 — 무엇을 언제 쓸까
| 항목 | R DESeq2 (원조) | pydeseq2 (파이썬) |
|---|---|---|
| 통계 모델 | 음이항 GLM · 분산 shrinkage · Wald | 동일 (재구현) |
| 언어·생태계 | R / Bioconductor | 파이썬(pandas·numpy·scanpy) |
| 입력 | DESeqDataSetFromMatrix | DeseqDataSet (pandas) |
| 결과 | results() 표 | results_df (동일 컬럼) |
| 성숙도 | 2014~ · 사실상 표준 | 2023~ · 빠르게 성숙 중 |
| 강점 | 방대한 문서·확장(apeglm 등) | 파이썬 파이프라인 통합·단일세포 연계 |
결론부터 말하면, 결과는 거의 같으니 생태계로 고르면 됩니다. 파이프라인이 이미 파이썬(scanpy·scikit-learn) 위에 있다면 pydeseq2가 R을 오가는 수고를 없애 줍니다. 반대로 apeglm 축소나 복잡한 상호작용 설계처럼 DESeq2 생태계의 성숙한 확장이 필요하면 R이 아직 앞섭니다. pydeseq2는 원논문에서 R DESeq2와의 일치를 정량적으로 검증했으니, 파이썬 사용자에게는 "R을 안 깔아도 되는 DESeq2"로 충분히 신뢰할 만합니다.
• (R 원조) DESeq2로 RNA-seq 차등발현(DEG) 분석 — 이 글의 R 판, 같은 카운트→볼케이노 흐름을 R로
• (결과 시각화) 볼케이노 플롯 시각화 · matplotlib·seaborn 생물 데이터 시각화 — 결과표를 그림으로
• (통계 기초) p-value와 FDR, 다중검정보정 — padj가 왜 필요한지
• (전처리) RNA-seq 정규화란? — 왜 원시 카운트를 넣고 정규화는 도구에 맡기는지
References
- Love, M. I., Huber, W., & Anders, S. (2014). Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biology, 15(12), 550. (DESeq2 원논문)
- Muzellec, B., Teleńczuk, M., Cabeli, V., & Andreux, M. (2023). PyDESeq2: a python package for bulk RNA-seq differential expression analysis. Bioinformatics, 39(9), btad547. (pydeseq2 원논문)
- Benjamini, Y., & Hochberg, Y. (1995). Controlling the false discovery rate: a practical and powerful approach to multiple testing. Journal of the Royal Statistical Society B, 57(1), 289–300. (BH 다중검정보정 원논문)
- PyDESeq2 documentation. (2024). A simple PyDESeq2 workflow · Step-by-step workflow. pydeseq2.readthedocs.io
- Anders, S., & Huber, W. (2010). Differential expression analysis for sequence count data. Genome Biology, 11(10), R106. (음이항 모델·크기 인자 기초)
Pipette & Pipeline · A bio portfolio journal
'Code Bench > Python for Bio' 카테고리의 다른 글
| gseapy로 GO·GSEA 농축분석 — 파이썬으로 유전자셋 분석 (0) | 2026.07.09 |
|---|---|
| numpy·scipy로 생물 데이터 수치·통계 — 배열 연산·가설검정 (0) | 2026.07.01 |
| scikit-learn으로 생물 데이터 머신러닝 — 분류·교차검증·ROC (0) | 2026.06.29 |
| matplotlib·seaborn으로 생물 데이터 시각화 — 발현 데이터 플롯 실전 (0) | 2026.06.25 |
| scanpy로 단일세포 RNA-seq 분석하기 — QC부터 UMAP·클러스터까지 (Python for Bio) (0) | 2026.06.23 |
