DESeq2 결과표 읽기

DESeq2와 PyDESeq2의 차등 발현 결과표는 유전자마다 두 조건의 발현 차이와 그 불확실성을 계산한 통계 결과입니다. 정수 raw count와 샘플 정보를 입력받아, 각 유전자를 한 행으로 두고 baseMean, log2FoldChange, lfcSE, stat, pvalue, padj 여섯 값을 출력합니다.

이 문서는 기본적인 Wald 검정 결과를 설명합니다. 예시는 Treatment를 Control과 비교하는 다음 방향을 사용합니다.

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

따라서 양의 log2FoldChange는 Treatment에서 높다는 뜻이고, 음수는 Control에서 높다는 뜻입니다. 비교 방향이 바뀌면 부호도 바뀝니다.

먼저 결과표 한 행을 읽어 보자

아래 값은 설명을 위해 만든 가상 결과입니다.

gene	baseMean	log2FoldChange	lfcSE	stat	pvalue	padj
`gene_A`	500	2.00	0.50	4.00	$6.3 \times 10^{-5}$	0.003

이 한 행은 다음처럼 읽습니다.

gene_A는 전체 샘플에서 정규화 count가 평균 500 정도 관측됐습니다. Treatment에서 Control보다 약 4배 높다고 추정됐고, 변화량 2.00의 표준오차는 0.50입니다. 변화량을 표준오차로 나눈 Wald 통계량은 4.00이며, 유전자 하나의 p-value는 약 0.000063입니다. 모든 유전자를 함께 검사한 것을 보정한 p-value는 0.003입니다.

여섯 값은 서로 독립적인 이름 목록이 아닙니다. 관계를 먼저 나누면 표가 쉬워집니다.

역할	열	다음 값과의 관계
전체 관측 규모 요약	`baseMean`	다른 다섯 열을 직접 계산하는 출발값이 아님
조건 효과의 크기	`log2FoldChange`	`lfcSE`와 함께 `stat` 계산에 사용
효과 추정의 불확실성	`lfcSE`	`log2FoldChange`와 함께 `stat` 계산에 사용
효과와 불확실성의 비율	`stat`	`pvalue` 계산에 사용
유전자 하나의 검정 결과	`pvalue`	모든 유전자의 값과 함께 보정
다중검정 보정 결과	`padj`	후보 선택에 사용하는 통계적 근거

핵심은 baseMean과 log2FoldChange가 입력과 출력 관계가 아니라는 점입니다. 둘 다 같은 count 데이터에서 계산되지만 서로 다른 질문에 답합니다.

`baseMean`: 이 유전자가 전반적으로 얼마나 관측됐나

baseMean은 모든 샘플의 정규화 count 평균입니다. 샘플마다 시퀀싱 깊이가 다르므로 raw count를 그대로 평균하지 않고, 각 샘플의 size factor로 나눈 값을 사용합니다.

유전자 $i$ , 샘플 $j$ 의 raw count를 $K_{ij}$ , 샘플의 size factor를 $s_j$ , 전체 샘플 수를 $m$ 이라고 하면 다음처럼 이해할 수 있습니다.

\operatorname{baseMean}_i = \frac{1}{m} \sum_{j=1}^{m} \frac{K_{ij}}{s_j}

baseMean = 500은 특정 샘플에서 read가 500개였다는 뜻이 아닙니다. 모든 조건의 샘플을 합쳐 시퀀싱 깊이를 보정한 뒤 평균한 관측 규모가 약 500이라는 뜻입니다.

`baseMean`만으로 방향은 알 수 없다

다음 두 유전자는 baseMean이 같아도 조건 차이는 전혀 다를 수 있습니다. 이해를 위한 단순 평균 예시입니다.

유전자	Control 평균	Treatment 평균	전체 평균	단순 fold change
`gene_A`	50	50	50	1배
`gene_B`	10	90	50	9배

두 유전자의 전체 평균은 모두 50입니다. 하지만 gene_A는 조건 차이가 없고 gene_B는 Treatment에서 높습니다. 따라서 baseMean으로 log2FoldChange를 계산하거나 방향을 판단할 수 없습니다.

baseMean은 다음 질문에 답할 때 유용합니다.

이 유전자가 전체적으로 충분히 관측됐는가?
매우 낮은 count에서 큰 배수 변화가 나온 것은 아닌가?
MA plot의 가로축에서 어느 발현 구간에 있는가?

반대로 baseMean은 Treatment 평균이나 Control 평균을 따로 보여 주지 않습니다. 조건별 값과 환자별 패턴은 정규화 count를 별도로 확인해야 합니다.

`log2FoldChange`: 어느 방향으로 몇 배 달라졌나

log2FoldChange는 비교하려는 조건 효과를 로그 2 단위로 표현한 값입니다. 부호로 방향을 보고, 절댓값으로 변화의 크기를 읽습니다.

log2FoldChange	Treatment / Control	해석
`2`	4	Treatment에서 약 4배 높음
`1`	2	Treatment에서 약 2배 높음
`0`	1	차이 없음
`-1`	1/2	Treatment가 Control의 약 절반
`-2`	1/4	Treatment가 Control의 약 4분의 1

값을 원래 배수로 되돌릴 때는 다음 식을 사용합니다.

\text{fold change}=2^{\text{log2FoldChange}}

가상 결과의 log2FoldChange = 2는 $2^2=4$ 이므로 Treatment에서 약 4배 높다는 뜻입니다.

왜 로그 2로 표현하나

배수 변화는 원래 비대칭입니다. 2배 증가는 2지만 2배 감소는 0.5입니다. 로그 2로 바꾸면 같은 크기의 증가와 감소가 1과 -1이 되어 0의 양쪽에서 대칭적으로 비교할 수 있습니다.

로그는 곱셈 관계를 덧셈으로 바꾸기도 합니다.

\log_2(a \times b)=\log_2(a)+\log_2(b)

시퀀싱 깊이, 환자별 기준선과 조건 효과처럼 발현량에 곱으로 작용하는 요소를 통계 모델 안에서는 더하는 항으로 다룰 수 있습니다. 밑이 반드시 2여야 하는 것은 아니지만, 1을 2배로 바로 읽을 수 있어 결과표에는 로그 2가 편리합니다.

두 그룹 평균을 직접 나눈 값이 아니다

두 숫자만 비교한다면 다음 식으로 log2 fold change를 구할 수 있습니다.

\log_2\left( \frac{\text{Treatment 발현량}} {\text{Control 발현량}} \right)

하지만 DESeq2와 PyDESeq2는 다음과 같은 단순 계산을 사용하지 않습니다.

log2((treatment_mean + 1) / (control_mean + 1))

유전자마다 모든 샘플의 정수 raw count를 음이항분포 모델에 넣고, design formula에 적은 조건 효과를 추정합니다. paired design이 다음과 같다면:

design = "~ patient + condition"

환자마다 다른 기준선을 허용한 상태에서 여러 환자에게 반복되는 condition 효과를 구합니다. 그 모델 계수를 로그 2 단위로 나타낸 것이 log2FoldChange입니다. baseMean은 이 모델 계수의 계산 재료가 아니라 별도로 제공되는 관측 규모 요약입니다.

음이항분포 모델의 식

유전자 $i$ , 샘플 $j$ 의 count $K_{ij}$ 는 평균 $\mu_{ij}$ 와 유전자별 dispersion $\alpha_i$ 를 가진 음이항분포로 모델링합니다.

K_{ij} \sim \operatorname{NB}(\mu_{ij}, \alpha_i)

예상 count는 샘플의 size factor $s_j$ 와 유전자의 조건별 발현 수준 $q_{ij}$ 로 나눌 수 있습니다.

\mu_{ij}=s_jq_{ij}

로그를 취하면 곱셈이 덧셈으로 바뀝니다.

\log(\mu_{ij}) = \log(s_j) + \beta_{i0} + x_{j1}\beta_{i1} + \cdots

$\log(s_j)$ 는 샘플별 시퀀싱 깊이를 반영하는 항입니다. $x_{j1}$ 같은 값은 metadata에서 만든 design matrix의 열이고, $\beta_{i1}$ 은 해당 열의 효과입니다.

condition이 Control과 Treatment 두 범주이고 Control이 기준 범주라면, condition 계수가 Treatment 대 Control의 로그 fold change가 됩니다. 모델 내부 계산은 자연로그를 사용할 수 있지만 DESeq2와 PyDESeq2 결과표는 계수와 표준오차를 로그 2 단위로 제공합니다.

LFC shrinkage를 적용했는지 확인한다

낮은 count나 큰 dispersion을 가진 유전자는 log2 fold change가 불안정하게 커질 수 있습니다. LFC shrinkage는 정보가 부족한 유전자의 변화량을 0 쪽으로 줄여 유전자 순위와 시각화를 안정화하는 선택 단계입니다.

PyDESeq2 0.5.4에서는 lfc_shrink()를 별도로 호출합니다.

stats.lfc_shrink(coeff="condition[T.Treatment]")

이 호출 전후에는 log2FoldChange와 lfcSE가 달라질 수 있습니다. 공식 예제에서는 기존 Wald 검정의 stat, pvalue, padj는 유지됩니다. 결과를 재현하려면 shrinkage 적용 여부와 사용한 계수를 기록해야 합니다.

`lfcSE`: log2 fold change를 얼마나 정확하게 추정했나

lfcSE는 log2 fold change standard error, 즉 log2 fold change의 표준오차입니다. log2FoldChange가 추정한 조건 효과라면 lfcSE는 그 추정값이 얼마나 불확실한지를 나타냅니다.

표준오차는 표준편차와 다릅니다.

값	답하는 질문
표준편차(SD)	개별 샘플 값들이 얼마나 퍼져 있는가?
표준오차(SE)	그 샘플들로 계산한 효과 추정값이 얼마나 불확실한가?

단순한 평균에서는 $SE=SD/\sqrt{n}$ 관계가 있지만, lfcSE는 이 공식을 count에 바로 적용한 값이 아닙니다. 음이항분포 모델, 유전자별 dispersion, size factor와 design matrix를 함께 고려해 조건 계수의 표준오차를 계산합니다.

일반적으로 다음 조건은 lfcSE를 크게 만들 수 있습니다.

샘플 수가 적다.
같은 조건 안에서 count가 크게 흔들린다.
유전자의 count가 매우 낮다.
design에 추정해야 할 효과가 많다.
두 효과가 데이터에서 서로 잘 분리되지 않는다.

반대로 반복 샘플이 많고 같은 방향의 변화가 안정적으로 관찰되면 표준오차가 작아지는 경향이 있습니다.

변화량과 비교해서 읽는다

lfcSE에는 모든 분석에 통하는 단독 기준이 없습니다. 0.5보다 작으면 좋다처럼 읽지 않고 log2FoldChange와 비교합니다.

가상 결과에서는 log2FoldChange = 2.00, lfcSE = 0.50입니다. 변화량이 표준오차의 4배입니다.

2.00 \div 0.50=4.00

Wald 근사를 사용하면 대략적인 95% 구간은 다음처럼 확인할 수 있습니다.

2.00 \pm 1.96 \times 0.50 = 1.02 \sim 2.98

배수로 바꾸면 약 $2^{1.02}=2.0$ 배에서 $2^{2.98}=7.9$ 배입니다. 정확한 배수에는 불확실성이 있지만 구간 전체가 0보다 크므로 Treatment에서 증가했다는 방향은 비교적 안정적입니다.

이 구간은 기본 Wald 근사를 이해하기 위한 값입니다. shrinkage나 다른 검정 설정을 사용했다면 해당 방법이 제공하는 추론 절차를 확인해야 합니다.

`stat`: 변화량이 표준오차의 몇 배인가

기본 Wald 검정에서 stat은 다음과 같이 계산합니다.

\operatorname{stat} = \frac{ ext{log2FoldChange}-0}{\text{lfcSE}}

분자의 0은 조건 효과가 없다, 즉 log2FoldChange = 0이라는 귀무가설입니다. 가상 결과에서는 2.00 ÷ 0.50 = 4.00이므로 stat = 4.00입니다.

큰 양수: Treatment에서 증가했고 0에서 멀리 떨어져 있음
큰 음수: Treatment에서 감소했고 0에서 멀리 떨어져 있음
0에 가까운 값: 추정 효과가 작거나 불확실성이 큼

stat의 절댓값이 크면 차이가 없다는 가설과 데이터가 잘 맞지 않는다는 통계적 근거가 강해집니다. 하지만 생물학적으로 중요한 유전자라는 뜻은 아닙니다. 작은 변화도 표준오차가 매우 작으면 큰 stat을 가질 수 있습니다.

lfc_null을 0이 아닌 값으로 지정하면 검정 기준도 바뀝니다. 예를 들어 단순히 차이가 있는지가 아니라 |log2FoldChange| > 1인지 검사할 수 있습니다. 이때는 stat을 항상 log2FoldChange ÷ lfcSE로 읽으면 안 됩니다.

`pvalue`: 차이가 없다는 모델과 얼마나 맞지 않나

기본 양측 Wald 검정의 귀무가설은 다음과 같습니다.

H_0:\text{log2FoldChange}=0

pvalue는 이 귀무가설이 맞다고 가정했을 때, 현재 stat만큼 또는 그보다 극단적인 결과가 관찰될 확률입니다. Wald 통계량을 표준정규분포와 비교해 양쪽 꼬리의 확률을 계산합니다.

p = 2\left(1-\Phi\left(|\operatorname{stat}|\right)\right)

$\Phi$ 는 표준정규분포의 누적분포함수입니다. 가상 결과의 stat = 4.00은 약 $6.3 \times 10^{-5}$ 의 p-value가 됩니다.

CSV에서 6.3E-05처럼 보이면 $6.3 \times 10^{-5}$ 를 줄여 쓴 과학적 표기법입니다. 소수로는 0.000063입니다.

p-value가 뜻하지 않는 것

pvalue = 0.01은 다음 뜻이 아닙니다.

결과가 우연일 확률이 1%다.
귀무가설이 참일 확률이 1%다.
이 유전자가 원인일 확률이 99%다.
변화량이 크거나 생물학적으로 중요하다.

p-value는 조건 효과가 없다는 특정 통계 모델 아래에서 데이터가 얼마나 극단적인지를 나타냅니다. 모델과 design이 연구 질문에 맞는지, count 품질과 이상치가 괜찮은지는 별도로 확인해야 합니다.

`padj`: 수천 유전자를 검사한 영향을 보정한다

RNA-seq에서는 유전자 하나가 아니라 수천 개를 동시에 검사합니다. 실제로 모든 유전자에 차이가 없다고 가정한 단순한 장면에서 20,000개 유전자에 pvalue < 0.05만 적용하면 평균적으로 약 1,000개가 우연히 기준을 통과할 수 있습니다.

padj는 여러 p-value를 함께 고려해 이런 거짓 양성의 누적을 다루는 adjusted p-value입니다. 기본적으로 Benjamini-Hochberg 방식으로 거짓 발견율(false discovery rate, FDR)을 제어합니다. 이때의 거짓 발견은 분류 혼동행렬의 FP 한 건과 구분해야 합니다.

FDR은 선택된 결과 집합 안에서 기대되는 거짓 발견의 비율을 제한합니다. padj < 0.05는 각 유전자마다 틀릴 확률이 정확히 5%라는 뜻이 아닙니다.

`padj`는 한 행만으로 계산할 수 없다

pvalue는 해당 유전자의 stat에서 계산하지만, padj는 분석에 포함된 모든 유전자의 p-value 순위에 영향을 받습니다. 같은 유전자의 p-value가 같아도 함께 검사한 유전자 집합이나 필터 설정이 달라지면 padj가 달라질 수 있습니다.

그래서 사전 필터 기준과 분석에 들어간 유전자 수를 기록해야 합니다. 결과를 본 뒤 마음에 들지 않는 유전자를 임의로 제외하고 padj를 다시 계산하면 선택 과정 자체가 결과에 영향을 줄 수 있습니다.

independent filtering이 함께 적용될 수 있다

PyDESeq2 0.5.4의 DeseqStats는 기본적으로 independent filtering을 사용합니다. baseMean이 매우 낮아 검정력이 부족한 유전자를 다중검정 보정 대상에서 제외하면서, 설정한 alpha에서 발견 수를 늘릴 수 있는 cutoff를 찾습니다.

이 경우 pvalue는 있지만 padj가 NaN인 행이 생길 수 있습니다. 이것은 padj = 1과 같지 않습니다. 해당 유전자가 독립 필터를 통과하지 않아 보정값이 계산되지 않았다는 뜻입니다.

`NaN`은 어떤 경우에 생기나

빈 값은 모두 같은 이유로 생기지 않습니다. PyDESeq2 설정과 실행 로그를 함께 확인해야 합니다.

보이는 결과	가능한 이유	먼저 확인할 것
`pvalue`와 `padj`가 모두 `NaN`	Cook’s distance로 검출된 count 이상치	일부 샘플 하나가 결과를 지배하는가
`pvalue`는 있고 `padj`만 `NaN`	independent filtering에서 낮은 `baseMean` 행 제외	독립 필터 설정과 cutoff
모든 count가 0	효과와 분산을 추정할 정보가 없음	사전 필터와 원본 count
매우 큰 LFC와 큰 SE	한 조건의 낮은 count로 비율이 불안정함	조건별·샘플별 count

NaN을 일괄적으로 0이나 1로 바꾸기 전에 왜 비어 있는지 확인합니다. 그래프에서 편의를 위해 padj 결측을 1로 취급했다면, 통계 결과를 변경한 것이 아니라 시각화 분류를 위한 처리였다고 명시해야 합니다.

여섯 값을 함께 읽는 순서

한 열만 보고 DEG 후보를 고르지 않습니다. 다음 순서로 읽으면 관측 규모, 효과 크기와 통계적 근거를 분리할 수 있습니다.

비교 방향 확인: contrast의 tested level과 reference level은 무엇인가?
baseMean 확인: 충분히 관측됐는가, 지나치게 낮지 않은가?
log2FoldChange 확인: 어느 방향으로 몇 배 달라졌는가?
lfcSE 확인: 효과 크기의 불확실성은 어느 정도인가?
stat 확인: 변화량이 표준오차에 비해 얼마나 큰가?
padj 확인: 다중검정 뒤에도 통계적 근거가 남는가?
샘플별 count 확인: 일부 샘플이나 한 환자가 결과를 만들지는 않았는가?

자주 만나는 조합은 다음처럼 해석할 수 있습니다.

결과 조합	해석
큰 LFC, 작은 SE, 작은 `padj`	크고 안정적으로 추정된 조건 차이
작은 LFC, 작은 SE, 작은 `padj`	작지만 반복 샘플에서 일관된 조건 차이
큰 LFC, 큰 SE, 큰 `padj`	변화량은 커 보이지만 불확실함
낮은 `baseMean`, 큰 LFC	read 몇 개 차이가 큰 비율이 됐는지 확인 필요
작은 `padj`, 일부 샘플만 극단적	이상치와 모델 적합을 다시 확인해야 함

후보 선택에서 padj < 0.05, |log2FoldChange| >= 1 같은 기준을 사용할 수 있지만 보편적인 생물학적 경계는 아닙니다. 연구 질문, 표본 수, 후속 검증 비용과 필요한 효과 크기에 맞춰 정합니다.

결과표만으로 말할 수 없는 것

여섯 값은 조건 차이의 크기와 통계적 근거를 보여 주지만 다음 결론을 직접 증명하지 않습니다.

해당 유전자가 질환의 원인이다.
발현 변화가 단백질 양이나 활성을 그대로 바꿨다.
bulk 조직(tissue)의 변화가 특정 세포 종류 안에서 일어났다.
다른 데이터셋이나 인구집단에서도 같은 효과가 재현된다.
padj가 가장 작은 유전자가 생물학적으로 가장 중요하다.

bulk RNA-seq의 DEG는 세포 안의 조절 변화와 조직의 세포 구성 변화를 함께 반영할 수 있습니다. 상위 후보는 환자별 정규화 count, 원시 read 품질, 조직 정보와 독립 데이터에서 다시 확인해야 합니다.

재현할 때 함께 기록할 것

같은 결과표를 다시 만들려면 CSV만 저장하지 말고 다음 조건도 기록합니다.

DESeq2 또는 PyDESeq2의 정확한 버전
raw count를 만든 annotation과 정량 도구 버전
사전 저발현 필터 기준
design formula
contrast의 tested level과 reference level 순서
Wald test 또는 LRT 선택
alpha, lfc_null, alt_hypothesis
independent filtering과 Cook’s filtering 설정
LFC shrinkage 적용 여부와 계수
결과를 만든 전체 샘플과 제외한 샘플

도구의 입력과 실행 흐름은 PyDESeq2에서, p-value와 FDR의 일반적인 의미는 통계적 검정과 다중검정에서 이어집니다. 여러 샘플로 DEG를 검정하는 전체 흐름은 여러 샘플의 차등발현 검정에서 설명합니다.

DESeq2 결과표 읽기

먼저 결과표 한 행을 읽어 보자

baseMean: 이 유전자가 전반적으로 얼마나 관측됐나

baseMean만으로 방향은 알 수 없다

log2FoldChange: 어느 방향으로 몇 배 달라졌나