## scRNA-seq 분석 워크플로우
Single-cell RNA-seq(scRNA-seq) 데이터 분석은 원시 데이터로부터 생물학적 인사이트를 얻기까지 다음 6단계로 이루어집니다.
---
### 1. 데이터 전처리 (Preprocessing)
- **Quality Control (QC)**
- 세포별 총 UMI count, 검출 유전자 수, 미토콘드리아 유전자 비율 계산 후 임계값 기반 저품질 셀 제거
- 예: counts < 200, genes < 200, mito% > 5%
- **Doublet** 탐지: Scrublet, DoubletFinder 사용 가능
- **Filtering**
- 세포·유전자 필터링 (min_cells, min_genes)
- ambient RNA 보정, 빈 드롭릿 제거
- **Normalization & Log-변환**
- `normalize_total(target_sum=1e4)` → `log1p()`
```python
import scanpy as sc
adata = sc.read_10x_mtx("data/filtered_feature_bc_matrix/")
sc.pp.filter_cells(adata, min_genes=200)
sc.pp.filter_genes(adata, min_cells=3)
adata.var['mt'] = adata.var_names.str.startswith('MT-')
sc.pp.calculate_qc_metrics(adata, qc_vars=['mt'], inplace=True)
adata = adata[adata.obs.pct_counts_mt < 5]
sc.pp.normalize_total(adata, target_sum=1e4)
sc.pp.log1p(adata)
```
---
### 2. 클러스터링 (Clustering)
1. **Highly Variable Genes (HVGs)** 선택
- 상위 1,000–3,000개 유전자
2. **차원 축소**
- PCA (`n_comps=50`)
3. **그래프 생성**
- kNN 그래프 기반 (`pp.neighbors`)
4. **Leiden** 또는 Louvain 알고리즘으로 군집화
```python
sc.pp.highly_variable_genes(adata, n_top_genes=2000)
adata = adata[:, adata.var.highly_variable]
sc.pp.pca(adata, n_comps=50)
sc.pp.neighbors(adata, n_neighbors=10, n_pcs=50)
sc.tl.leiden(adata, resolution=0.5)
print(adata.obs['leiden'].value_counts())
```
---
### 3. 세포 타입 주석 (Annotation)
- **마커 유전자 식별**
- `rank_genes_groups(groupby='leiden', method='wilcoxon')`
- **수동/자동 주석**
- CD3D → T cell, CD19 → B cell 등
- SingleR, CellTypist, Azimuth 등 자동 주석 도구 활용 가능
```python
sc.tl.rank_genes_groups(adata, groupby='leiden', method='wilcoxon')
# 상위 5개 마커 출력
print(adata.uns['rank_genes_groups']['names']['0'][:5])
# 수동 매핑
mapping = {'0':'T cell','1':'B cell','2':'Monocyte'}
adata.obs['cell_type'] = adata.obs['leiden'].map(mapping).fillna('Unknown')
```
---
### 4. 시각화 (Visualization)
- **PCA plot**: `sc.pl.pca` (PC1 vs PC2)
- **UMAP**: 빠르고 구조 보존 우수
- `sc.tl.umap()` → `sc.pl.umap(color=['leiden','cell_type'])`
- **t-SNE**: 파라미터 민감도 주의
```python
sc.tl.umap(adata)
sc.pl.umap(adata, color=['leiden'])
```
---
### 5. 통합 분석 (Integration)
1. **배치 통합 (Multi-sample)**
- Seurat CCA, MNN, Harmony, BBKNN 등
- **Scanpy 예시 (BBKNN)**: `sce.pp.bbknn(adata, batch_key='batch')`
2. **다중 모달 (Multi-modal)**
- scRNA+scATAC: WNN, MOFA+, TotalVI, LIGER 등
```python
import scanpy.external as sce
adata1 = sc.read_h5ad('sample1.h5ad')
adata2 = sc.read_h5ad('sample2.h5ad')
ad = adata1.concatenate(adata2, batch_key='batch')
sce.pp.bbknn(ad, batch_key='batch')
ad.obs['batch']
sc.tl.umap(ad); sc.tl.leiden(ad)
sc.pl.umap(ad, color=['batch','leiden'])
```
---
### 6. 추가 분석: RNA Velocity
`RNA velocity`는 각 세포의 **미래 전사체 상태**를 예측하는 기법으로, 시간 정보(time dynamics)를 갖지 않는 scRNA-seq 데이터에 **동적(flow) 정보**를 부여합니다. 주요 아이디어는 **미성숙(unspliced) 전사체**와 **성숙(spliced) 전사체** 비율의 차이를 통해 **유전자 발현 변화 방향과 속도**를 추정하는 것입니다.
- **기본 개념**
- 전사 초기에는 pre-mRNA(미성숙)가 우선 생성되고, 이후 splicing 과정을 거쳐 mature mRNA(성숙)로 전환됩니다.
- 세포별 spliced vs unspliced transcript abundance 비교를 통해, 특정 유전자의 발현이 **증가 중인지 감소 중인지**를 파악합니다.
- 각 유전자별 속도를 종합하여 세포 단위의 **벡터(velocity vector)** 를 계산하고, 이를 UMAP 같은 저차원 공간에 시각화합니다.
- **주요 도구**
- **Velocyto**: `.loom` 포맷으로 spliced/unspliced count matrix를 생성하는 초기 분석 파이프라인 툴
- **scVelo**: Python 기반 확장 라이브러리로, stochastic model 및 dynamical model을 통해 **더 정교한 velocity 추정**과 **latent time** 계산 기능을 제공
```python
import scvelo as scv
# loom 파일 로드 (velocyto로 생성된 spliced/unspliced 정보 포함)
adata_velo = scv.read('sample.loom', cache=True)
# 전처리: 유전자 필터링 및 정규화
scv.pp.filter_and_normalize(adata_velo, min_shared_counts=20, n_top_genes=2000)
# 모멘트 계산 (주성분 + 이웃 그래프 기반)
scv.pp.moments(adata_velo, n_pcs=30, n_neighbors=30)
# Velocity 추정 (stochastic 모델)
scv.tl.velocity(adata_velo, mode='stochastic')
# Velocity 그래프 생성
scv.tl.velocity_graph(adata_velo)
# UMAP embedding 위에 스트림 플롯으로 시각화
scv.pl.velocity_embedding_stream(
adata_velo, basis='umap', color='clusters', legend_loc='right'
)
# (옵션) Dynamical model로 속도 복구 및 latent time 계산
scv.tl.recover_dynamics(adata_velo) # 유전자별 dynamics 파라미터 추정
scv.tl.velocity(adata_velo, mode='dynamical')
scv.tl.velocity_graph(adata_velo)
scv.tl.velocity_pseudotime(adata_velo) # latent time 계산
scv.pl.scatter(
adata_velo, color='velocity_pseudotime', cmap='gnuplot', size=40
)
```
- **분석 팁**
- **Gene selection**: velocity 분석 시에도 HVG 기반 필터링이 필요하며, 너무 low-count 유전자는 제거
- **모델 선택**: `stochastic`은 빠르지만, 데이터에 동적 변화를 상세히 반영하려면 `dynamical` 모드를 권장
- **시각화 옵션**: `velocity_embedding_stream` 외에 `velocity_embedding_grid`, `velocity_embedding` 함수도 사용 가능
- **Latent time**: 세포의 분화 궤적(pseudotime)과 비교하여 **time-resolved trajectory** 분석에 활용
---
`RNA velocity`를 통해 단일 세포의 **분화 흐름**, **상태 전이 경로**를 동적으로 파악할 수 있습니다. Velocyto와 scVelo를 조합하면, 비정형적 scRNA-seq 데이터에도 **시간적 축**을 추가해 더 풍부한 생물학적 해석이 가능합니다.