주성분분석(PCA)으로 RNA-seq 샘플 QC하기 — vst·plotPCA·prcomp·배치효과 (R 실전)

TL;DR — RNA-seq 분석을 시작하면 차등발현보다 먼저 해야 하는 게 샘플 품질 점검(QC)입니다. 그 1번 도구가 주성분분석(principal component analysis, PCA) — 수만 개 유전자를 PC1·PC2 두 축으로 압축해, 같은 그룹끼리 모이는지 한눈에 봅니다.

함정 두 가지만 기억하세요 — ① 원시 카운트(raw count)로 PCA 하면 안 됩니다. 먼저 vst()로 변환합니다(안 그러면 PC1이 '시퀀싱 깊이'가 됩니다). ② prcomp()는 행을 관측치로 봅니다 — 유전자×샘플 행렬은 t()로 전치해 샘플을 행으로 넣어야 합니다. (코드는 공개 데이터 airway 기준, DESeq2 1.46·factoextra 2.0.0 공식 문서로 작성했습니다.)

🔗 관련 글  ·  R 분석 워크플로: DESeq2로 차등발현 분석  ·  ComplexHeatmap 발현 히트맵  ·  Seurat 단일세포 분석

그림 1. PC1·PC2 산점도 — 처리(dex)로 색, 세포주로 모양. 축에 설명 분산 %를 표시.

0. 주성분분석(PCA)이 뭔가 — 1분 직관

유전자가 2만 개면, 샘플 하나는 2만 차원의 점입니다. 사람은 그걸 못 봅니다. PCA는 이 점들이 가장 많이 흩어진 방향을 찾아 새 축으로 삼습니다.

• PC1 = 데이터가 가장 크게 변하는 방향. PC2 = PC1과 직각이면서 그다음으로 크게 변하는 방향.
• 각 축은 '전체 분산의 몇 %를 설명하는지' 값을 가집니다. PC1이 80%면, 샘플 차이의 대부분이 한 방향으로 정리된다는 뜻입니다.

RNA-seq에서 PCA를 쓰는 이유는 단순합니다 — "같은 조건 샘플끼리 모이나? 튀는 샘플은 없나? 그룹이 갈리는 게 생물학 때문인가 배치(batch) 때문인가?"를 분석 시작 전에 확인하기 위해서입니다.

1. 준비물 — 설치

Bioconductor 패키지 DESeq2 (변환·PCA)와 예제 데이터 airway, 그리고 ggplot2를 씁니다.

if (!require("BiocManager", quietly = TRUE)) install.packages("BiocManager")
BiocManager::install(c("DESeq2", "airway"))   # 최초 1회
install.packages("ggplot2")

library(DESeq2)
library(airway)
library(ggplot2)

2. 예제 데이터 — airway 불러오기

airway는 기도 평활근 세포에 천식약(dexamethasone)을 처리한 RNA-seq로, 샘플 8개(세포주 4개 × 처리/비처리)입니다. 복붙으로 재현됩니다.

data("airway")
airway$dex <- relevel(airway$dex, ref = "untrt")   # 비처리(untrt)를 기준으로

dds <- DESeqDataSet(airway, design = ~ cell + dex)
dds <- dds[rowSums(counts(dds)) > 1, ]             # 전부 0인 유전자 제거(가벼운 사전 필터)
dds

class: DESeqDataSet
dim: 63677 8 ...
colData names(9): ... cell dex ...

colData(dds)에 dex (처리 여부)와 cell (세포주)이 들어 있습니다. 이게 나중에 점의 색·모양이 됩니다.

3. 첫 함정 — 원시 카운트로 PCA 하지 마라

RNA-seq 카운트는 평균이 크면 분산도 커지는 성질(음이항 분포)이 있습니다. 이대로 PCA를 돌리면 발현량 큰 유전자 몇 개와 시퀀싱 깊이가 PC1을 지배해, 생물학이 아니라 '총 read 수'를 보게 됩니다. 그래서 먼저 분산안정화변환(variance stabilizing transformation, vst)으로 분산을 평탄화합니다.

vsd <- vst(dds, blind = FALSE)   # 분산안정화 변환 (로그 스케일 비슷)
# 샘플이 적거나(n<30) 깊이 편차가 크면 rlog(dds, blind=FALSE)가 더 안전

그림 2. 원시 카운트는 PC1이 '시퀀싱 깊이'에 끌려가지만, vst 변환 후엔 처리 효과가 드러난다.

• vst() — 빠르고 대규모(수백 샘플)에 적합. rlog() — 샘플이 적거나 라이브러리 크기 편차가 클 때 더 안정적(대신 느림).
• blind=FALSE는 변환의 분산 추정에만 설계를 참고할 뿐, 여전히 비지도(unsupervised) 변환입니다.

⚠️ 핵심 — PCA·히트맵·샘플 거리 같은 시각화/QC에는 항상 vst (또는 rlog) 값을 씁니다. 차등발현 분석 자체는 원시 카운트로 DESeq()에 넣습니다(둘은 입력이 다릅니다).

4. 가장 빠른 길 — plotPCA() 한 줄

DESeq2에는 PCA 함수가 내장돼 있습니다. 분산이 큰 상위 500개 유전자(ntop=500)로 자동 계산해 ggplot 객체를 돌려줍니다.

plotPCA(vsd, intgroup = "dex")               # 처리 여부로 색
# plotPCA(vsd, intgroup = c("dex", "cell"))   # 두 변수 함께
# plotPCA(vsd, intgroup = "dex", pcsToUse = 3:4)  # PC3 vs PC4도 가능

한 줄이면 끝나지만, 축 라벨이 밋밋하고 색·모양을 더 손보고 싶을 때가 많습니다. 그래서 데이터만 받아 직접 그립니다.

5. 축에 '몇 % 분산'을 — returnData + ggplot2

returnData = TRUE로 좌표를 받고, percentVar 속성으로 축 라벨을 만듭니다.

pcaData <- plotPCA(vsd, intgroup = c("dex", "cell"), returnData = TRUE)
percentVar <- round(100 * attr(pcaData, "percentVar"))   # 각 축의 설명 분산 %

ggplot(pcaData, aes(PC1, PC2, color = dex, shape = cell)) +
  geom_point(size = 4) +
  xlab(paste0("PC1: ", percentVar[1], "% variance")) +
  ylab(paste0("PC2: ", percentVar[2], "% variance")) +
  coord_fixed() +            # 두 축의 스케일을 같게(거리 왜곡 방지)
  theme_bw()

• 색 = dex (처리), 모양 = cell (세포주)로 두 변수를 한 그림에 담았습니다.
• coord_fixed()가 중요합니다 — PC1·PC2의 눈금 비율을 같게 해야 거리가 왜곡되지 않습니다.

6. scree plot — 축을 몇 개나 볼까

각 PC가 분산의 몇 %를 설명하는지 막대로 보는 게 scree plot입니다. prcomp()로 직접 계산해 봅니다(plotPCA와 같은 상위 500개 유전자 기준).

# plotPCA 내부와 동일하게: 분산 상위 500개 유전자 선택
rv     <- matrixStats::rowVars(assay(vsd))     # rowVars는 base R이 아니다 — 네임스페이스 주의!
select <- order(rv, decreasing = TRUE)[seq_len(min(500, length(rv)))]

pr <- prcomp(t(assay(vsd)[select, ]))   # t()로 '샘플을 행'으로! vst 값엔 scale. 불필요
pv <- pr$sdev^2 / sum(pr$sdev^2)        # 설명 분산 비율

barplot(100 * pv[seq_len(min(8, length(pv)))],
        names.arg = paste0("PC", seq_len(min(8, length(pv)))),
        ylab = "% variance explained", main = "Scree plot (airway)")

# pr$x      : 점수(scores) — 샘플 × PC  → 산점도에 찍는 값
# pr$rotation: 적재(loadings) — 유전자 × PC → 어떤 유전자가 그 축을 만드는가
# pr$sdev    : 각 PC의 표준편차

• 샘플이 8개면 0이 아닌 PC는 최대 7개입니다(그래서 min(8, …)).
• scale. 인자: vst 값엔 안 씁니다(FALSE). vst가 이미 분산을 맞춰 놨기 때문입니다. 원시·로그 값에 PCA할 때만 scale.=TRUE로 표준화합니다.

7. 해석 — PC1이 '생물학'이냐 '배치'냐

그림 3. PC1이 조건으로 갈리면 호재, 배치로 갈리면 배치효과. scree로 '몇 개 축이 의미 있나'를 본다.

PCA 그림은 이렇게 읽습니다.

• 그룹끼리 모이고, 조건이 PC1에서 갈린다 → 처리 효과가 가장 큰 변동. 좋은 신호.
• 그룹과 무관하게, 배치·레인·실험 날짜로 갈린다 → 배치효과(batch effect). PC1이 생물학이 아니라 기술적 잡음입니다. 설계에 배치를 넣거나(~ batch + condition), 시각화용으로 limma::removeBatchEffect()를 씁니다.
• 한 샘플만 멀리 떨어져 있다 → 이상치(outlier) 의심. 라이브러리 크기·매핑률·RNA 분해를 함께 확인합니다.
• PC1이 80~90% → 한 축이 압도(강한 처리 효과거나 이상치). scree가 평평 → 뚜렷한 구조 없음(잡음 위주).

💡 실전 팁 — PCA를 관심 변수(처리)로도, 의심 변수(배치)로도 각각 색칠해 두 장 그려 보세요. 배치로 칠했을 때 깔끔하게 갈리면 그게 바로 배치효과입니다.

8. 자주 나는 에러 & 해결

• could not find function "rowVars" — base R이 아닙니다. matrixStats::rowVars() 또는 MatrixGenerics::rowVars()로 호출하세요. (복붙했을 때 가장 흔히 깨지는 지점)
• 점이 수천 개로 흩어진다 — t()를 빠뜨려 유전자로 PCA한 것입니다. prcomp(t(mat))로 샘플을 행에 두세요.
• PC1이 그냥 시퀀싱 깊이 같다 — 원시 카운트로 돌렸을 때입니다. vst()/rlog() 먼저.
• 그룹이 안 갈린다 — 정상일 수 있습니다(효과가 약하거나, PC1·PC2 밖에 있음). pcsToUse로 PC3·PC4도 보세요.

9. 확장 — PCA vs UMAP·t-SNE

PCA는 선형·결정론적이라 축에 '몇 % 분산'이라는 해석이 붙고, 벌크 QC와 전체 구조 파악에 이상적입니다. 반면 UMAP·t-SNE는 비선형·확률적이라 미세한 국소 구조를 보기 좋지만, 점 사이 거리·군집 크기를 정량적으로 해석하면 안 됩니다.

그래서 단일세포 RNA-seq에서는 보통 PCA를 먼저 돌려 30~50개 PC로 잡음을 줄인 뒤, 그 결과를 UMAP에 넣습니다. 둘은 경쟁이 아니라 순서대로 쓰는 짝입니다.

정리 & 다음 학습

• ✅ RNA-seq는 차등발현 전에 PCA로 샘플 QC — 그룹이 모이나, 이상치는 없나.
• ✅ 원시 카운트 금지 → vst()/rlog() 변환 후 PCA.
• ✅ prcomp(t(...)) — 전치로 샘플을 행에. vst 값엔 scale.=FALSE.
• ✅ 가장 빠른 길은 plotPCA(), 다듬으려면 returnData=TRUE + ggplot2.
• ✅ PC1이 조건이면 호재, 배치면 배치효과 — 두 변수로 각각 칠해 확인.

🎓 다음 학습
• (심화) DESeq2로 차등발현(DEG) 분석 — PCA로 QC 끝났으면 다음은 DEG
• (비교) ComplexHeatmap 발현 히트맵 — 샘플 거리·발현 패턴 시각화
• (참고) Seurat 단일세포 분석 — PCA→UMAP 차원축소

References

1. Love, M. I., Huber, W., Anders, S. (2014). Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biology, 15, 550.
2. DESeq2 공식 vignette — Data transformations and visualization (Bioconductor 3.20, DESeq2 1.46). plotPCA(ntop=500, returnData, pcsToUse).
3. Himes, B. E., et al. (2014). RNA-seq transcriptome profiling identifies CRISPLD2 ... (airway 데이터셋). PLoS ONE, 9, e99625.
4. Kassambara, A., Mundt, F. factoextra: PCA 시각화(fviz_eig·fviz_pca_ind·fviz_pca_biplot), v2.0.0 (2026).
5. McInnes, L., et al. (2018). UMAP. arXiv:1802.03426. / van der Maaten & Hinton (2008). t-SNE. JMLR.

Pipette & Pipeline · A bio portfolio journal