# DESeq2 결과표 읽기

> DESeq2와 PyDESeq2의 유전자별 결과표에서 baseMean, log2FoldChange, lfcSE, stat, pvalue, padj가 무엇을 계산하며 어떻게 연결되는지 설명한다.

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

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

```python
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는 다음과 같은 단순 계산을 사용하지 않습니다.

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

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

```python
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()`를 별도로 호출합니다.

```python
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`로 읽으면 안 됩니다.

:::note[Likelihood ratio test에서는 stat의 뜻이 다르다]
이 문서의 식은 기본 Wald 검정에 해당합니다. Likelihood ratio test(LRT)를 사용하면 `stat`은 축소 모델과 전체 모델의 deviance 차이이며, 카이제곱분포로 p-value를 계산합니다.
:::

## `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를 함께 고려해 이런 [거짓 양성의 누적](/reference/classification-metrics/#통계-검정의-거짓-양성과는-무엇이-다른가)을 다루는 **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 후보를 고르지 않습니다. 다음 순서로 읽으면 관측 규모, 효과 크기와 통계적 근거를 분리할 수 있습니다.

1. **비교 방향 확인**: contrast의 tested level과 reference level은 무엇인가?
2. **`baseMean` 확인**: 충분히 관측됐는가, 지나치게 낮지 않은가?
3. **`log2FoldChange` 확인**: 어느 방향으로 몇 배 달라졌는가?
4. **`lfcSE` 확인**: 효과 크기의 불확실성은 어느 정도인가?
5. **`stat` 확인**: 변화량이 표준오차에 비해 얼마나 큰가?
6. **`padj` 확인**: 다중검정 뒤에도 통계적 근거가 남는가?
7. **샘플별 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](/reference/pydeseq2/)에서, p-value와 FDR의 일반적인 의미는 [통계적 검정과 다중검정](/reference/statistical-testing/)에서 이어집니다. 여러 샘플로 DEG를 검정하는 전체 흐름은 [여러 샘플의 차등발현 검정](/lessons/bulk-rna-deg/)에서 설명합니다.

### 공식 자료

- [PyDESeq2 0.5.4 기본 workflow](https://pydeseq2.readthedocs.io/en/stable/auto_examples/plot_minimal_pydeseq2_pipeline.html)
- [PyDESeq2 0.5.4 DeseqStats API](https://pydeseq2.readthedocs.io/en/stable/api/docstrings/pydeseq2.ds.DeseqStats.html)
- [DESeq2 공식 패키지 매뉴얼](https://bioconductor.org/packages/release/bioc/manuals/DESeq2/man/DESeq2.pdf)
- [DESeq2 원 논문](https://doi.org/10.1186/s13059-014-0550-8)