Major Study./Bioinformatics

Single Cell Analysis Best Practice 정리해보기

sosal 2023. 2. 6. 15:56
반응형

BIML, single cell 강의 들으면서 정리해본 내용입니다.

 

1. Data Format

https://anndata.readthedocs.io/en/latest/index.html

Annotated data: Single cell data를 효율적으로 구성한 데이터 format

 

obsp: (n_obs, n_vars)인 sparse matrix dictionary

일반적으로 n_obs는 Cell의 수이고, n_vars는 Gene의 수


obsm: (n_obs, n_comps)인 sparse matrix dictionary

여기서 n_comps는 구성 요소의 수. -> 차원 감소 또는 클러스터링 알고리즘의 결과를 저장하는 데 사용

(PCA 또는 t-SNE 시각화 등의 2차원 정보 등을 저장)

 

varm: (n_vars, n_vars)인 sparse matrix dictionary

여기서 n_vars는 변수(예: 유전자)의 수

거리 또는 유사성과 같은 변수 간의 관계에 대한 정보를 저장

varp: (n_vars, n_comps)인 sparse matrix dictionary

n_comps: 구성요소의 수, Dimensional reduction 관련 정보 저장

var는 이름, 주석 또는 기타 메타데이터와 같은 변수에 대한 정보를 저장
obs는 세포의 유형, batch, 또는 기타 메타데이터와 같은 관찰에 대한 정보를 저장

 

기본적으로
X축: Gene or Gene-set index

Y축: Cell index

 

 

 

ex) scanpy 패키지를 활용하여, Single cell data import 하기

import scanpy as sc
import anndata as ad

pbmc_5k = sc.read_10x_mtx(
    '/pbmc_5k_v3_filtered_feature_bc_matrix',  
    var_names='gene_symbols',
    cache=True)



scanpy 라이브러리의 read_10x_mtx 함수를 사용하여,
10x Genomics 형식 데이터 파일을 읽는 함수입니다.
그 결과를 pbmc_5k라는 AnnData 개체에 저장합니다.
 
 

 

 

2. Basic QC (Quality Control)

Single cell은 단일 세포 수준에서 유전자 발현을 측정하는 만큼, 정확하게 정량화하기가 어렵다.

발현량이 낮은 Gene 들은 적은 수의 세포로부터 증폭을 하므로, 아예 count가 안되는 유전자들이 매우 많다.

이를 'Drop out' 이라고 한다.

 

=> 따라서 Lowly expressed gene들은 아예 0으로 count가 되는 경우가 많음.

 

하지만, 유전자가 단일 세포에서 발현이 낮다고, 혹은 Drop-out이 되었다고 해서, 전혀 발현되지 않는 것은 아닙니다.

드롭아웃은 단순히 측정 노이즈 또는 제한된 시퀀싱 깊이의 결과일 수 있기 때문에,

분석 결과에 대한 신뢰도를 높이기 위해 단일 실험에서 여러 세포를 분석해야 하는 경우가 많습니다.

 

실제로 Drop out이라고 하더라도, 그 유전자가 아예 발현이 안된다고는 할 수는 없다.

=> 어떤 Signal이 나왔다고 하면, 그것에 대해 확신을 가지고 분석할 수 있지만

=> Drop-out이 낮다고 하여 그것이 진짜로 발현하지 않는다고 확신할 수는 없다.

 

Single cell 하나만 보고 분석을 하면, False가 많을 수 밖에 없음

따라서 Clustering을 통하여 Neighbors를 찾아서

(KNN graph, K-nearest neighbors Graph: 가장 optimial은 아님, 빠르기 떄문에 사용)

 

 

Initial QC

 

Doublets

세포가 많이 읽혔다고 꼭 좋다고는 할 수는 없음

Artificial doublet으로 설명이 가능한 부분에 대해, doublets를 발견할 수 있음

adata.obs['doublet_scores'], adata.obs['predicted_doublets'] = scr.Scrublet(adata.X).scrub_doublets(verbose=False
adata.obs['doublet_scores'].plot.kde()
 
 

Combinatorial filtering 

500~1000개 정도로 적은 수로 읽힌 cell (세포) 같은 경우, Low-coverage dying cells or soup

=> 제거해주는 filtering이 필요함

 

Gene이 6000개 이상 발현되는 유전자도 거의 없음

 

High Mitochondrial contents: 죽어가는 세포 or stressed cells (버리는 경우 많음)

=> 이건 sequencing 기법에 따라 Mitochondrial coverage가 다르기 때문에, cut-off를 잘 설정해야함.

 
 

하지만, 조금 heuristic 하게 해야하는 부분..

 

 

ambient RNAs

세포가 터져 흘러나와, 다른세포랑 합쳐지는 경우

ex: CD3G는 오직 T-cell에서만 발견됨

하지만, 다른 세포에서도 이 유전자가 발현된다? => Ambient RNAs

 

가장 어려운 난제 중 하나임

Linear regression approach를 기반으로, 이를 해결하는 알고리즘들이 나와있음

 

Yang, Shiyi, et al. "Decontamination of ambient RNA in single-cell RNA-seq with DecontX." Genome biology
하지만 Over-correction 문제를 일으킬 수도 있음


이러한 문제 혹은 결과에 대한 판단을 위해서는 여러 Organ들에 대한 지식, 그리고 biology를 조금 더 잘 알아야함..

 

 

Normalization

세포별로 read count 자체가 크게 다르기 때문에,

TMM 같이 library-size, library-composition 등을 고려하는

Bulk-sequencing의 normalization을 사용할 수 없음

 

일반적인 분석의 경우, RPKM 등에서 쓰였던 count normalization 만으로도 충분

발현값이 높을 수록 normalization artifact로 발생하는 경우가 많음

 

단순 Log-transformation을 많이 사용

Xlog = log(Xnorm + 1)

 

Selecting - Highly variable genes

Not informative: 모든 cell에서 발견되는 or 어디에서도 발견되지 않는 expression,

Informative: 특정 Cell type에서 특이적으로 발견되는 expression

 

Scaling: Centralise gene expression with zero mean

한 데이터라도 scaling이 잘못 되어있으면, 전체 데이터가 망가지므로

꼭 데이터의 range를 확인하여 어떤 작업을 수행하여 만들어진 데이터인지 체크해야함

 

 

 

 

Dimension reduction

Genes: 30,000

Cell: up to 1,000,000

 

차원이 너무 크기에, 모든 세포간의 정확한 Distance를 구하는건 사실상 불가능

따라서 Dimensional reduction이 필수가 됨

 

PCA -> PC50개만 사용하여, T-cell에 zoom-in 하여 본다? 잘 안될 수 있음

따라서 얼마나 많은 Dimension을 활용하여야 하느냐 -> heuristic하게 search 해야함

 

Dimension reduction 알고리즘에 대한 뭔가 추상적인 설명?

t-SNE: 밀어내는걸 강조

Diffusion map: 당기는걸 강조

UMAP: 그 가운데 쯤 있음.. 그래서 많이 씀

 

# PCA
sc.tl.pca(adata)
sc.pl.pca_variance_ratio(adata, log=True)
 
# UMAP
sc.pp.neighbors(adata, n_pcs=10)
sc.tl.umap(adata)
sc.pl.umap(adata, color = 'dataset')
 
sc.pp.neighbors(adata, n_pcs=50)
sc.tl.umap(adata)
sc.pl.umap(adata, color = 'dataset')

 

UMAP 또한, 몇개의 PCA component를 활용하느냐에 따라 다른 결과가 나옴.

(근데 이걸 원하는 그림이 잘 나올때까지 계속 반복하면, 그게 사실 데이터에 대한 overfit 아닌가 싶네요)

 

 

3. Annotation Using Marker Genes

Marker gene의 expression을 UMAP상에서 혹은 Heatmap이나 dotplot을 통해서 확인 한 후 annotation을 진행
Marker gene의 경우 reference paper를 참고하거나 Azimuth, PanglaoDB 등의 marker gene database를 활용
Azimuth: https://azimuth.hubmapconsortium.org/
PanglaoDB: https://panglaodb.se/

 

PanglaoDB - A Single Cell Sequencing Resource For Gene Expression Data

PanglaoDB is a database for the scientific community interested in exploration of single cell RNA sequencing experiments from mouse and human. We collect and integrate data from multiple studies and present them through a unified framework. Adapted from th

panglaodb.se

 

Cell annotation은 계층화된 Annotation이 필요함.

(전체를 들여봤다가, 세부적으로 zoom-in 하여 들어가는 방식)

 

예를 들어, T-Cell vs B-Cell 구분과, T-cell 안에서의 subtype 구분은 서로 다른 계층에서의 Annotation이 필요함

marker_genes = {
    'CD4T' : 'IL7R, MAL, LTB, CD4, TRAC, CD3D, CD3G'.split(', '),
    .... (생략)
}
 
sc.pl.umap(adata, color = marker_genes['CD4T'], cmap='OrRd'

 

- Leiden: Community Detection Algorithm for Graphs

 

sc.tl.leiden(adata, resolution = 0.3, key_added='leiden_0.3')
sc.tl.leiden(adata, resolution = 0.5, key_added='leiden_0.5')
sc.tl.leiden(adata, resolution = 1.0, key_added='leiden_1.0')
sc.tl.leiden(adata, resolution = 2.0, key_added='leiden_2.0')
 
위에서 언급했 듯, 어떤 Resolution level에서 clustering 할 지에 따라서
세포 군집들의 결과는 당연히 다른 양상을 띄게 된다.

 

동일한 UMAP 결과에서 visualization 결과를 그려준다.

 

sc.pl.dotplot(adata, marker_genes, groupby = 'leiden_0.3', swap_axes = True)
마커에와 Cluster 군집에 따른 발현 양상

 

 

** Graph based cell clustering 좋은 자료

https://nbisweden.github.io/workshop-scRNAseq/labs/compiled/scanpy/scanpy_04_clustering.html

 

scanpy_04_clustering

**Your turn** By now you should know how to plot different features onto your data. Take the QC metrics that were calculated in the first exercise, that should be stored in your data object, and plot it as violin plots per cluster using the clustering meth

nbisweden.github.io

 

 

 

 

rank_genes_groups function은 cluste와 나머지에 그룹에 대한 Differential Expressed Gene (DEG) 분석을 해줌으로써, 각 cluster에 특이적 발현 유전자를 꼽아줍니다.

 

sc.tl.rank_genes_groups(adata, 'leiden_0.3', method='wilcoxon')
sc.pl.rank_genes_groups(adata, n_genes=25, sharey=False)
sc.pl.rank_genes_groups_heatmap(adata, n_genes=25, groupby="leiden_0.3", show_gene_labels=True)
 

 

=> 하지만 이렇게 하면 안된다!!

pseudobulk DE methods 가 제일 좋다고 알려져있음

과거 Bulk RNASeq의 distribution을 잘 이해하여, library-size, library-composition을 고려한 DEG 분석이 매우 잘 정립되어 있음. 하지만 Single-cell RNASeq은 Drop-out 등 데이터의 분포가 불안정 하기에, cell 단위의 resolution에서는 에러가 많을 수 밖에 없음.

따라서, Cell 단위로 나뉘어져 있는 Raw count를 Bulk-seq처럼 합쳐서 분석하는 것이 가장 robust하다고 위 논문에서 이미 알려져있음.

 

 

 

4. Trajectory inference

Clustering은 expression을 기반으로, 비슷한 cell 들을 clustering 하는 것이 목적이었다면,

Trajectory inference는 세포들의 발달에 대한 궤적 (trajectory)을

분자 이벤트의 '순서'를 추론하여 경로를 생성하는 방법입니다.

 

어디까지나 이는 Computational method를 통해 reconstructed 된 것이기 때문에,

'Exact representation of the developmental biology'가 아니라는 것을 명심해야합니다.

 

특히, 위 그림을 직접 생성했던 박종은 교수님도,

위 발달과정에서 세포들이 만들어지는 biological sequence를 따르는 'UMAP'을 그리기 위해

수백번의 위 그림을 그렸다는 이야기를 말씀하셨습니다.

 

 

 

5. Data integration, Batch effect

- HUMAN CELL ATLAS - Data portal

Public database로 공개되어있는 Single cell data는 gene set이 굉장히 다르기에, 합쳐서 쓰기가 쉽지 않습니다.

 

Intersection이 매우 좋지만, 이는 너무 적은 Gene set이 될 수 있고,

Union -> 없는 데이터를 0으로 채우는 방식으로 합치는 것은 최악이라고 한다.

 

따라서 적절하게 없는 Gene set은 발현이 없어서 날렸다는 가정하에,

특정 데이터 셋으로 filter하고, 없다면 0으로 채우는 걸 적절하게 하는것이 좋다고 한다.

 

HUMAN CELL ATLAS, sanger institute에서 Raw data에서부터 mapping을 처음부터 다시하여

(STAR alignment) gene set을 맞춰주는 작업을 하고있다고 합니다.

 

 

Single cell - Batch effect

 

Batch effect는 어떤 실험 환경이냐에 따라서, biology에 의해 발생하는 차이가 아닌,

실험을 하면서 만들어지는 차이를 뜻합니다.

(예: Reagents, Protocols, Meaturement platforms)

 

이것은 systematic difference를 발생시키므로, Batch effect에 의한 차이는

실제 우리가 scientific하게 발견하고자 하는 결과에 매우 악영향을 끼칠 것입니다.

 

Bulk는 많은 수의 average 결과이기에, single cell analysis이 bulk보다 큰 경향을 보이긴 하는데,

실제로 Single-cell analysis의 batch effect는 무지막지하게 크지는 않습니다.

Graph-based clustering에서 (주로 Graph) 주로 한 세포의 기준으로 주변에 대해 살펴보기 때문에,

batch effect가 굉장히 커보이게 됩니다.

 

 

Batch-effect는 여러 단계에서 수행할 수 있습니다.

Gene count 자체에서 batch-effect correction (Seurat Integrate)

PCA에서 batch-effect correction (Harmony)

...

 

Batch correction이 잘 합쳐졌는지, 평가하기가 굉장히 어렵다는게 challenging.

무엇을 봐야하냐?

 

- Measuring integration

서로 다른 batch끼리 잘 섞여있느냐?

(한 점으로 모두다 뭉쳐버리면 안됨, 따라서 이것만으로는 좋은 metric이라 할 수 없음)

 

- Measuring Cell type preservation

 

cell type이 없는데 어떻게 이걸 하나요?

-> Individual data 안에서 clustering 하면서, 서로간의 mixture가 잘 생기느냐? 로 간접적으로 볼 수 있다.

 

요즘 외부데이터가 많아서, 외부 데이터를 기반으로 clusteing 및 cell type prediction 이후,

정해져있는 cell type을 기반으로 합칠 수 있다..

 

등등...

 

rough한 정답을 만들어 놓고 batch-correction을 하는 작업들이 필요함

 

 

* Bulk에서 주로 쓰던 Linear regression, 근데 이거 scRNA에서는 잘 쓰지 않음.

single cell에서는 batch effect가 생각보다 크지 않아서,

coeff를 조금 낮추는 Regularized linear regression을 활용하여 사용

 

- BBKNN (UMAP based)

sc.external.pp.bbknn(adata, batch_key='dataset')
sc.tl.umap(adata)
sc.pl.umap(adata, color='dataset')
adata.obsm['X_umap_bbknn'] = adata.obsm['X_umap'].copy()
 

- Harmony (Cluster level에서의 batch correction)

sc.external.pp.harmony_integrate(adata, key = 'dataset', max_iter_harmony = 50)
sc.pp.neighbors(adata, n_pcs=15, use_rep = 'X_pca_harmony')
sc.tl.umap(adata)
sc.pl.umap(adata, color = 'dataset')

 

- SCVI (Single-Cell Variational Inference, 딥러닝 기반)

import scvi
scvi.model.SCVI.setup_anndata(adata, layer="counts", batch_key="dataset")
vae = scvi.model.SCVI(adata, n_layers=2, n_latent=30)
vae.train(max_epochs = 50)
adata.obsm["X_scVI"] = vae.get_latent_representation()

좌: Batch correction 이전의 UMAP,    우: BBKNN & Harmony 알고리즘으로 Batch correction

 

딥러닝 기반의 Batch correction: SCVI 메소드의 결과

 

박종은 교수님의 Best practice

Linear regression인데, 외부 Label을 가져와서 사용

(Gene expression ~ Batch + Cell Type) += 추가적인 다른단계에서 Batch effect correction

=> 리비전할때 고생할 확률이 있다? ㅋㅋ

 

딥러닝 기반

- Deep Generative modeling for single cell transcriptomics

https://docs.scvi-tools.org/ 

 

Documentation

scvi-tools(single-cell variational inference tools) is a package for end-to-end analysis of single-cell omics data primarily developed and maintained by the Yosef Lab at UC Berkeley. scvi-tools has...

docs.scvi-tools.org

 

가장 중요한것, Garbage in -> Garbage out.

QC를 잘 끝낸 이후에서야 이런 것들을 적용할 수 있다.

 

 

 

6. Cell Annotation

import celltypist
from celltypist import models

 

celltypist는 tutorial이 굉장히 잘 되어있음

https://www.celltypist.org/

 

Automated cell type annotation for scRNA-seq datasets

CellTypist provides automated cell type annotation for scRNA-seq datasets

www.celltypist.org

 

Celltypist는 Bulk-seq에서의 CPM처럼, 전체 sum이 10^4로 되어있음

마지막 Spatial transcriptomics 등에 대한 썰 이후 끝.

 

- Pre-built model

models.download_models(force_update = True)
models.models_description()
 
- pre-built 된 모델들의 목록
Immune_All_Low.pkl immune sub-populations combined from 20 tissue...
Immune_All_High.pkl immune populations combined from 20 tissues of...
Adult_Mouse_Gut.pkl cell types in the adult mouse gut combined fro...
Autopsy_COVID19_Lung.pkl cell types from the lungs of 16 SARS-CoV-2 inf...
COVID19_Immune_Landscape.pkl immune subtypes from lung and blood of COVID-1...
model = models.Model.load(model = 'Immune_All_Low.pkl')
predictions = celltypist.annotate(adata, model = 'Immune_All_Low.pkl', majority_voting = True)
adata = predictions.to_adata()
 
 

Batch correction evaluation & Celltypist

sc.pl.embedding(adata, basis = 'X_umap_scvi', color = ['dataset','majority_voting'], ncols=1)

 

* Problem

환자 데이터와 정상 샘플을 분류할 때, 실제로 위와 같이 완벽하게 겹치게 되는것이

사실 Over-correction이 될 수 있다..

과연 Biological truth는 무엇일까? 에 대해 계속 고민해야함..