학습 목표

1. 모델 적합 과정 이해
2. 두 가지 가설 검정 방법 비교 (Wald 검정 vs. LRT)
3. 다중 검정 보정의 중요성 이해
4. 다중 검정 보정에 사용되는 다양한 기법 이해

1. 모델 적합과 가설 검정

DESeq2 분석의 마지막 단계는 각 유전자의 카운트 데이터를 모델에 적합하고 차등 발현 여부를 검정하는 것입니다.

2. 일반화 선형 모델

RNA-seq에서 생성된 카운트 데이터는 과분산(분산 > 평균) 특성을 보이므로, 이를 모델링하기 위한 통계적 분포가 필요합니다. DESeq2는 음이항 분포를 사용하여 RNA-seq 카운트를 다음과 같이 모델링합니다:

K_ij ~ NB(mu_ij, alpha_i)
mu_ij = s_j * q_ij
log2(q_ij) = x_j * beta_i

여기서 필요한 두 가지 주요 매개변수는 크기 인자(size factor)와 분산 추정치(dispersion estimate)입니다. 이후 일반화 선형 모델(GLM)을 사용하여 데이터를 적합합니다. 모델링은 주어진 매개변수 집합下에서 데이터가 어떻게 생성되는지를 수학적으로 형식화하는 과정입니다. 모델 적합 후, 각 샘플 그룹에 대한 계수와 그 표준 오차가 추정됩니다. 계수는 로그2 배수 변화(log2 fold change)의 추정치이며, 이는 가설 검정의 입력값으로 사용됩니다.

3. 가설 검정

가설 검정의 첫 단계는 각 유전자에 대한 귀무 가설(null hypothesis)을 설정하는 것입니다. 일반적인 경우 귀무 가설은 두 샘플 그룹 간에 차등 발현이 없다(LFC == 0)는 것입니다. 그런 다음 통계 검정을 사용하여 관찰된 데이터를 바탕으로 귀무 가설이 참인지 판단합니다.

3.1. Wald 검정

DESeq2에서 Wald 검정은 두 그룹을 비교할 때 기본 가설 검정 방법입니다. Wald 검정은 최대 우도 추정(MLE)을 통해 얻은 매개변수에 대해 수행됩니다. 여기서는 각 유전자의 모델 계수(LFC)를 검정하며, 이 계수는 분산과 같은 매개변수를 사용하여 최대 우도 추정됩니다.

DESeq2는 Wald 검정을 다음 과정으로 수행합니다:

1. LFC를 표준 오차로 나누어 z 통계량을 계산
2. z 통계량을 표준 정규 분포와 비교하여 p-값 계산. 이 p-값은 관측된 값만큼 극단적인 z 통계량이 무작위로 선택될 확률을 나타냄
3. p-값이 작으면 귀무 가설을 기각하고, 유전자가 차등 발현된다는 증거가 있다고 판단

모델 적합과 Wald 검정은 이미 DESeq() 함수의 일부로 실행되었습니다:

# 이전 튜토리얼에서 이미 실행됨
dds <- DESeqDataSetFromTximport(txi, colData = meta, design = ~ sampletype)
dds <- DESeq(dds)

3.2. 우도비 검정

두 개 이상의 샘플 범주를 비교할 때, DESeq2는 Wald 검정 대신 우도비 검정(LRT)을 제공합니다. LRT는 특정 범주에서 다른 범주에 비해 유전자 발현이 증가 또는 감소하는지를 평가하는 대신, 서로 다른 샘플 범주 간에 발현이 변화하는 유전자를 식별합니다.

• Wald 검정과의 차이는?

Wald 검정(기본값)은 각 유전자에 대해 하나의 모델만 추정하고 LFC == 0이라는 귀무 가설을 평가합니다.

우도비 검정은 최대 우도 추정된 매개변수에 대해 수행됩니다. 이 검정에서는 각 유전자에 대해 두 개의 모델을 추정하고, 한 모델의 적합도를 다른 모델과 비교합니다.

• m1은 축소 모델(주요 요인 항이 제거된 설계 공식)
• m2는 완전 모델(dds 객체 생성 시 제공한 전체 설계 공식)

우리는 완전 모델이 축소 모델만큼 적합하다는 귀무 가설을 평가합니다. 귀무 가설을 기각하면, 완전 모델(우리가 관심 있는 주요 요인)이 상당한 변동을 설명하며, 따라서 해당 유전자가 여러 수준에서 차등 발현됨을 의미합니다. DESeq2는 편차 분석(ANODEV)을 사용하여 두 모델 적합을 비교함으로써 LRT를 구현합니다. 결과적으로 LR은 카이제곱 분포를 따르며, 이를 사용하여 관련 p-값을 계산합니다.

LRT를 사용하려면 DESeq() 함수에 두 개의 추가 매개변수를 지정합니다:

1. LRT 검정 사용 지정
2. "축소된" 모델

# 우도비 검정
dds_lrt <- DESeq(dds, test="LRT", reduced = ~ 1)

우리의 "완전" 모델에는 하나의 요인(샘플 유형)만 있으므로, "축소" 모델(해당 요인 제거)에는 설계 공식에 아무것도 남지 않습니다. DESeq2는 설계 공식에 아무것도 없는 모델을 적합할 수 없으므로, 절편은 ~ 1 구문을 사용하여 모델링됩니다.

4. 다중 검정 보정

Wald 검정 또는 LRT 중 어떤 것을 사용하든, 각 검정된 유전자는 p-값과 연결됩니다. 이 결과를 사용하여 어떤 유전자가 통계적으로 유의하게 차등 발현되는지 결정합니다. 그러나 p-값을 직접 사용할 수는 없습니다.

4.1. p-값의 문제

유의성 임계값 p < 0.05를 가진 유전자는 5%의 확률로 위양성(false positive)일 가능성이 있습니다. 예를 들어, 20,000개의 유전자에 대해 차등 발현을 검정하는 경우 p < 0.05에서는 우연히 1,000개의 유전자가 위양성으로 나타날 것으로 예상됩니다. 총 3,000개의 유전자가 차등 발현된다면, 약 1/3의 유전자가 위양성일 수 있습니다!

각 p-값은 단일 검정(단일 유전자)의 결과이기 때문입니다. 검정하는 유전자가 많을수록 위양성률은 증가합니다. 이것이 다중 검정 문제입니다.

4.2. 보정 방법

다중 검정에는 몇 가지 일반적인 보정 방법이 있습니다:

• 본페로니(Bonferroni) 보정: 조정된 p-값 = p-값 × m (m = 총 검정 수). 매우 보수적인 방법으로, 위음성(false negative) 가능성이 높아 일반적으로 권장되지 않습니다.
• FDR/벤자미니-호흐버그(Benjamini-Hochberg) 방법: Benjamini와 Hochberg(1995)는 위발견율(FDR) 개념을 정의하고, 독립적인 p-값 목록이 주어졌을 때 예상 FDR을 지정된 수준 이하로 제어하는 알고리즘을 개발했습니다.
• Q-값/스토리(Storey) 방법: 특성을 유의하다고 판단할 때 달성 가능한 최소 FDR. 예를 들어, 유전자 X의 q-값이 0.013이면, 유전자 X의 p-값만큼 작은 p-값을 가진 유전자 중 1.3%가 위양성임을 의미합니다.

DESeq2는 검정 전에 유의미한 차등 발현 가능성이 낮은 유전자(예: 낮은 카운트 또는 이상치 샘플을 가진 유전자)를 제거하여 검정할 유전자 수를 줄이는 데 도움을 줍니다. 또한, 벤자미니-호흐버그 절차를 사용하여 위발견율을 낮추는 다중 검정 보정이 구현됩니다.

4.3. FDR < 0.05의 의미

FDR 임계값을 < 0.05로 설정하면, 차등 발현 유전자 중 위양성 비율이 5%일 것으로 예상한다는 의미입니다. 예를 들어, 500개의 유전자를 차등 발현이라고 판단하고 FDR 임계값이 0.05인 경우, 그중 약 25개가 위양성일 것으로 예상됩니다.