PyDESeq2: Python 차등 발현 분석

PyDESeq2는 bulk RNA-seq의 조건 사이에서 발현량이 달라진 유전자를 찾는 Python 패키지입니다. 유전자에 배정된 read 수인 raw count와 실험 조건을 적은 샘플 정보(metadata)를 입력받아 시퀀싱 깊이와 생물학적 변동을 모델링하고, 조건 효과의 크기와 통계적 근거를 계산합니다.

예를 들어 실험 대상(subject)마다 Control과 Treatment 샘플이 하나씩 있다면, count 두 개를 직접 나누는 대신 여러 대상에서 반복된 변화를 함께 모델링합니다. 결과는 유전자마다 log2FoldChange, pvalue, padj 같은 열이 있는 표로 나옵니다.

DESeq2와 무엇이 다른가

DESeq2는 R과 Bioconductor에서 사용하는 bulk RNA-seq 차등 발현 분석 패키지입니다. PyDESeq2는 그 방법을 Python 사용자가 이용할 수 있도록 처음부터 다시 구현한 패키지입니다. R이나 DESeq2를 내부에서 호출하는 래퍼가 아닙니다.

두 도구는 정규화, dispersion 추정, 음이항 회귀와 Wald test라는 큰 흐름을 공유합니다. 하지만 PyDESeq2 공식 문서도 재구현 과정과 지원 기능의 차이 때문에 R DESeq2와 결과값이 조금 다를 수 있다고 설명합니다. 분석을 재현하려면 패키지 이름뿐 아니라 버전도 기록해야 합니다.

재현 가능한 분석 프로젝트에서는 다음처럼 사용할 버전을 고정합니다. 아래 예시는 0.5.4를 사용합니다.

uv add "pydeseq2==0.5.4"

입력은 표 두 개다

PyDESeq2에는 count matrix와 metadata가 필요합니다. 둘 다 pandas DataFrame으로 전달할 수 있습니다.

1. count matrix

sample	gene_A	gene_B	gene_C
`S1_Control`	120	8	520
`S1_Treatment`	210	5	490
`S2_Control`	980	12	430
`S2_Treatment`	1,150	7	470

행 하나는 샘플 하나입니다.
열 하나는 유전자 하나입니다.
값은 음수가 아닌 정수 raw count입니다.

HTSeq 파일처럼 유전자가 행이고 샘플이 열인 표는 전치해야 합니다. __no_feature처럼 유전자가 아닌 특수 카운터도 입력 전에 분리합니다. raw count가 만들어지는 과정은 HTSeq raw count에서 설명합니다.

2. metadata

sample	subject	condition
`S1_Control`	S1	Control
`S1_Treatment`	S1	Treatment
`S2_Control`	S2	Control
`S2_Treatment`	S2	Treatment

metadata의 행 이름은 count matrix의 샘플 이름과 정확히 맞아야 합니다. subject와 condition은 통계 모델에서 사용할 샘플 정보입니다. 대상 ID는 숫자의 크기를 비교하는 측정값이 아니라 대상을 구분하는 범주형 변수로 저장합니다.

분석 설계가 비교 질문을 정한다

design = "~ subject + condition"

이 식은 대상마다 다른 발현 기준선을 고려한 뒤 condition 효과를 추정하라는 뜻입니다. condition만 넣을지, batch나 subject도 함께 넣을지는 실험 구조와 분석 질문에 따라 달라집니다. ~와 +를 읽는 방법은 design formula에서 이어집니다.

어느 방향의 차이를 구할지는 contrast로 명시합니다.

contrast = ["condition", "Treatment", "Control"]

이 contrast는 Treatment - Control을 뜻합니다. 따라서 양의 log2FoldChange는 Treatment에서 더 높고, 음수는 Treatment에서 더 낮다는 뜻입니다.

PyDESeq2가 내부에서 하는 일

PyDESeq2의 분석은 다음 순서로 진행됩니다.

size factor 추정: 샘플마다 다른 시퀀싱 깊이를 보정합니다.
dispersion 추정: 같은 조건의 생물학적 반복 사이에서 count가 얼마나 흔들리는지 유전자별로 추정합니다.
모델 적합: raw count를 음이항분포로 모델링하고 design의 각 효과를 추정합니다.
Wald test: contrast로 지정한 조건 효과가 0과 구분될 만큼 큰지 검사합니다.
다중검정 보정: 수천 유전자의 p-value를 함께 보정해 padj를 계산합니다.

raw count를 하나의 정규화된 표로 바꾼 다음 단순 t-test를 반복하는 방식이 아닙니다. 정규화와 대상·조건 효과, 유전자별 변동을 하나의 count model 안에서 다룹니다.

코드에서는 두 객체가 일을 나눈다

from pydeseq2.dds import DeseqDataSet
from pydeseq2.ds import DeseqStats

dds = DeseqDataSet(
    counts=counts,
    metadata=metadata,
    design="~ subject + condition",
    n_cpus=1,
)
dds.deseq2()

stats = DeseqStats(
    dds,
    contrast=["condition", "Treatment", "Control"],
)
stats.summary()

results = stats.results_df

객체	역할
`DeseqDataSet`	count, metadata와 design을 받아 size factor, dispersion과 fold change를 적합
`DeseqStats`	지정한 contrast에 대해 Wald test와 다중검정 보정을 실행하고 결과표 생성

dds.deseq2()가 끝나지 않은 객체로 DeseqStats를 만들면 필요한 모델 추정이 완료되지 않은 상태입니다. 반대로 하나의 적합된 DeseqDataSet에서 서로 다른 contrast를 지정해 여러 비교 결과를 만들 수 있습니다.

결과표의 열 읽기

stats.results_df에는 유전자마다 다음 값이 들어 있습니다.

열	뜻
`baseMean`	모든 샘플에서 size factor로 정규화한 count의 평균
`log2FoldChange`	tested level이 reference level보다 얼마나 높은지 나타내는 log2 변화량
`lfcSE`	추정한 log2 fold change의 표준오차
`stat`	Wald test 통계량
`pvalue`	조건 효과가 없다는 모델 아래에서 현재 결과를 평가한 p-value
`padj`	여러 유전자를 동시에 검사한 것을 보정한 p-value

log2FoldChange = 1이면 tested level이 약 2배, -1이면 약 절반입니다. 방향은 contrast 순서에 따라 바뀌므로 결과를 볼 때 비교 이름을 함께 기록해야 합니다.

후보를 고를 때는 padj만 보지 않습니다. 변화량인 log2FoldChange, 불확실성인 lfcSE, 평균 발현량인 baseMean과 샘플별 패턴을 함께 확인합니다. 여섯 열의 계산 과정과 연결 관계는 DESeq2 결과표 읽기에서, 두 축으로 후보를 좁히는 방법은 볼케이노 플롯 읽기에서 설명합니다. p-value와 다중검정의 일반적인 의미는 통계적 검정과 다중검정을 참고하세요.

DEG 검정용 값과 시각화용 값은 다르다

차등 발현 검정에는 정수 raw count를 입력합니다. PCA나 heatmap에는 큰 count의 분산을 완화한 VST 같은 변환값을 사용합니다.

목적	사용하는 값
DEG 모델 적합	정수 raw count
PCA·샘플 거리·heatmap	VST 또는 다른 시각화용 변환값
결과 해석	`log2FoldChange`, `padj`와 샘플별 발현 패턴

VST 값으로 DEG 검정을 다시 실행하거나 TPM을 PyDESeq2의 count 입력으로 넣지 않습니다. 같은 데이터라도 통계 검정과 시각화가 요구하는 값의 성질이 다르기 때문입니다.

실행 전에 확인할 것

count matrix의 행이 샘플이고 열이 유전자인가
count가 음수가 아닌 정수인가
count와 metadata의 샘플 이름과 순서가 일치하는가
design에 적은 열이 metadata에 존재하는가
subject ID와 batch ID 같은 식별자가 범주형인가
각 비교군에 생물학적 반복이 있는가
저발현 유전자를 정한 기준으로 필터링했는가
contrast의 tested level과 reference level 순서가 의도와 맞는가
사용한 PyDESeq2 버전과 필터 기준을 기록했는가

PyDESeq2는 잘못 구성한 metadata의 생물학적 의미까지 판단해 주지는 않습니다. 계산이 성공해도 design이 연구 질문과 맞지 않으면 다른 질문에 대한 정확한 결과를 만들 수 있습니다.