Guided clustering tutorial – 2,700 PBMCs

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)
https://satijalab.org/seurat/articles/pbmc3k_tutorial_files/figure-html/qc2-1.png
# 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
https://satijalab.org/seurat/articles/pbmc3k_tutorial_files/figure-html/qc2-2.png
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
https://satijalab.org/seurat/articles/pbmc3k_tutorial_files/figure-html/var_features-1.png

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")
https://satijalab.org/seurat/articles/pbmc3k_tutorial_files/figure-html/pca_viz-1.png
DimPlot(pbmc, reduction = "pca")
https://satijalab.org/seurat/articles/pbmc3k_tutorial_files/figure-html/pca_viz-2.png

DimHeatmap()은 데이터세트에서 heterogeneity 가 어디서 오는지 보여주고 후의 downstream analysis 에서 어떤 PCs 를 포함할지 정할때 유용합니다. 세포와 유전자 둘 다 PCA 스코어에 따라 정리됩니다. cells 를 숫자로 설정하면 ‘extreme’ 한 세포들을 spectrum 의 두 끝에 plot 하게 되는데, 이는 큰 데이터 세트의 plotting 을 dramatically 하게 빨리되게 합니다.

DimHeatmap(pbmc, dims = 1, cells = 500, balanced = TRUE)
https://satijalab.org/seurat/articles/pbmc3k_tutorial_files/figure-html/single-heatmap-1.png
DimHeatmap(pbmc, dims = 1:15, cells = 500, balanced = TRUE)
https://satijalab.org/seurat/articles/pbmc3k_tutorial_files/figure-html/multi-heatmap-1.png

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)
https://satijalab.org/seurat/articles/pbmc3k_tutorial_files/figure-html/jsplots-1.png

다른 방법으론 ‘Elbow plot’ 이 있습니다. 이건 variance 의 percentage 에 따라서 principal components 들이 ranking 됩니다 (ElbowPlot()). 이 예시에선 ‘elbow’를 PC9-10에서 볼 수 있고, 이건 대부분의 true signal 들이 첫 10개의 PCs 에 들어간 것을 의미합니다.

ElbowPlot(pbmc)
https://satijalab.org/seurat/articles/pbmc3k_tutorial_files/figure-html/elbow_plot-1.png

데이터 세트의 실제 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")
https://satijalab.org/seurat/articles/pbmc3k_tutorial_files/figure-html/tsneplot-1.png

여기서 객체를 저장해서 위에서 나온 계산을 다시 할 필요 없이 쉽게 로드할 수 있고, 공동작업자와 쉽게 공유할 수 있습니다

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"))
https://satijalab.org/seurat/articles/pbmc3k_tutorial_files/figure-html/markerplots-1.png
# you can plot raw counts as well
VlnPlot(pbmc, features = c("NKG7", "PF4"), slot = "counts", log = TRUE)
https://satijalab.org/seurat/articles/pbmc3k_tutorial_files/figure-html/markerplots-2.png
FeaturePlot(pbmc, features = c("MS4A1", "GNLY", "CD3E", "CD14", "FCER1A", "FCGR3A", "LYZ", "PPBP",
    "CD8A"))
https://satijalab.org/seurat/articles/pbmc3k_tutorial_files/figure-html/markerplots-3.png

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()
https://satijalab.org/seurat/articles/pbmc3k_tutorial_files/figure-html/clusterHeatmap-1.png

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()
https://satijalab.org/seurat/articles/pbmc3k_tutorial_files/figure-html/labelplot-1.png
saveRDS(pbmc, file = "../output/pbmc3k_final.rds")