scvelo

scVelo — RNA 속도 분석

개요

scVelo는 단일 세포 RNA-seq 데이터에서 RNA 속도 분석을 위한 선도적인 Python 패키지입니다. mRNA 스플라이싱의 동역학을 모델링하여 세포 상태 전이를 추론합니다. 이를 위해 unspliced(전구 mRNA)와 spliced(성숙 mRNA) 존재비의 비율을 사용하여 각 세포에서 유전자가 상향 또는 하향 조절되는지 여부를 결정합니다. 이를 통해 시간 경과 데이터 없이도 발생 궤적을 재구성하고 세포 운명 결정을 식별할 수 있습니다.

설치: uv pip install scvelo

주요 자료:

• 문서: https://scvelo.readthedocs.io/
• GitHub: https://github.com/theislab/scvelo
• 논문: Bergen et al. (2020) Nature Biotechnology. PMID: 32747759

이 스킬을 사용해야 하는 경우

다음과 같은 경우 scVelo를 사용하세요:

• 스냅샷 데이터로부터 궤적 추론: 세포가 어떤 방향으로 분화하는지 결정
• 세포 운명 예측: 전구 세포와 그 하류 운명 식별
• 드라이버 유전자 식별: 관찰된 궤적을 가장 잘 설명하는 동역학을 가진 유전자 발견
• 발달 생물학: 조혈, 신경발생, 상피-간엽 전이 모델링
• 잠재 시간 추정: 스플라이싱 동역학에서 파생된 가성시간을 따라 세포 정렬
• Scanpy 보완: UMAP 임베딩에 방향 정보 추가

전제 조건

scVelo는 unspliced 및 spliced RNA의 카운트 매트릭스가 필요합니다. 이는 다음을 통해 생성됩니다:

1. STARsolo 또는 kallisto|bustools (lamanno 모드)
2. velocyto CLI: velocyto run10x / velocyto run
3. alevin-fry / simpleaf (spliced/unspliced 출력)

데이터는 layers["spliced"] 및 layers["unspliced"]가 있는 AnnData 객체에 저장됩니다.

표준 RNA 속도 워크플로우

1. 설정 및 데이터 로드

import scvelo as scv
import scanpy as sc
import numpy as np
import matplotlib.pyplot as plt

# 설정 구성
scv.settings.verbosity = 3       # 계산 단계 표시
scv.settings.presenter_view = True
scv.settings.set_figure_params('scvelo')

# 데이터 로드 (spliced/unspliced 레이어가 있는 AnnData)
# 옵션 A: loom 파일에서 로드 (velocyto 출력)
adata = scv.read("cellranger_output.loom", cache=True)

# 옵션 B: velocyto loom을 Scanpy 처리된 AnnData와 병합
adata_processed = sc.read_h5ad("processed.h5ad")  # UMAP, 클러스터 포함
adata_velocity = scv.read("velocyto.loom")
adata = scv.utils.merge(adata_processed, adata_velocity)

# 레이어 확인
print(adata)
# obs × var: N × G
# layers: 'spliced', 'unspliced' (필수)
# obsm['X_umap'] (시각화에 필수)

2. 전처리

# 필터링 및 정규화. scVelo 0.3부터 filter_and_normalize()는 유전자 필터링 및
# 세포별 정규화만 수행합니다. n_top_genes를 더 이상 사용하지 않으며
# 로그 변환을 수행하지 않으므로 로그 단계 및 HVG 선택은 Scanpy에서 수행합니다.
scv.pp.filter_and_normalize(
    adata,
    min_shared_counts=20    # spliced+unspliced 최소 카운트
)
sc.pp.log1p(adata)
sc.pp.highly_variable_genes(adata, n_top_genes=2000, subset=True)

# 1차 및 2차 모멘트 계산 (평균 및 분산)
# knn_connectivities를 먼저 계산해야 합니다.
sc.pp.neighbors(adata, n_neighbors=30, n_pcs=30)
scv.pp.moments(
    adata,
    n_pcs=30,
    n_neighbors=30
)

3. 속도 추정 — 확률론적 모델

확률론적 모델은 빠르며 탐색적 분석에 적합합니다:

# 확률론적 속도 (빠르지만 정확도는 낮음)
scv.tl.velocity(adata, mode='stochastic')
scv.tl.velocity_graph(adata)

# 시각화
scv.pl.velocity_embedding_stream(
    adata,
    basis='umap',
    color='leiden',
    title="RNA 속도 (확률론적)"
)

4. 속도 추정 — 동역학 모델 (권장)

동역학 모델은 전체 스플라이싱 동역학을 피팅하며 더 정확합니다:

# 동역학 복원 (계산 집약적; 10K 세포 기준 약 10-30분 소요)
scv.tl.recover_dynamics(adata, n_jobs=4)

# 동역학 모델에서 속도 계산
scv.tl.velocity(adata, mode='dynamical')
scv.tl.velocity_graph(adata)

5. 잠재 시간

동역학 모델을 사용하면 공유 잠재 시간(pseudotime)을 계산할 수 있습니다:

# 잠재 시간 계산
scv.tl.latent_time(adata)

# UMAP에 잠재 시간 시각화
scv.pl.scatter(
    adata,
    color='latent_time',
    color_map='gnuplot',
    size=80,
    title='잠재 시간'
)

# 잠재 시간으로 정렬된 상위 유전자 식별
top_genes = adata.var['fit_likelihood'].sort_values(ascending=False).index[:300]
scv.pl.heatmap(
    adata,
    var_names=top_genes,
    sortby='latent_time',
    col_color='leiden',
    n_convolve=100
)

6. 드라이버 유전자 분석

# 가장 높은 속도 피팅을 가진 유전자 식별
scv.tl.rank_velocity_genes(adata, groupby='leiden', min_corr=0.3)
df = scv.DataFrame(adata.uns['rank_velocity_genes']['names'])
print(df.head(10))

# 속도 및 일관성
scv.tl.velocity_confidence(adata)
scv.pl.scatter(
    adata,
    c=['velocity_length', 'velocity_confidence'],
    cmap='coolwarm',
    perc=[5, 95]
)

# 특정 유전자에 대한 위상 초상화
scv.pl.velocity(adata, ['Cpe', 'Gnao1', 'Ins2'],
               ncols=3, figsize=(16, 4))

7. 속도 화살표 및 가성시간

# UMAP에 화살표 그림
scv.pl.velocity_embedding(
    adata,
    arrow_length=3,
    arrow_size=2,
    color='leiden',
    basis='umap'
)

# 스트림 그림 (더 깔끔한 시각화)
scv.pl.velocity_embedding_stream(
    adata,
    basis='umap',
    color='leiden',
    smooth=0.8,
    min_mass=4
)

# 속도 가성시간 (잠재 시간의 대안)
scv.tl.velocity_pseudotime(adata)
scv.pl.scatter(adata, color='velocity_pseudotime', cmap='gnuplot')

8. PAGA 궤적 그래프

# 속도 정보를 활용한 PAGA 그래프
scv.tl.paga(adata, groups='leiden')
df = scv.get_df(adata, 'paga/transitions_confidence', precision=2).T
df.style.background_gradient(cmap='Blues').format('{:.2g}')

# 속도를 사용한 PAGA 그림
scv.pl.paga(
    adata,
    basis='umap',
    size=50,
    alpha=0.1,
    min_edge_width=2,
    node_size_scale=1.5
)

전체 워크플로우 스크립트

import scvelo as scv
import scanpy as sc

def run_rna_velocity(adata, n_top_genes=2000, mode='dynamical', n_jobs=4):
    """
    전체 RNA 속도 워크플로우.

    Args:
        adata: 'spliced' 및 'unspliced' 레이어, obsm에 UMAP이 있는 AnnData
        n_top_genes: 속도 계산을 위한 상위 HVG 수
        mode: 'stochastic' (빠름) 또는 'dynamical' (정확)
        n_jobs: 동역학 모델을 위한 병렬 작업 수

    Returns:
        속도 정보가 추가된 처리된 AnnData
    """
    scv.settings.verbosity = 2

    # 1. 전처리 (scVelo 0.3에서는 filter_and_normalize에서 log/HVG 제거됨)
    scv.pp.filter_and_normalize(adata, min_shared_counts=20)
    sc.pp.log1p(adata)
    sc.pp.highly_variable_genes(adata, n_top_genes=n_top_genes, subset=True)

    if 'neighbors' not in adata.uns:
        sc.pp.neighbors(adata, n_neighbors=30)

    scv.pp.moments(adata, n_pcs=30, n_neighbors=30)

    # 2. 속도 추정
    if mode == 'dynamical':
        scv.tl.recover_dynamics(adata, n_jobs=n_jobs)

    scv.tl.velocity(adata, mode=mode)
    scv.tl.velocity_graph(adata)

    # 3. 하위 분석
    if mode == 'dynamical':
        scv.tl.latent_time(adata)
        scv.tl.rank_velocity_genes(adata, groupby='leiden', min_corr=0.3)

    scv.tl.velocity_confidence(adata)
    scv.tl.velocity_pseudotime(adata)

    return adata

AnnData의 주요 출력 필드

워크플로우 실행 후 다음 필드가 추가됩니다:

위치	키	설명
adata.layers	velocity	세포당 유전자당 RNA 속도
adata.layers	fit_t	세포당 유전자당 피팅된 잠재 시간
adata.obsm	velocity_umap	UMAP 상의 2D 속도 벡터
adata.obs	velocity_pseudotime	속도로부터 계산된 가성시간
adata.obs	latent_time	동역학 모델로부터 계산된 잠재 시간
adata.obs	velocity_length	각 세포의 속력
adata.obs	velocity_confidence	세포별 신뢰도 점수
adata.var	fit_likelihood	유전자 수준 모델 피팅 품질
adata.var	fit_alpha	전사율
adata.var	fit_beta	스플라이싱율
adata.var	fit_gamma	분해율
adata.uns	velocity_graph	세포-세포 전이 확률 행렬

속도 모델 비교

모델	속도	정확도	사용 시기
stochastic	빠름	중간	탐색적 분석; 대규모 데이터셋
deterministic	중간	중간	단순 선형 동역학
dynamical	느림	높음	출판 품질 분석; 드라이버 유전자 식별

모범 사례

• 탐색적 분석은 확률론적 모드로 시작하고 최종 분석은 동역학 모드로 전환
• Unspliced read의 좋은 커버리지 필요: 짧은 read(< 100 bp)는 인트론 커버리지를 놓칠 수 있음
• 최소 2,000개 세포: 더 적은 세포에서는 RNA 속도가 noisy함
• 속도는 일관되어야 함: 화살표는 알려진 생물학을 따라야 함; 무작위성은 문제를 나타냄
• k-NN 대역폭 중요: 너무 적은 이웃 → noisy 속도; 너무 많은 이웃 → 과도하게 평활화됨
• 건전성 확인: 루트 세포(전구 세포)는 마커 유전자에 대해 높은 unspliced/spliced 비율을 가져야 함
• 동역학 모델은 구별되는 동역학 상태가 필요함: 명확한 분화 과정에 가장 잘 작동

문제 해결

문제	해결책
Unspliced 레이어 없음	velocyto 재실행 또는 STARsolo를 --soloFeatures Gene Velocyto와 함께 사용
속도 유전자 매우 적음	min_shared_counts 낮추기; 시퀀싱 깊이 확인
무작위로 보이는 화살표	다른 n_neighbors 또는 속도 모델 시도
동역학 모델에서 메모리 오류	n_jobs=1 설정; n_top_genes 줄이기
모든 곳에서 음수 속도	spliced/unspliced 레이어가 바뀌지 않았는지 확인

추가 자료

• scVelo 문서: https://scvelo.readthedocs.io/
• 튜토리얼 노트북: https://scvelo.readthedocs.io/tutorials/
• GitHub: https://github.com/theislab/scvelo
• 논문: Bergen V et al. (2020) Nature Biotechnology. PMID: 32747759
• velocyto (전처리): http://velocyto.org/
• CellRank (운명 예측, scVelo 확장): https://cellrank.readthedocs.io/
• dynamo (대사 표지 대안): https://dynamo-release.readthedocs.io/