Setup the Seurat Object
이 튜토리얼에서는 10X Genomics 에 올라와있는 Peripheral Blood Mononuclear Cells (PBMC) 를 분석할 것 입니다. 2,700개의 세포들이 Illumina NextSeq 500으로 sequence 되었고, raw data 는 여기서 찾으실 수 있습니다.
데이터를 읽는거로 시작을 할 것 인데, [Read10x()](<https://satijalab.org/seurat/reference/read10x>) 함수는 10X 의 cellranger pipeline 에서 나온 값을 읽고, UMI count matrix 를 돌려줍니다. 이 matrix 의 값은 각 세포(열)에서 나온 각 feature 의 molecule 수 (=유전자, 행) 을 나타냅니다.
다음으론 count matrix 을 이용해서 Seurat 객체 (object) 를 생성합니다. 싱글셀 데이터 세트에서 객체는 데이터 (count matrix 와 같은)와 분석 (PCA, clustering results)을 모두 포함하는 container 역할을 합니다. Count matrix 는 pbmc[["RNA"]]@counts 로 저장됩니다.
library(dplyr)
library(Seurat)
library(patchwork)
# PBMC 데이터 로드하기
pbmc.data <- Read10X(data.dir = "../data/pbmc3k/filtered_gene_bc_matrices/hg19/")
# raw 데이터로 Seurat object 만들기
pbmc <- CreateSeuratObject(counts = pbmc.data, project = "pbmc3k", min.cells = 3, min.features = 200)
pbmc
## An object of class Seurat
## 13714 features across 2700 samples within 1 assay
## Active assay: RNA (13714 features, 0 variable features)
Standard pre-processing workflow
여기서부턴 Seurat 에서 standard 한 scRNA-seq 데이터 전처리 workflow 를 포함합니다. 이것은 QC metrics 에 따른 세포 selection 과 filtration, 데이터 normalization 과 scaling, 그리고 highly variable features 를 찾는 것입니다.
QC and selecting cells for further analysis
Seurat 에서는 쉽게 사용자가 정한 기준으로 QC metrics 와 세포를 필터할 수 있습니다. 자주 사용되는 QC metrics 는 아래와 같습니다.
- 각 세포에서 unique 한 유전자가 detect 된 수
- 퀄리티가 낮은 세포나, 빈 droplet 은 유전자가 적습니다
- Cell doublets (droplet 하나에 세포 2개가 들어간 것)와 multiplets (세포 여러개가 들어간 것)은 유전자의 수가 훨씬 높게 나옵니다
- 하나의 세포에서 나온 molecule 의 수
- Mitochondrial genome 에 map 되는 reads 의 퍼센트
- 퀄리티가 낮거나 죽어가고 있는 세포는 mitochondrial contamination 이 높습니다
- mitochondrial QE metrics 는
PercentageFeatureSet()함수를 사용해서 계산할 수가 있고, 이 함수는 percentage of counts originating from a set of features 를 계산합니다 MT-로 시작하는 모든 유전자 세트는 mitochondria 유전자 세트입니다
# The [[ operator can add columns to object metadata. This is a great place to stash QC stats
pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")
아래의 예시에서는 QC metrics 를 visualize 하고, 그걸 사용해서 세포를 filter 했습니다.
- unique feature 의 수가 2,500 보다 높거나 200 보다 낮은 세포를 filter 합니다
- mitochondrial counts 가 >5% 인 세포들도 filter 합니다
# Visualize QC metrics as a violin plot
VlnPlot(pbmc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
# FeatureScatter is typically used to visualize feature-feature relationships, but can be used
# for anything calculated by the object, i.e. columns in object metadata, PC scores etc.
plot1 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "percent.mt")
plot2 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
plot1 + plot2
pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)
Normalizing the data
데이터 세트에서 필요없는 세포들을 지운 후 다음 단계는 데이터를 normalize 하는 것 입니다. 기본적으로 각 세포의 feature expression measurements 을 total expression 으로 normalize 하고, scale factor (기본적으로 10,000) 를 곱하고 이걸 log-transforms 합니다.
여기서 Normalize 된 값은 pbmc[["RNA"]]@data 에 저장됩니다.
pbmc <- NormalizeData(pbmc, normalization.method = "LogNormalize", scale.factor = 10000)
pbmc <- NormalizeData(pbmc)
Identification of highly variable features (feature selection)
다음으로는 데이터세트에서 높은 cell-to-cell variation 이 나오는 features 들의 subset 을 계산합니다 (i.e, 어떤 세포에서는 많이 발현되고, 다른 세포에서는 적게 발현되는 유전자들). Downstream analysis 를 할때 이 유전자들에 초점을 맞추는 것이 싱글셀 데이터세트에서 생물학적 signal 을 찾는거에 도움이 됩니다.
Seurat 의 자세한 절차는 여기서 나오고, 이 버전은 싱글셀 데이터의 mean-variance 관계를 직접 모델링함으로써 이전 버전에서 개선되었으며 이는 FindVariableFeatures() 함수로 구현됩니다. 이 함수에선 기본적으로 하나에 데이터 세트당 2,000개의 유전자가 반환됩니다. 이것은 PCA 와 같은 downstream analysis 에서 쓰입니다.
pbmc <- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 2000)
# Identify the 10 most highly variable genes
top10 <- head(VariableFeatures(pbmc), 10)
# plot variable features with and without labels
plot1 <- VariableFeaturePlot(pbmc)
plot2 <- LabelPoints(plot = plot1, points = top10, repel = TRUE)
plot1 + plot2
Scaling the data
PCA와 같은 dimensional reduction technique 을 하기 전에 필요한 전처리 스텝인 linear transformation (‘scaling”)을 합니다. ScaleData() 함수는:
- 세포에서 평균 expression 이 0이 되도록 각 유전자의 expression 을 옮깁니다
- 각 유전자의 expression 을 scale 해서 세포들의 variance 가 1이 되도록 조정합니다
- 이 스텝은 downstream analysis 에서 equal weight 을 주어서 highly-expressed 된 유전자가 dominate 하지 않도록 합니다
- 여기서 나온 결과는
pbmc[["RNA"]]@scale.data에 저장됩니다
all.genes <- rownames(pbmc)
pbmc <- ScaleData(pbmc, features = all.genes)
Perform linear dimensional reduction
다음으론 scale 된 데이터에 PCA 를 합니다. 기본적으론 위에 나온 variable features 들이 input 으로 사용되지만 features argument 으로 다른 subset 을 고를 수 있습니다.
pbmc <- RunPCA(pbmc, features = VariableFeatures(object = pbmc))
Seurat 은 VizDimReduction(), [DimPlot()](<https://satijalab.org/seurat/reference/DimPlot.html>), [DimHeatmap()](<https://satijalab.org/seurat/reference/DimHeatmap.html>) 등 PCA 를 정의하는 세포와 유전자를 visualization 할 수 있는 유용한 방법을 제공합니다.
# Examine and visualize PCA results a few different ways
print(pbmc[["pca"]], dims = 1:5, nfeatures = 5
## PC_ 1
## Positive: CST3, TYROBP, LST1, AIF1, FTL
## Negative: MALAT1, LTB, IL32, IL7R, CD2
## PC_ 2
## Positive: CD79A, MS4A1, TCL1A, HLA-DQA1, HLA-DQB1
## Negative: NKG7, PRF1, CST7, GZMB, GZMA
## PC_ 3
## Positive: HLA-DQA1, CD79A, CD79B, HLA-DQB1, HLA-DPB1
## Negative: PPBP, PF4, SDPR, SPARC, GNG11
## PC_ 4
## Positive: HLA-DQA1, CD79B, CD79A, MS4A1, HLA-DQB1
## Negative: VIM, IL7R, S100A6, IL32, S100A8
## PC_ 5
## Positive: GZMB, NKG7, S100A8, FGFBP2, GNLY
## Negative: LTB, IL7R, CKB, VIM, MS4A7
VizDimLoadings(pbmc, dims = 1:2, reduction = "pca")
DimPlot(pbmc, reduction = "pca")
DimHeatmap()은 데이터세트에서 heterogeneity 가 어디서 오는지 보여주고 후의 downstream analysis 에서 어떤 PCs 를 포함할지 정할때 유용합니다. 세포와 유전자 둘 다 PCA 스코어에 따라 정리됩니다. cells 를 숫자로 설정하면 ‘extreme’ 한 세포들을 spectrum 의 두 끝에 plot 하게 되는데, 이는 큰 데이터 세트의 plotting 을 dramatically 하게 빨리되게 합니다.
DimHeatmap(pbmc, dims = 1, cells = 500, balanced = TRUE)
DimHeatmap(pbmc, dims = 1:15, cells = 500, balanced = TRUE)
Determine the ‘dimensionality’ of the dataset
scRNA-seq 데이터에서 유전자 하나에서 extensive 하게 technical noise 가 나오는것을 넘어서기위해, Seurat 은 세포들을 PCA 스코어에 따라 cluster를 하고, 각 PC 는 correlated feature set 의 정보가 다 들어간 ‘metafeature’ 입니다. 따라서, top principal components 들은 데이터세트의 robust compression 을 나타냅니다. 그런데 몇개의 components 를 사용해야할까요? 10개? 20개? 100개?
저희는 Macosko et al 에서 JackStraw procedure 에서 inspire 된 resampling test 를 사용했습니다. 데이터의 subset 을 랜덤하게 바꾸고 (1%만) PCA 를 다시 실행해서 features scores 의 ‘null distribution’ 을 만들고 이 process 를 반복했습니다. ‘significant’ 한 PCs 는 p-value features 가 낮게 나온것으로 확인했습니다.
# NOTE: This process can take a long time for big datasets, comment out for expediency. More
# approximate techniques such as those implemented in ElbowPlot() can be used to reduce
# computation time
pbmc <- JackStraw(pbmc, num.replicate = 100)
pbmc <- ScoreJackStraw(pbmc, dims = 1:20)
JackStrawPlot() 함수는 uniform distribution 으로 각 PC 의 p-values 의 distribution (점선) 을 비교하는 visualization 툴입니다. ‘Significant’ 한 PC 들은 낮은 p-value 들로 많이 나옵니다 (점선위에 있는 커브). 여기에선 10-12 PCs 뒤에 확 떨어지는것 같아 보입니다.
JackStrawPlot(pbmc, dims = 1:15)
다른 방법으론 ‘Elbow plot’ 이 있습니다. 이건 variance 의 percentage 에 따라서 principal components 들이 ranking 됩니다 (ElbowPlot()). 이 예시에선 ‘elbow’를 PC9-10에서 볼 수 있고, 이건 대부분의 true signal 들이 첫 10개의 PCs 에 들어간 것을 의미합니다.
ElbowPlot(pbmc)
데이터 세트의 실제 dimensionality 를 확인하는것은 어려울 수 있습니다. 그래서 저희는 이 3개의 방법을 추천드립니다. 첫번째는 더 supervised 한 방법으로 PCs 를 통해서 heterogeneity 의 relevant 한 sources 를 찾는것이고, GSEA 와 함께 사용될 수 있습니다. 두번째는 random null model 에서 나온 statistical test 를 사용하는데 이 방법은, 큰 데이터 세트에서 사용했을때 시간이 많이 들고 clear 한 PC cutoff 를 반환하지 않을 수도 있습니다. 세번째는 자주 쓰이는 방법이고 바로 계산될 수 있습니다. 이 예시에서 3가지 방법 모두 비슷한 결과가 나왔고 PC 7-12 사이에서 cutoff 를 정하는것이 정당했을 수 있습니다.
여기서는 10을 선택했지만, 사용자가 다음 사항을 고려할 것을 권장합니다.
- Dendritic cell 과 NK aficionados 들은 PCs 12 와 13 연관있는 유전자들이 rare한 immue subsets 들이 정의합니다 (i.e. MZB1 은 plasmacytoid DCs 의 marker). 그러나 이 그룹들은 너무 rare해서 prior knowledge 없이 데이터세트의 background noise 에서 구분하기 어렵습니다
- 사용자가 다양한 수의 PC(10, 15 또는 50!)로 다운스트림 분석을 반복할 것을 권장합니다. 보시다시피 결과는 크게 다르지 않습니다
- 사용자에게 이 parameter 를 선택할 때 더 높은 수를 고르는 것으로 실수할 것을 권장합니다. 예를 들어 PC 5개만으로 Downstream analysis를 수행하면 결과에 큰 영향을 미치고 부정적인 영향을 미칩니다.
Cluster the cells
Seurat v3 는 Macosko et al 에서의 initial strategies 을 기반으로 graph-based clustering approach 를 사용합니다. 중요한 것은 clustering analysis 를 실행하는 distance metric (이미 식별된 PCs를 기반으로) 은 그대로 유지된다는 점입니다. 그러나 cellular distance matrix 를 cluster 로 분활하는 방식은 크게 개선되었습니다. 저희의 접근 방식은 scRNA-seq 데이터에 graph-based clustering 을 사용한 최근 manuscript와 [SNN-Cliq, Xu and Su, Bioinformatics, 2015] CyTOF 에 사용한 manuscript [PhenoGraph, Levine et al., Cell, 2015] 에서 크게 inspired 되었습니다. 간단히 말해서, 이러한 방법들은 K-nearest neighbor (KNN) 와 같은 그래프에 비슷한 feature expression 패턴이 나오는 세포끼리 나누고, 이 그래프를 상호 연관성이 높은 ‘quasi-cliques’ 이나 ‘communities’ 로 분활하고 세포를 embed 합니다.
PhenoGraph 에서와 같이, 저희는 우선 PCA 공간에서의 euclidean distance 를 기반으로 KNN 그래프를 만들고, 두 세포 사이의 local neighborhoods 에서의 shared overlap 을 기반으로 edge weights 를 조절합니다 (Jaccard similarity). 이 단계는 FindNeighbors() 함수로 수행되며 이전에 정의된 데이터 세트의 dimensionality (처음 10개 PCs) 가 input 으로 사용됩니다.
세포를 cluster 하기 위해서, standard modularity function 을 목표로 Louvain algorithm (default) 또는 SLM [SLM, Blondel et al., Journal of Statistical Mechanics] 과 같은 modularity optimization techniques 을 적용합니다.
FindClusters() 함수는 이 절차를 구현하며 downstream clustering 의 ‘granularity’ 를 설정하는 resolution parameter 를 포함하며, 이 값이 증가하면 cluster 의 수가 증가합니다. 3000개의 세포의 싱글셀 데이터세트에서는 이 parameter 의 값을 0.4-1.2 사이로 설정하면 좋은 결과가 나옵니다. 데이터 세트가 커지면 optimal resolution 이 높아지는 경우가 많습니다. Idents() 를 사용해서 cluster 를 찾을 수 있습니다.
pbmc <- FindNeighbors(pbmc, dims = 1:10)
pbmc <- FindClusters(pbmc, resolution = 0.5)
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 2638
## Number of edges: 95965
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8723
## Number of communities: 9
## Elapsed time: 0 seconds
# Look at cluster IDs of the first 5 cells
head(Idents(pbmc), 5)
## AAACATACAACCAC-1 AAACATTGAGCTAC-1 AAACATTGATCAGC-1 AAACCGTGCTTCCG-1
## 2 3 2 1
## AAACCGTGTATGCG-1
## 6
## Levels: 0 1 2 3 4 5 6 7 8
Run non-linear dimensional reduction (UMAP/tSNE)
Seurat 은 이러한 tSNE 와 UMAP 같은 여러가지 non-linear dimensional reduction technique 을 제공합니다. 이러한 알고리즘의 목표는 데이터의 기본 manifold 를 학습하고 low-dimensional 한 공간에 유사한 세포를 배치하는 것입니다. 위에서 결정한 graph-based clusters 들의 세포는 dimension reduction plot 에서 co-localize 되어야합니다. UMAP 과 tSNE 의 input 으로는 clustering analysis 에서의 PC를 사용하는 것이 좋습니다.
# If you haven't installed UMAP, you can do so via reticulate::py_install(packages =
# 'umap-learn')
pbmc <- RunUMAP(pbmc, dims = 1:10)
# note that you can set `label = TRUE` or use the LabelClusters function to help label
# individual clusters
DimPlot(pbmc, reduction = "umap")
여기서 객체를 저장해서 위에서 나온 계산을 다시 할 필요 없이 쉽게 로드할 수 있고, 공동작업자와 쉽게 공유할 수 있습니다
saveRDS(pbmc, file = "../output/pbmc_tutorial.rds")
Finding differentially expressed features (cluster biomarkers)
Seurat 을 사용하면 differential expression 을 통해 cluster 를 정의하는 markers를 찾는 데 도움이 됩니다. 기본적으론 다른 모든 세포들과 비교해서 하나의 cluster의 positive와 negative markers를 식별합니다 (ident.1). FindAllMarkers() 는 이 과정을 모든 cluster 에서 자동화하지만, cluster 그룹을 서로간에, 또는 모든 세포에 대해 테스트 할 수도 있습니다.
min.pct argument 는 두 세포 그룹중 하나에서 하나의 feature 가 minimum percentage 로 나와야하고, thresh.test argument 는 이 feature 가 두개의 그룹에서 다르게 발현되어야합니다. 이 두가지를 모두 0으로 설정 할 수 있지만 시간이 증가합니다. 이렇게 하면 highly discriminatory 하지 않은 많은 수의 feature 들도 다 테스트 하기 때문입니다.
이러한 computation 을 빨리하기 위한 또 다른 옵션으로는 max.cell.per.ident 를 설정할 수 있습니다. 이렇게 하면 설정된 것보다 더 많은 세포가 없도록 각 identity class 를 downsample 합니다. 일반적으로 출력 손실이 발생하지만, 속도 증가는 상당할 수 있으며 가장 많이 differentially expressed 된 features 들은 여전히 잘 나옵니다.
# find all markers of cluster 2
cluster2.markers <- FindMarkers(pbmc, ident.1 = 2, min.pct = 0.25)
head(cluster2.markers, n = 5)
## p_val avg_log2FC pct.1 pct.2 p_val_adj
## IL32 2.593535e-91 1.2154360 0.949 0.466 3.556774e-87
## LTB 7.994465e-87 1.2828597 0.981 0.644 1.096361e-82
## CD3D 3.922451e-70 0.9359210 0.922 0.433 5.379250e-66
## IL7R 1.130870e-66 1.1776027 0.748 0.327 1.550876e-62
## LDHB 4.082189e-65 0.8837324 0.953 0.614 5.598314e-61
# find all markers distinguishing cluster 5 from clusters 0 and 3
cluster5.markers <- FindMarkers(pbmc, ident.1 = 5, ident.2 = c(0, 3), min.pct = 0.25)
head(cluster5.markers, n = 5)
## p_val avg_log2FC pct.1 pct.2 p_val_adj
## FCGR3A 2.150929e-209 4.267579 0.975 0.039 2.949784e-205
## IFITM3 6.103366e-199 3.877105 0.975 0.048 8.370156e-195
## CFD 8.891428e-198 3.411039 0.938 0.037 1.219370e-193
## CD68 2.374425e-194 3.014535 0.926 0.035 3.256286e-190
## RP11-290F20.3 9.308287e-191 2.722684 0.840 0.016 1.276538e-186
# find markers for every cluster compared to all remaining cells, report only the positive
# ones
pbmc.markers <- FindAllMarkers(pbmc, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)
pbmc.markers %>%
group_by(cluster) %>%
top_n(n = 2, wt = avg_log2FC)
## # A tibble: 18 x 7
## # Groups: cluster [9]
## p_val avg_log2FC pct.1 pct.2 p_val_adj cluster gene
## <dbl> <dbl> <dbl> <dbl> <dbl> <fct> <chr>
## 1 1.74e-109 1.07 0.897 0.593 2.39e-105 0 LDHB
## 2 1.17e- 83 1.33 0.435 0.108 1.60e- 79 0 CCR7
## 3 0 5.57 0.996 0.215 0 1 S100A9
## 4 0 5.48 0.975 0.121 0 1 S100A8
## 5 7.99e- 87 1.28 0.981 0.644 1.10e- 82 2 LTB
## 6 2.61e- 59 1.24 0.424 0.111 3.58e- 55 2 AQP3
## 7 0 4.31 0.936 0.041 0 3 CD79A
## 8 9.48e-271 3.59 0.622 0.022 1.30e-266 3 TCL1A
## 9 1.17e-178 2.97 0.957 0.241 1.60e-174 4 CCL5
## 10 4.93e-169 3.01 0.595 0.056 6.76e-165 4 GZMK
## 11 3.51e-184 3.31 0.975 0.134 4.82e-180 5 FCGR3A
## 12 2.03e-125 3.09 1 0.315 2.78e-121 5 LST1
## 13 1.05e-265 4.89 0.986 0.071 1.44e-261 6 GZMB
## 14 6.82e-175 4.92 0.958 0.135 9.36e-171 6 GNLY
## 15 1.48e-220 3.87 0.812 0.011 2.03e-216 7 FCER1A
## 16 1.67e- 21 2.87 1 0.513 2.28e- 17 7 HLA-DPB1
## 17 7.73e-200 7.24 1 0.01 1.06e-195 8 PF4
## 18 3.68e-110 8.58 1 0.024 5.05e-106 8 PPBP
Seurat 에는 test.use parameter 로 설정할 수 있는 여러가지 differential expression 테스트가 있습니다. 예를 들어, ROC 테스트는 개별 marker 에 대한 ‘classification power’ 를 반환합니다 (0-랜덤에서 1-perfect 까지).
cluster0.markers <- FindMarkers(pbmc, ident.1 = 0, logfc.threshold = 0.25, test.use = "roc", only.pos = TRUE)
Seurat 은 marker expression 을 visualization 하기 위한 도구도 포함하고 있습니다. VlnPlot() 은 cluster 간 expression percentage 를 보여주고, FeaturePlot() 은 tSNE 또는 PCA plot 에 feature expression 을 visualize 하는 것인데 이 두가지가 visualization에 가장 많이 사용됩니다. 다른 visualization 방법으로는 RidgePlot() , CellScatter() , DotPlot() 도 있습니다.
VlnPlot(pbmc, features = c("MS4A1", "CD79A"))
# you can plot raw counts as well
VlnPlot(pbmc, features = c("NKG7", "PF4"), slot = "counts", log = TRUE)
FeaturePlot(pbmc, features = c("MS4A1", "GNLY", "CD3E", "CD14", "FCER1A", "FCGR3A", "LYZ", "PPBP",
"CD8A"))
DoHeatmap() 은 주어진 세포와 features 에 대해 heatmap 을 생성합니다. 이 경우에는 각 cluster마다 상위 20개의 markers (20개보다 적은경우는 모든 markers) 를 plot 합니다.
pbmc.markers %>%
group_by(cluster) %>%
top_n(n = 10, wt = avg_log2FC) -> top10
DoHeatmap(pbmc, features = top10$gene) + NoLegend()
Assigning cell type identity to clusters
다행히 이 데이터세트의 경우에는, 표준 markers 들로 unbiased clustering 에 쉽게 match 해서 세포 type 을 알 수 있습니다.
Cluster ID, Markers and Cell Type
new.cluster.ids <- c("Naive CD4 T", "CD14+ Mono", "Memory CD4 T", "B", "CD8 T", "FCGR3A+ Mono",
"NK", "DC", "Platelet")
names(new.cluster.ids) <- levels(pbmc)
pbmc <- RenameIdents(pbmc, new.cluster.ids)
DimPlot(pbmc, reduction = "umap", label = TRUE, pt.size = 0.5) + NoLegend()
saveRDS(pbmc, file = "../output/pbmc3k_final.rds")