pseudotime은 순서를 줬지만 방향은 주지 못했습니다

S32에서 우리는 pseudotime을 통해 세포의 순서를 복원했습니다. 그런데 한 가지 치명적인 한계가 있었습니다 — 이 세포가 지금 앞으로 분화하고 있는지, 아니면 그 상태에 머물러 있는지 알 수 없다는 점입니다. pseudotime은 정적인 거리이지, 동적인 방향이 아닙니다.

2018년, La Manno 등이 발표한 RNA velocity는 이 문제에 놀랍도록 우아한 답을 줬습니다. 아이디어의 핵심은 이것입니다: 세포 안에는 아직 스플라이싱되지 않은 미성숙 mRNA(unspliced)와 완성된 성숙 mRNA(spliced)가 동시에 존재하며, 이 둘의 비율이 그 유전자가 지금 켜지고 있는지 꺼지고 있는지를 말해 줍니다. 유전자가 갓 켜졌다면 unspliced가 많고 spliced는 아직 적습니다. 반대로 유전자가 꺼지는 중이라면 새 전사가 멈추었으므로 unspliced는 줄고 spliced만 남아 분해됩니다.

이 편에서는 이 비율이 어떻게 세포 하나하나에 속도 벡터를 만들어 주는지 수학적으로 유도하고, scVelo로 직접 벡터장을 그려 봅니다.

스플라이싱 동역학: 미분 방정식 한 쌍

유전자 하나를 생각합시다. 전사 속도를 $\alpha$, 스플라이싱 속도를 $\beta$, 분해 속도를 $\gamma$라 하면, unspliced mRNA의 양 $u$와 spliced mRNA의 양 $s$는 다음 미분 방정식을 따릅니다.

정상 상태(steady state)에서는 $du/dt = 0$, $ds/dt = 0$이므로 $u_{\text{ss}} = \alpha / \beta$, $s_{\text{ss}} = \alpha / \gamma$가 됩니다. 이를 정리하면 정상 상태에서의 관계는 다음과 같습니다.

이것은 $(s, u)$ 평면에서 원점을 지나는 직선입니다. 기울기는 $\gamma / \beta$이며, 이 직선이 바로 "정상 상태 선"입니다.

정상 상태 선 위와 아래

세포가 이 직선 위에 정확히 놓여 있다면, 그 유전자는 전사와 분해가 평형을 이루고 있는 것입니다. 하지만 실제 세포는 직선에서 벗어나 있습니다.

• 직선 위쪽(unspliced가 예상보다 많음): 유전자가 유도 단계 — 전사가 갓 시작되어 unspliced가 쌓이는 중입니다. $ds/dt > 0$이므로 spliced가 앞으로 증가할 것입니다.
• 직선 아래쪽(unspliced가 예상보다 적음): 유전자가 억제 단계 — 전사가 꺼져 unspliced 공급이 끊겼고 spliced가 분해되는 중입니다. $ds/dt < 0$이므로 spliced가 줄어들 것입니다.

각 유전자의 "현재 위치에서 정상 상태까지의 잔차"가 그 유전자의 velocity입니다. 수천 개 유전자의 velocity를 종합하면, 세포 하나에 대해 유전자 공간에서의 속도 벡터가 만들어집니다. 이 벡터를 UMAP 좌표로 투영하면 — 화살표가 됩니다.

Velocyto와 scVelo: 두 세대

Velocyto (1세대, 2018)

La Manno et al.의 원래 방법입니다. BAM 파일에서 인트론 리드(unspliced)와 엑손 리드(spliced)를 분리하고, 유전자별로 $(s, u)$ 산점도에 직선을 피팅한 뒤 잔차를 velocity로 사용합니다. 단, 모든 세포가 동일한 $\gamma/\beta$를 공유한다고 가정하므로(정상 상태 가정), 과도기적 상태가 포착되지 않습니다.

scVelo (2세대, 2020)

Bergen et al.은 정상 상태 가정을 버리고, **동역학 모델(dynamical model)**을 도입했습니다. 각 유전자에 대해 전사 속도 $\alpha$가 세포 상태별로 변하도록 EM 알고리즘으로 $\alpha$, $\beta$, $\gamma$와 잠재 시간(latent time)을 동시에 추정합니다. 이 덕분에 유도/억제 위상 전환을 명시적으로 모델링할 수 있고, 결과적으로 더 정확한 방향 추정이 가능합니다.

실습: 췌장 내분비 분화 데이터로 화살표 그리기

scVelo 공식 튜토리얼의 췌장 데이터(Bastidas-Ponce et al. 2019)를 사용합니다. 전구세포(Ductal/Ngn3)에서 알파/베타/델타/엡실론 세포로 분화하는 과정이 담겨 있어 velocity의 방향성을 확인하기에 이상적입니다.

import scvelo as scv

# 데이터 로드 (내장 pancreas 데이터셋)
adata = scv.datasets.pancreas()

# 전처리: 필터 + 정규화 + 로그 + HVG
scv.pp.filter_and_normalize(adata, min_shared_counts=20, n_top_genes=2000)

# 모멘트 계산 (kNN 이웃 기반 1·2차 모멘트)
scv.pp.moments(adata, n_pcs=30, n_neighbors=30)

# 동역학 모델 피팅 (유전자별 alpha, beta, gamma + latent time)
scv.tl.recover_dynamics(adata)
scv.tl.velocity(adata, mode="dynamical")

# velocity 그래프 + UMAP 투영
scv.tl.velocity_graph(adata)
scv.pl.velocity_embedding_stream(adata, basis="umap", color="clusters")

velocity_embedding_stream이 UMAP 위에 물 흐르는 듯한 화살표 스트림을 그려 줍니다. Ductal 전구세포에서 출발해 Ngn3를 거쳐 각 내분비 세포 유형으로 갈라지는 흐름이 보이면 성공입니다.

개별 유전자의 위상 초상화

# 특정 유전자의 (s, u) 위상도 — 유도/억제 판별
scv.pl.velocity(adata, var_names=["Ins2", "Gcg"], color="clusters")

Ins2(인슐린)의 위상도에서는 베타 세포 클러스터가 유도 상태(직선 위쪽)에 몰려 있고, Gcg(글루카곤)에서는 알파 세포가 유도 상태에 있을 것입니다. 이 그림 하나가 "이 유전자가 이 세포군에서 지금 켜지고 있는가"에 대한 직관적 답입니다.

잠재 시간 vs pseudotime

scv.tl.latent_time(adata)
scv.pl.scatter(adata, color="latent_time", cmap="gnuplot")

scVelo의 latent time은 동역학 모델이 추정한 절대적 시간축이므로, S32의 diffusion pseudotime보다 생물학적으로 더 해석 가능한 경우가 많습니다. 다만 이것도 상대적 순서이지 실제 시간(시·일)은 아닙니다.

unspliced/spliced 카운트는 어디서 오는가

velocity 분석의 전제는 unspliced와 spliced를 분리할 수 있어야 한다는 것입니다. 두 가지 주요 경로가 있습니다.

• Velocyto CLI: BAM(10x CellRanger 출력)에서 인트론 영역에 매핑된 리드를 unspliced로 분류합니다. 결과는 .loom 파일로 나옵니다.
• kallisto bustools (kb-python): FASTQ에서 직접 스플라이싱 상태를 구분합니다. CellRanger 없이도 작동하며, 속도도 빠릅니다.

# kb-python으로 spliced/unspliced 동시 정량
kb count -i index.idx -g t2g.txt -x 10xv3 -o output \
   --workflow lamanno \
   R1.fastq.gz R2.fastq.gz

--workflow lamanno가 인트론 참조를 자동으로 생성하고, spliced/unspliced/ambiguous 세 행렬을 출력합니다. 이 행렬을 AnnData 레이어(adata.layers["spliced"], adata.layers["unspliced"])에 넣으면 scVelo 입력이 됩니다.

CS 매핑

• 미분(도함수): velocity의 핵심은 $ds/dt$, 즉 spliced mRNA의 시간 변화율입니다. 연속 함수의 도함수를 이산 데이터에서 추정하는 것이 수치 미분의 바이오 버전입니다.
• 벡터장(vector field): 각 세포에 속도 벡터가 할당되면 유전자 공간(또는 UMAP 2D)에 벡터장이 형성됩니다. 유체역학의 흐름도와 동일한 시각화입니다.
• ODE 안정점: 정상 상태 선은 ODE 시스템의 안정점(fixed point)입니다. 유도/억제는 안정점에서 벗어난 과도 응답으로, 제어 이론의 과도 응답 분석과 같은 구조입니다.

자주 만나는 결함

• unspliced 카운트 부족 — 3' 라이브러리(10x Chromium 3')는 인트론 커버리지가 낮아 unspliced가 희소합니다. 필터링이 너무 엄격하면 velocity를 추정할 유전자가 남지 않습니다. min_shared_counts를 낮춰 보거나, 5' 라이브러리를 고려합시다.
• 모든 유전자에 velocity를 기대 — velocity가 의미 있으려면 유도/억제 전환이 존재해야 합니다. 항상 켜져 있는 하우스키핑 유전자는 정상 상태 선 위에 놓이므로 velocity가 0에 가깝습니다. 이것은 결함이 아니라 정상입니다.
• 배치 효과를 velocity로 오해 — 배치 간 기술적 차이가 spliced/unspliced 비율에 영향을 줄 수 있습니다. S31의 배치 보정이 선행돼야 합니다.
• 화살표 방향만 보고 결론 — 스트림 플롯은 시각적으로 인상적이지만, velocity confidence(scv.tl.velocity_confidence)를 반드시 확인합시다. confidence가 낮은 영역의 화살표는 노이즈입니다.

더 깊게 파고 싶다면

본문은 BPD가 직접 재구성한 서술입니다. 심화는 아래 정통 자료를 활용합시다.

• 원 논문(RNA velocity): La Manno et al. (2018), RNA velocity of single cells, Nature 560:494. velocity 개념의 원천입니다.
• 원 논문(scVelo): Bergen et al. (2020), Generalizing RNA velocity to transient cell states through dynamical modeling, Nature Biotechnology 38:1408. 동역학 모델의 정수입니다.
• 참고 무료 웹북: Single-cell best practices (sc-best-practices.org) — RNA velocity 챕터가 방법 비교와 진단을 상세히 다룹니다.
• 후속 발전(UniTVelo): Gao et al. (2022), UniTVelo: temporally unified RNA velocity reinforces single-cell trajectory inference, Nature Communications 13:6586. 통합 시간 추정으로 정확도를 높였습니다.
• 후속 발전(veloVI): Gayoso et al. (2024), Deep generative modeling of transcriptional dynamics for RNA velocity analysis in single cells, Nature Methods 21:50. 변분 추론 기반 접근입니다.

velocity는 각 세포에 화살표를 줬습니다. 하지만 세포는 혼자 분화하지 않습니다 — 주변 세포와 리간드-수용체 신호를 주고받으며 운명을 결정합니다. 다음 편 S34에서는 이 세포간 통신 네트워크를 CellChat과 CellPhoneDB로 재구성합니다.