DIA-NN 파이프라인 코드 리뷰 — FDR 보정 안 하면 생기는 일

Code review session with statistical charts

코드 리뷰를 하게 된 배경

BioAI Market에 DIA-NN 분석 파이프라인을 통합하면서, 기존에 사용하던 R 코드를 전수 리뷰할 기회가 있었다. 코드 위치는 /Users/kang/work2026/dia_results/r_code/였고, DIA-NN의 MaxLFQ 정량 결과를 받아서 전처리, 정규화, DE 분석, 시각화까지 수행하는 파이프라인이었다.

전반적으로 4/5 stars 평가를 줄 만큼 잘 짜여진 코드였다. 변수명이 명확하고, 각 단계가 함수로 분리되어 있었으며, 주석도 적절했다. 하지만 한 줄이 모든 것을 망칠 수 있었다.

핵심 버그: DEP_USE_ADJUSTED_P <- FALSE

코드 상단의 설정 변수들 중에 이런 게 있었다:

# === Configuration ===
DEP_PVALUE_THRESHOLD <- 0.05
DEP_LOG2FC_THRESHOLD <- 1.0
DEP_USE_ADJUSTED_P <- FALSE  # ← 이게 문제

DEP_USE_ADJUSTED_P <- FALSE라는 건, DE 분석 결과에서 raw p-value를 기준으로 유의미한 단백질을 필터링한다는 뜻이다. adjusted p-value (FDR-corrected p-value)를 무시하는 것이다.

이게 왜 문제인지, 숫자로 보여주겠다.

FDR 보정이 왜 필수인가

프로테오믹스 DE 분석에서는 보통 수천 개의 단백질을 동시에 검정한다. 가령 5,000개 단백질에 대해 각각 t-test를 수행한다고 하자.

p < 0.05를 기준으로 하면, 단 한 개도 실제로 차이가 없더라도 기대되는 거짓 양성 수는:

5,000 × 0.05 = 250개

5,000개 중 250개가 "유의미하다"고 나오지만, 이 중 상당수가 순전히 우연에 의한 결과다. 이게 **다중 검정 문제(multiple testing problem)**다.

BH(Benjamini-Hochberg) Correction

BH correction은 False Discovery Rate를 제어하는 가장 널리 쓰이는 방법이다. 원논문(Benjamini & Hochberg, 1995)에서 제안되었다.

알고리즘은 직관적이다:

# BH correction 원리
# 1. p-value를 오름차순 정렬
# 2. 각 p-value에 보정 계수 적용
# 순위(i) / 전체 검정 수(m) × 목표 FDR(q)

p_sorted <- sort(p_values)
m <- length(p_values)
for (i in 1:m) {
  adjusted_p[i] <- p_sorted[i] * m / i
}

R에서는 p.adjust(p_values, method = "BH")로 한 줄이면 된다. 이미 파이프라인에서 adjusted p-value를 계산은 하고 있었다. 다만 사용하지 않았을 뿐이다.

실제 데이터에서의 영향

같은 데이터셋에 raw p vs adjusted p를 적용했을 때의 차이:

# Raw p-value < 0.05 기준
dep_raw <- results %>% 
  filter(p_value < 0.05, abs(log2FC) > 1)
nrow(dep_raw)  # 487개

# Adjusted p-value < 0.05 기준  
dep_adj <- results %>% 
  filter(adj_p_value < 0.05, abs(log2FC) > 1)
nrow(dep_adj)  # 127개

487개 vs 127개. 약 3.8배 차이다. raw p-value 기준으로는 487개가 "유의미"하다고 나오지만, FDR 보정을 하면 127개로 줄어든다. 나머지 360개는 거짓 양성일 가능성이 높은 것들이다.

이 360개의 "가짜" DEP가 downstream 분석에 포함되면:

바이오마커 검증에서 false positive 증가
Pathway enrichment 분석이 왜곡
논문에 틀린 결과가 실림

Volcano Plot으로 시각화

library(ggplot2)

volcano_raw <- ggplot(results, aes(x = log2FC, y = -log10(p_value))) +
  geom_point(aes(color = significant_raw), alpha = 0.5) +
  geom_hline(yintercept = -log10(0.05), linetype = "dashed", color = "red") +
  ggtitle("Raw p-value (487 DEPs)") +
  theme_minimal()

volcano_adj <- ggplot(results, aes(x = log2FC, y = -log10(adj_p_value))) +
  geom_point(aes(color = significant_adj), alpha = 0.5) +
  geom_hline(yintercept = -log10(0.05), linetype = "dashed", color = "red") +
  ggtitle("Adjusted p-value (127 DEPs)") +
  theme_minimal()

Volcano plot을 나란히 놓으면 차이가 극명하다. raw p 기준 plot에서는 점들이 빽빽하게 유의미 영역에 몰려 있고, adjusted p 기준에서는 확실히 유의미한 것들만 남는다.

권장 수정사항

가장 중요한 수정:

# 수정 전
DEP_USE_ADJUSTED_P <- FALSE

# 수정 후
DEP_USE_ADJUSTED_P <- TRUE  # FDR correction 필수 적용

이 한 줄이 분석 결과의 신뢰도를 결정한다.

기타 개선점

코드 리뷰에서 발견한 다른 개선사항들:

1. 재현성을 위한 seed 설정

# 추가 필요
set.seed(42)  # imputation, bootstrapping 등에서 재현성 보장

2. 시각화 코드 정리

# 기존: 각 plot마다 theme 설정 반복
# 개선: 공통 theme 함수화
my_theme <- function() {
  theme_minimal() +
    theme(
      plot.title = element_text(size = 14, face = "bold"),
      axis.text = element_text(size = 10),
      legend.position = "bottom"
    )
}

3. 로깅 추가

# 각 단계별 처리 결과 로깅
cat(sprintf("[%s] Filtered %d proteins → %d significant DEPs\n",
    Sys.time(), nrow(results), nrow(dep_filtered)))

교훈: 설정 하나가 연구 전체를 바꾼다

이번 코드 리뷰에서 얻은 가장 큰 교훈은, 통계적 설정 하나가 연구 결론을 완전히 바꿀 수 있다는 것이다. 코드 자체는 깔끔했지만, FALSE 하나가 수백 개의 거짓 양성을 만들어내고 있었다.

BioAI Market에서는 이 경험을 반영해서, DE 분석 시 FDR 보정을 기본값으로 강제했다. 사용자가 명시적으로 raw p-value를 선택할 수는 있지만, 경고 메시지를 표시한다.

DIA-NN 공식 GitHub에서 DIA-NN의 최신 버전과 문서를 확인할 수 있다. 프로테오믹스 분석에서 통계적 엄밀성의 중요성은 아무리 강조해도 지나치지 않다.

💡 데이터의 정규성 검정과 적절한 통계 검정 방법 선택에 대해서는 다음 글: 데이터가 정규분포인지 어떻게 아나에서 다룬다.

프로테오믹스 분석 파이프라인의 전체 구조는 sbmlab.com에서 확인할 수 있다.