Seurat로 시작하는 단일세포 RNA-seq 분석 — pbmc3k로 QC부터 셀타입 주석까지 (R 실전)

TL;DR — DESeq2 같은 벌크 RNA-seq가 "조직 전체의 평균 발현"을 본다면, 단일세포 RNA-seq (scRNA-seq)는 세포 하나하나의 발현을 읽어 "이 조직에 어떤 세포가 몇 종류 있고 각자 뭘 하는지"를 밝힙니다. 그 표준 도구가 R 패키지 Seurat입니다.

이 글은 가장 유명한 입문 데이터 pbmc3k (말초혈 단핵세포 2,700개)로, 데이터 로딩 → 품질관리(QC) → 정규화 → 차원축소(PCA) → 클러스터링 → UMAP → 셀타입 주석까지 복붙하면 돌아가는 Seurat v5 코드로 따라갑니다.

그리고 v4 코드를 그대로 쓰면 깨지는 v5의 함정 3가지 — FindAllMarkers 기본값 변경, JoinLayers, 레이어 문법($counts) — 도 짚습니다. 마지막으로 가장 중요한 한 가지 — "클러스터는 알고리즘의 산물이지, 그 자체가 정답 셀타입이 아니다."

🔗 관련 글  ·  같은 R 분석 실전: DESeq2로 벌크 RNA-seq 차등발현  ·  clusterProfiler 농축분석

그림 1. scRNA-seq 분석 파이프라인 — 10x 행렬 → QC → 정규화 → PCA → 클러스터 → UMAP → 셀타입 주석.

0. 준비물 — 설치와 데이터

Seurat는 CRAN에서 바로 설치되고, 지금 설치하면 v5 (현재 5.4 계열)가 깔립니다.

install.packages("Seurat")   # v5 설치
library(Seurat)
library(dplyr)
library(patchwork)

데이터는 10x Genomics의 pbmc3k를 씁니다. pbmc3k_filtered_gene_bc_matrices.tar.gz를 받아 풀면 filtered_gene_bc_matrices/hg19/ 안에 세 파일(matrix.mtx, genes.tsv, barcodes.tsv)이 있습니다. 이게 세포 × 유전자 카운트 행렬입니다.

💡 SeuratData 패키지로 InstallData("pbmc3k")처럼 받을 수도 있지만, 그렇게 받는 객체는 이미 처리가 끝난 구버전 객체라 "원시 카운트부터 시작" 흐름과 어긋납니다. 입문에는 아래 Read10X 경로를 권합니다.

1. 데이터 로딩과 Seurat 객체 만들기

pbmc.data <- Read10X(data.dir = "filtered_gene_bc_matrices/hg19/")
pbmc <- CreateSeuratObject(counts = pbmc.data, project = "pbmc3k",
                           min.cells = 3, min.features = 200)
pbmc

min.cells = 3은 3개 미만 세포에서만 잡힌 유전자를, min.features = 200은 200개 미만 유전자만 검출된 세포(거의 빈 방울)를 처음부터 걸러 줍니다.

⚠️ v5 함정 ① — 카운트에 접근할 때 v4의 pbmc[["RNA"]]@counts (슬롯)는 v5에서 pbmc[["RNA"]]$counts (레이어)로 바뀌었습니다. v5 객체는 counts·data·scale.data를 '레이어'로 관리합니다.

2. 품질관리(QC) — 죽은 세포와 더블렛 거르기

단일세포 데이터에는 죽어가는 세포(미토콘드리아 유전자 비율↑)와 두 세포가 한 방울에 들어간 더블렛(검출 유전자 수↑)이 섞여 있습니다. 이걸 걸러야 합니다.

# 미토콘드리아 유전자 비율 (사람: "^MT-", 마우스: "^mt-")
pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")

# 분포 확인
VlnPlot(pbmc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")

# 필터링
pbmc <- subset(pbmc, subset = nFeature_RNA > 200 &
                              nFeature_RNA < 2500 &
                              percent.mt   < 5)

⚠️ 흔한 실수 — 사람은 "^MT-", 마우스는 "^mt-" (소문자)입니다. 패턴을 잘못 쓰면 percent.mt가 전부 0이 되어 죽은 세포를 못 거릅니다.

⚠️ 임계값은 상수가 아니다 — 200/2500/5%는 pbmc3k에 맞춘 값입니다. 조직·플랫폼마다 다르므로, 반드시 VlnPlot으로 분포를 먼저 보고 정하세요. 그리고 더블렛 제거는 Seurat 기본 기능이 아니라 DoubletFinder·scDblFinder 같은 별도 도구가 필요합니다(입문에선 생략).

3. 정규화 → 변이 유전자 → 스케일 → PCA

세포마다 읽힌 총량이 다르니 정규화하고, 세포를 구분하는 데 중요한 고변이 유전자(HVG)만 추려, 스케일한 뒤 PCA로 차원을 압축합니다.

pbmc <- NormalizeData(pbmc, normalization.method = "LogNormalize", scale.factor = 10000)
pbmc <- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 2000)

all.genes <- rownames(pbmc)
pbmc <- ScaleData(pbmc, features = all.genes)

pbmc <- RunPCA(pbmc, features = VariableFeatures(object = pbmc))
ElbowPlot(pbmc)   # 그래프가 꺾이는 지점 → pbmc3k는 약 10개 PC

ElbowPlot에서 곡선이 평평해지는 지점까지를 '의미 있는 차원'으로 봅니다. pbmc3k에서는 보통 앞쪽 10개 PC를 씁니다.

💡 대안 — SCTransform. 위 세 줄(정규화·HVG·스케일)을 한 함수로 대체할 수 있습니다.

pbmc <- SCTransform(pbmc, vars.to.regress = "percent.mt", verbose = FALSE)

기술적 노이즈 보정에 강해 실전 파이프라인에서 많이 씁니다. 단 이 경로에서는 보통 차원을 더 넉넉히(dims = 1:30) 잡습니다. 입문·해석 투명성에는 위의 LogNormalize가 직관적입니다.

4. 클러스터링과 UMAP

이제 비슷한 세포끼리 묶고(클러스터링), 2차원 지도(UMAP)에 펼칩니다.

pbmc <- FindNeighbors(pbmc, dims = 1:10)
pbmc <- FindClusters(pbmc, resolution = 0.5)   # 기본 Louvain 알고리즘

pbmc <- RunUMAP(pbmc, dims = 1:10)             # 기본은 uwot (R 네이티브)
DimPlot(pbmc, reduction = "umap", label = TRUE)

resolution = 0.5면 pbmc3k에서 약 9개 클러스터(0~8)가 나옵니다.

⚠️ resolution은 '클러스터 수 조절 노브'다 — 높이면 잘게 쪼개지고(과분할), 낮추면 뭉쳐집니다(과병합). 정답 resolution은 없습니다. 알려진 마커로 생물학적 타당성을 확인하며 정해야 합니다(clustree로 안정성 비교 가능). UMAP의 거리도 '대략의 이웃 관계'일 뿐, 축 자체에 절대적 의미는 없습니다.

그림 2. UMAP 클러스터링 — 비슷한 세포끼리 묶인 약 9개 클러스터(0~8). (예시 도식)

5. 마커 유전자와 셀타입 주석

클러스터마다 '특이적으로 켜진 유전자(마커)'를 찾고, 알려진 표지로 셀타입을 붙입니다.

# v5 함정 ②: 여러 샘플을 합치거나(merge) 통합한 객체는 레이어가 쪼개져 있어
# FindAllMarkers 전에 JoinLayers가 '필수'. 단일 pbmc3k엔 불필요하나 호출해도 무해.
pbmc[["RNA"]] <- JoinLayers(pbmc[["RNA"]])

# v5 함정 ③: FindAllMarkers 기본값이 v5에서 min.pct=0.01, logfc.threshold=0.1로
# '완화'됐다(v4는 0.1/0.25). 그대로 두면 마커가 폭증하므로, 고전처럼 엄격히 명시:
pbmc.markers <- FindAllMarkers(pbmc, only.pos = TRUE,
                               min.pct = 0.25, logfc.threshold = 0.25)

pbmc.markers %>% group_by(cluster) %>% slice_max(n = 5, order_by = avg_log2FC)

마커를 알려진 면역세포 표지와 맞추면 다음과 같이 정리됩니다.

클러스터	대표 마커	셀타입
0	IL7R, CCR7	Naive CD4 T
1	CD14, LYZ	CD14+ 단핵구
2	IL7R, S100A4	Memory CD4 T
3	MS4A1	B 세포
4	CD8A	CD8 T
5	FCGR3A, MS4A7	FCGR3A+ 단핵구
6	GNLY, NKG7	NK 세포
7	FCER1A, CST3	수지상세포(DC)
8	PPBP	혈소판

new.cluster.ids <- c("Naive CD4 T", "CD14+ Mono", "Memory CD4 T", "B", "CD8 T",
                     "FCGR3A+ Mono", "NK", "DC", "Platelet")
names(new.cluster.ids) <- levels(pbmc)
pbmc <- RenameIdents(pbmc, new.cluster.ids)

DimPlot(pbmc, reduction = "umap", label = TRUE, pt.size = 0.5) + NoLegend()
saveRDS(pbmc, file = "pbmc3k_final.rds")

마커 목록을 clusterProfiler로 농축분석에 넘기면, 각 셀타입이 어떤 경로를 쓰는지까지 볼 수 있습니다.

그림 3. 셀타입 주석 — 각 클러스터를 마커로 면역세포 종류에 매핑(Naive CD4 T·B·NK·단핵구 등).

6. v4 → v5, 알아둘 함정 정리

v4 시절 코드를 그대로 돌리면 막히는 대표 지점입니다.

• 레이어 문법 — @counts (슬롯) → $counts (레이어). 범용 접근자도 GetAssayData → LayerData 권장.
• JoinLayers — merge()·통합으로 쪼개진 레이어를 다시 합치는 v5 신설 함수. 멀티샘플 DE 전 필수. v4 코드엔 없어서 빠뜨리기 쉽습니다.
• FindAllMarkers 기본값 — logfc.threshold 0.25→0.1, min.pct 0.1→0.01로 완화. 모르고 쓰면 "마커가 너무 많이" 나옵니다.
• 통합 함수 — v4의 IntegrateData → v5는 IntegrateLayers() 한 함수(Harmony·CCA·RPCA 등 선택). 멀티샘플에서 배치 효과가 보이면 여기로.

7. 핵심 정리

• ✅ 단일세포 RNA-seq = 세포 하나하나의 발현을 읽어 군집 구조를 본다. 표준 도구는 Seurat (v5).
• ✅ 표준 흐름 — Read10X → QC (percent.mt·nFeature_RNA) → 정규화 → HVG → PCA → 이웃·클러스터 → UMAP → 마커 → 셀타입.
• ✅ v5 함정 — $counts 레이어 문법 · JoinLayers (멀티샘플) · FindAllMarkers 기본값 변경.
• ✅ QC 임계값·resolution은 데이터마다 다른 '선택'이지 정답 상수가 아니다.
• ⚠️ 가장 중요한 것 — 클러스터 ≠ 정답 셀타입. 클러스터는 알고리즘 산물이고, 셀타입 라벨은 마커 기반 사후 해석이다. resolution·마커 선택에 따라 바뀔 수 있으니, 자동 주석(SingleR·Azimuth)이나 알려진 생물학과 교차검증하라.

8. 다음 학습 추천

• 📊 DESeq2로 벌크 RNA-seq 차등발현 분석 — '조직 평균'을 보는 사촌 기법, 단일세포와 비교 (실전)
• 🧬 clusterProfiler로 농축분석 — 클러스터 마커 목록에 생물학적 의미 입히기 (실전)

References

1. Hao, Y., et al. (2024). Dictionary learning for integrative, multimodal, and scalable single-cell analysis (Seurat v5). Nature Biotechnology, 42, 293–304.
2. Satija Lab. Seurat — Guided Clustering Tutorial (pbmc3k). https://satijalab.org/seurat/articles/pbmc3k_tutorial
3. Satija Lab. Seurat v5 Command Cheat Sheet · SCTransform vignette · FindAllMarkers reference. https://satijalab.org/seurat/
4. 10x Genomics. 3k PBMCs from a Healthy Donor (filtered gene-barcode matrices).

Pipette & Pipeline · A bio portfolio journal — bridging research, data, and industry.