CITE-seq integration 어떤 embedding method 쓸것인가 to clustering results

Posted Aug 31, 2026

By Subin Cho 17 min read

scRNA + ADT + TCR/BCR 멀티오믹 데이터에서 WNN·totalVI·GEX-only 중 무엇을 고를지, 6단계 진단으로 결정한 기록.

1. 문제 상황

데이터

항목	내용
샘플	59개 (51개는 GEX+ADT, 8개는 GEX only)
Assay	scRNA-seq + CITE-seq(137-plex TotalSeq-C, isotype 7개 포함) + TCR/BCR (TCR/BCR 결과는 integration계산에서 제외)
셀 수	전체 약 359,000 / ADT 보유 51샘플 약 342,000
결측	8개 샘플은 GEX+TCR/BCR만 (ADT 없음). 한 환자데이터에 집중
기술적 변수	중간에 sequencing core 변경 (chemistry 동일, 처리 기관 상이)
Batch 보정	모든 방법에서 SampleID 기준

시도한 세 가지

WNN + Harmony — GEX·ADT를 각각 SampleID로 Harmony 통합 후 WNN으로 결합. 51샘플 object를 reference로 두고 ADT 없는 8샘플을 query로 projection.
totalVI — 1번과 동일한 순서, python/scvi-tools 기반.
GEX-only Harmony — ADT를 아예 쓰지 않음.

증상

WNN UMAP이 cell type을 가장 잘 반영하는 것처럼 보였다 (Azimuth 라벨 기준 육안 판단).
그런데 같은 그래프에서 나온 clustering이 cell type을 따라가지 않았다.
클러스터들이 UMAP 상에서 낱알처럼 흩어져 있었다 — 공간적으로 뭉치지 않음.
클러스터링 알고리즘을 MSC로 바꿔봐도 동일.

R1.WNN UMAP

R2.totalVI UMAP

R3.GEX only UMAP

시작하기 전에 잡고 간 두 가지 전제

전제 1 — 알고리즘을 바꿔도 입력 그래프가 같으면 같은 구조가 나온다. Leiden → Louvain → MSC를 돌려가며 비교하는 건 의미가 없다. 고칠 지점은 그래프 위쪽이다.

전제 2 — UMAP 외형은 증거가 아니다. 2D 투영이 예쁜 것과 embedding이 cell type을 잘 보존하는 것은 다른 문제다. “어느 그림이 나아 보이나”를 “어느 embedding의 NMI/ARI가 높나”로 바꿔야 한다.

또한 clustering과 visualization은 같은 공간에서 할 의무가 없다. 그래프는 A로, 그림은 B로 가는 조합이 정답일 수 있다.

2. 확인 항목

① 클러스터가 대체 무엇을 잡고 있는가

Cell type이 아니라면 무엇인지부터 특정한다. 클러스터 라벨과 각 공변량 사이의 normalized mutual information을 재서 가장 높은 것이 클러스터의 정체다.

        
      
library(aricode)
md <- NYU@meta.data
cl <- as.factor(md$harmony_wnn_clusters)

md$mt_bin  <- cut(md$percent.mt,  quantile(md$percent.mt,  0:10/10, na.rm = TRUE), include.lowest = TRUE)
md$cnt_bin <- cut(md$nCount_SCT,  quantile(md$nCount_SCT,  0:10/10, na.rm = TRUE), include.lowest = TRUE)

vars <- c("predicted.celltype.l2", "Seqcore", "SampleID", "PatientID",
          "Timepoint", "mt_bin", "cnt_bin")
sapply(vars, function(v) NMI(cl, as.factor(md[[v]])))

결과

변수	NMI
predicted.celltype.l2	0.233
SampleID	0.068
PatientID	0.058
Timepoint	0.017
Seqcore	0.002

해석

Seqcore 0.002 — sequencing core 변경 효과는 SampleID Harmony로 완전히 흡수됐다. 가장 걱정했던 변수인데 흔적이 없다.
SampleID 0.068, PatientID 0.058 — batch도 donor도 아니다. Harmony가 제 일을 했다. (순서가 SampleID > PatientID > Timepoint인 것은 nesting 구조상 자연스럽다.)
celltype 0.233 — cell type도 아니다. 잘 된 PBMC clustering이라면 l2 기준 0.55~0.75가 나온다.

→ cell type도 batch도 아닌 무언가를 잡고 있다. 첫 번째 가설(batch 오염)은 기각.

⚠️ global NMI는 “특정 클러스터 몇 개만 한 샘플 전용”인 경우를 놓친다. 같이 볼 것:
1 2 smax <- apply(prop.table(table(cl, md$SampleID), 1), 1, max) sort(smax, decreasing = TRUE)[1:10] # 0.5 넘으면 sample-specific 클러스터

② 정답 라벨 자체가 믿을 만한가

①의 분모가 되는 라벨이 뭉개져 있으면 진단 전체가 무효다. 반드시 먼저 걸러야 한다.

        
summary(NYU$predicted.celltype.l2.score)
table(NYU$predicted.celltype.l2 == "Doublet")

결과

   Min. 1st Qu.  Median    Mean  3rd Qu.    Max.
 0.1851  0.6849  0.8861  0.8207  0.9898  1.0000

 FALSE   TRUE
341835     36

해석

median 0.886, 3rd Qu 0.99, Doublet 36개/342,000. 라벨은 깨끗하다. → 두 번째 가설(라벨 품질) 기각. ①의 0.233은 진짜 신호다.

이후 모든 계산은 고신뢰 셀로 제한한다:

        
keep <- md$predicted.celltype.l2.score > 0.75 & md$predicted.celltype.l2 != "Doublet"
ct   <- droplevels(as.factor(md$predicted.celltype.l2[keep]))

③ 클러스터 크기 분포 — singleton 함정

purity나 crosstab을 보기 전에 반드시 확인한다.

        
      
sort(table(md$harmony_wnn_clusters), decreasing = TRUE)

tab1 <- prop.table(table(cl[keep], md$predicted.celltype.l1[keep]), 1)
tab2 <- prop.table(table(cl[keep], ct), 1)
data.frame(
  n      = as.vector(table(cl[keep])),          # ← n 없이는 읽을 수 없다
  pur_l1 = round(apply(tab1, 1, max), 2),
  top_l1 = colnames(tab1)[apply(tab1, 1, which.max)],
  pur_l2 = round(apply(tab2, 1, max), 2),
  top_l2 = colnames(tab2)[apply(tab2, 1, which.max)]
)[order(-as.vector(table(cl[keep]))), ]

결과 — 25개 클러스터 중

실질 클러스터 8개 (purity 0.26 ~ 0.50)
셀 1~2개짜리 16개 (purity 정확히 1.00)

crosstab에서 비율이 0.00 / 0.50 / 1.00만 나오면 그 클러스터는 셀이 1~2개다.

해석

FindClusters(algorithm = 4, resolution = 0.02) — Leiden이 거대 덩어리 몇 개와 singleton 무더기를 뱉은 것. UMAP에서 라벨 9~24가 중앙 한 점에 몰려 있던 이유이기도 하다 (라벨은 median 좌표에 찍히므로, 흩어진 클러스터나 미세 클러스터는 중앙으로 모인다).

purity를 l1과 l2 둘 다 봐야 하는 이유 27개 l2 라벨 기준 purity는 정상 클러스터도 깎는다. 순수한 B 섬이라도 B naive 45% / B memory 40% / B intermediate 15%로 갈리므로 max가 0.45다. l1 purity가 0.85+이면 lineage는 잡힌 것이고, 필요한 건 그래프 재구축이 아니라 lineage별 2단계 subclustering이다.

④ Modality weight 분포 (WNN 전용)

WNN은 셀마다 두 modality의 가중치를 계산한다. 이 분포가 건강한지가 WNN 품질의 핵심 지표다.

        
      
summary(NYU_wADT$ADT.weight)

NYU_wADT$w_bin <- cut(NYU_wADT$ADT.weight, c(-Inf, .05, .95, Inf),
                      labels = c("RNA-only", "mixed", "ADT-only"))
table(NYU_wADT$w_bin)
FeaturePlot(NYU_wADT, "ADT.weight", reduction = "harmony.wnn.umap")

결과

     Min.   1st Qu.    Median      Mean    3rd Qu.      Max.
0.000e+00  1.51e-06  2.731e-01  4.462e-01  9.963e-01  1.000e+00

RNA-only    mixed  ADT-only
  145019    89357    107495      (42% / 26% / 31%)

해석

정상이라면 0.3~0.7에 몰린 단봉 분포여야 한다. 실제로는 1사분위가 0, 3사분위가 1인 U자형. 셀의 73%가 한쪽 modality만 보고 이웃을 정했다.

→ WNN 그래프가 두 개의 서로 다른 그래프를 셀 단위로 기워 붙인 모자이크가 됐다. 같은 cell type의 두 세포가 어느 modality가 이겼느냐에 따라 다른 동네에 배치된다.

원인 후보를 좁힌다:

        
cor(NYU_wADT$ADT.weight, log1p(NYU_wADT$nCount_ADT), method = "spearman")   # 0.20
round(prop.table(table(NYU_wADT$w_bin, NYU_wADT$SampleID), 2), 2)           # 샘플 간 균일

spearman 0.20 → ADT depth 탓이 아니다 (세 번째 가설 기각)
샘플별 비율이 0.37~0.55로 균일 → staining batch 효과도 아니다 (네 번째 가설 기각)

배치 아티팩트가 아니라 이 데이터에서 WNN 계산이 체계적으로 만들어내는 결과다.

⑤ Weight가 cell type과 정렬되는가

        
      
round(prop.table(table(NYU$w_bin, NYU$predicted.celltype.l1), 2), 2)

결과

	B	CD4 T	CD8 T	DC	Mono	NK	other T
RNA-only	0.48	0.47	0.34	0.71	0.55	0.36	0.29
mixed	0.34	0.23	0.24	0.16	0.35	0.26	0.25
ADT-only	0.19	0.30	0.42	0.12	0.10	0.38	0.47

해석

Myeloid/DC는 RNA 쪽, T/NK는 ADT 쪽으로 기운다. 생물학적으로 타당하다 — 137-plex 패널은 CD4/CD8/CD45RA/CCR7로 T subset을 잘 가르지만 DC·mono subset은 transcriptome이 낫다. WNN이 설계 의도대로 작동하고는 있다.

문제는 어느 타입도 한쪽으로 쏠리지 않는다는 점이다. CD8 T 안에서 34%는 RNA로, 42%는 ADT로 이웃이 정해진다. 같은 타입인데 서로 다른 공간에서 이웃을 찾는다.

여기까지로는 “치명적”인지 “품질을 깎는 정도”인지 단정할 수 없다. 다음 단계가 필요하다.

⑥ 3-way 벤치마크 — 이것이 실제 답이다

같은 라벨, 같은 필터, 비슷한 클러스터 수로 세 방법을 나란히 잰다.

        
      
score <- function(cl) {
  c1 <- droplevels(as.factor(cl[keep]))
  c(k   = nlevels(c1),
    NMI = round(aricode::NMI(c1, ct), 3),
    ARI = round(aricode::ARI(c1, ct), 3))
}
rbind(WNN     = score(md$harmony_wnn_clusters),
      GEXonly = score(md$gex_harmony_clusters),
      totalVI = score(md$totalvi_clusters))

세 결과가 다른 object에 흩어져 있으면 barcode 기준으로 먼저 합친다:

        
      
md$cell   <- rownames(md)
md$gex_cl <- gex@meta.data[md$cell, "<컬럼명>"]
md$tvi_cl <- tvi@meta.data[md$cell, "<컬럼명>"]
keep <- keep & !is.na(md$gex_cl) & !is.na(md$tvi_cl)   # 세 방법 모두 있는 셀만

결과

방법	k	NMI	ARI
totalVI	16	0.688	0.493
GEX-only	15	0.528	0.386
WNN	20	0.269	0.143

해석 — 세 가지를 동시에 말한다

WNN은 클러스터가 가장 많은데 점수가 가장 낮다. k=20으로 GEX-only(k=15)보다 많은데 NMI는 절반이다. → 클러스터 수 아티팩트가 아니라 진짜로 나쁘다. ③에서 유보했던 판단이 여기서 확정된다.
WNN < GEX-only. ADT를 WNN 방식으로 넣으면 아예 안 넣느니만 못하다. 정보를 더한 게 아니라 파괴했다.
totalVI > GEX-only (0.688 vs 0.528). 그런데 ADT에는 진짜 정보가 들어 있다. 제대로 모델링하면 이득이 나오는데 WNN이 그걸 못 쓴 것이다.

이 세 줄이 ④의 U자형 분포와 정확히 맞물린다. WNN은 CLR 값을 액면가로 받아 셀마다 modality를 양자택일하고, totalVI는 protein background를 negative binomial mixture로 명시적으로 모델링해 신호와 배경을 분리한다. 진단과 결과가 일치한다.

이 비교는 오히려 GEX-only에 유리하게 편향돼 있다. 정답 라벨이 SCT+spca 기반 GEX-only 매핑이기 때문이다. 그 핸디캡을 안고도 totalVI가 이겼다.

3. 결론

무엇을 선택했나

totalVI를 visualization·clustering 양쪽의 기본 embedding으로 채택.

판단 근거	값
totalVI가 1위	NMI 0.688 / ARI 0.493 (k=16)
WNN이 GEX-only보다도 낮음	0.269 vs 0.528 — ADT를 WNN으로 넣으면 역효과
원인	ADT.weight U자형 포화 (Q1=0, Q3=0.996), 셀의 73%가 단일 modality
원인이 batch가 아님	Seqcore NMI 0.002, SampleID 0.068
라벨 신뢰 가능	prediction score median 0.886, Doublet 36개

WNN 계열 파이프라인은 정리한다. DSB normalization으로 ADT를 고쳐 WNN을 살리는 길은 남아 있지만, 이미 이긴 방법이 있는 상황에서 들일 노력은 아니다.

후속 작업

resolution sweep — k=16은 27개 l2 타입에 비해 부족하다. 클러스터 수가 라벨 수보다 적으면 NMI 상한이 눌린다. k≈30~35에서 NMI가 평평해지는 지점을 찾는다.
59샘플 joint totalVI — reference→query projection을 버린다. totalVI는 protein이 전혀 없는 batch를 mask로 native하게 처리하므로 한 번에 학습시키면 projection 아티팩트가 사라진다.
ADT gating으로 독립 검증 — Azimuth 라벨은 GEX 기반이라 순환 논리 위험이 있다. CD3/CD4/CD8a/CD19/CD14/CD16/CD56 gating 라벨로 재확인.
lineage별 2단계 subclustering — 27개 l2 타입을 global resolution 하나로 쪼개려는 것은 무리다. lineage를 먼저 나누고 각 lineage 안에서 재-embedding.

4. 재사용 체크리스트

새 멀티오믹 데이터에서 embedding을 고를 때 이 순서로 간다.

□ 0. 눈으로 고르지 않는다. UMAP 외형은 증거가 아니다.
□ 1. 정답 라벨의 신뢰도부터 확인 (prediction score, Doublet 비율)
      → 여기가 무너지면 이하 전부 무효
□ 2. 클러스터 크기 분포 확인 (singleton 걸러내기)
      → n 없는 purity 표는 읽을 수 없다
□ 3. 클러스터 vs 공변량 NMI (celltype / batch / donor / QC)
      → 클러스터의 정체를 특정
□ 4. (WNN인 경우) modality weight 분포가 단봉인지
      → U자형이면 그래프가 모자이크
□ 5. 후보 embedding들을 같은 라벨·비슷한 k로 NMI/ARI 비교
      → 이것이 실제 의사결정 근거
□ 6. 독립 라벨(ADT gating 등)로 교차 검증
      → 순환 논리 방지

클러스터 수는 반드시 맞춰서 비교한다. k가 라벨 수보다 크게 적으면 NMI가 구조적으로 깎인다.

6. 메모

이 과정에서 세운 가설 네 개 중 셋이 기각됐다.

가설	결과	기각 근거
클러스터가 batch/donor를 잡고 있다	✗	Seqcore NMI 0.002, SampleID 0.068
Azimuth 라벨이 뭉개져 있다	✗	median score 0.886
ADT depth가 weight를 결정한다	✗	spearman 0.20
WNN의 modality weighting이 포화됐다	✓	Q1=0, Q3=0.996 + 3-way 벤치마크

기각된 가설도 각각 하나씩 가능성을 지웠기 때문에 의미가 있다. 다만 ⑥의 3-way 벤치마크를 먼저 돌렸다면 훨씬 빨리 도착했을 것이다. “왜 안 되는가”를 파기 전에 “무엇이 이기는가”를 먼저 재는 편이 낫다. 원인 규명은 이긴 방법이 없을 때 하는 일이다.

Study, Single-Cell

seurat tool

This post is licensed under CC BY 4.0 by the author.