原文見Seurat - Guided Clustering Tutorial, Compiled: April 17, 2020
#1 Seurat安裝
install.packages("Seurat")
#2 數據下載
Peripheral Blood Mononuclear Cells (PBMC) 是10X Genomics dataset page提供的一個數據，包含2700個單細胞，出自Illumina NextSeq 500平臺。
PBMCs是來自健康供體具有相對少量RNA（around 1pg RNA/cell）的原代細胞。在Illumina NextSeq 500平臺，檢測到2700個單細胞，每個細胞獲得69000 reads。
tar -xvf pbmc3k_filtered_gene_bc_matrices.tar
文件夾下包含3個文件
barcodes.tsv
genes.tsv
matrix.mtx
matrix.mtx：matrix.mtx 是 MatrixMarket格式文件；更多內容見：http://math.nist.gov/MatrixMarket/formats.html
文件中儲存非零值;
注釋使用%標記；
第一行包含文件中總行數，總列數，總的記錄數
每行中提供記錄的所處的行號和列號，已經記錄的內容
?head filtered_gene_bc_matrices/hg19/matrix.mtx
%%MatrixMarket matrix coordinate real general
%
32738 2700 2286884
32709 1 4
32707 1 1
32706 1 10
32704 1 1
32703 1 5
32702 1 6
32700 1 10
#3 數據導入
##3.1 Read10X()函數可以自動讀入和解析數據。
library(dplyr)
library(Seurat)
library(patchwork)
#讀取PBMC數據
pbmc.data <- Read10X(data.dir = "../data/pbmc3k/filtered_gene_bc_matrices/hg19/")
#查看數據
dim(pbmc.data)
# 32738 2700
pbmc.data[c("CD3D", "TCL1A", "MS4A1"), 1:30]
# 3 x 30 sparse Matrix of class "dgCMatrix"
# CD3D 4 . 10 . . 1 2 3 1 . . 2 7 1 . . 1 3 . 2 3 . . . . . 3 4 1 5
# TCL1A . . . . . . . . 1 . . . . . . . . . . . . 1 . . . . . . . .
# MS4A1 . 6 . . . . . . 1 1 1 . . . . . . . . . 36 1 2 . . 2 . . . .
#.表示0
summary(colSums(pbmc.data))
# Min. 1st Qu. Median Mean 3rd Qu. Max.
# 548 1758 2197 2367 2763 15844
#查看每個細胞有多少基因被檢測到
at_least_one <- apply(pbmc.data, 2, function(x) sum(x>0))
hist(at_least_one, breaks = 100,
main = "Distribution of detected genes",
xlab = "Genes with at least one tag")
hist(colSums(pbmc.data),
breaks = 100, main = "Expression sum per cell",
xlab = "Sum expression")
##3.2 使用pbmc數據初始化Seurat對象
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)
head(pbmc$RNA@data[,1:5])
6 x 5 sparse Matrix of class "dgCMatrix"
AAACATACAACCAC-1 AAACATTGAGCTAC-1 AAACATTGATCAGC-1 AAACCGTGCTTCCG-1 AAACCGTGTATGCG-1
AL627309.1 . . . . .
AP006222.2 . . . . .
RP11-206L10.2 . . . . .
RP11-206L10.9 . . . . .
LINC00115 . . . . .
NOC2L . . . . .
#4 數據預處理
這部分是基于數據質控方法，標準化和檢測到的變化基因對數據進行篩選。
##4.1 對細胞的質控
可以參考文章：Classification of low quality cells from single-cell RNA-seq data
- 單個細胞中檢測到單個基因的數目
- 低質量的細胞以及空泡油滴中一般檢測到很少的基因
- 包含多個細胞的油滴會檢測到異常多的基因
- 類似，在一個單細胞中的基因總數
- 檢測到線粒體的基因數目百分比
- 低質量/垂死細胞通常會有線粒體污染
- 使用PercentageFeatureSet函數函數可以計算線粒體QC，計算線粒體基因所占百分比
- 基因名以MT- 開始的基因定義為線粒體基因
#線粒體基因占比計算
pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")
head(pbmc@meta.data, 5)
orig.ident nCount_RNA nFeature_RNA percent.mt
AAACATACAACCAC-1 pbmc3k 2419 779 3.0177759
AAACATTGAGCTAC-1 pbmc3k 4903 1352 3.7935958
AAACATTGATCAGC-1 pbmc3k 3147 1129 0.8897363
AAACCGTGCTTCCG-1 pbmc3k 2639 960 1.7430845
AAACCGTGTATGCG-1 pbmc3k 980 521 1.2244898
#畫圖查看基因數目, UMI數目, 線粒體基因占比
VlnPlot(pbmc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
#查看基因數目, 線粒體基因占比與UMI數目的關系
plot1 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "percent.mt")
plot2 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
plot1 + plot2
#質控
- 篩選檢測到基因數目超過2500或低于200的細胞
- 單個細胞中線粒體基因數目占比超過>5%
pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)
##4.2 數據標準化
默認使用數據標準化方法是LogNormalize, 每個細胞總的表達量都標準化到10000，然后log取對數；結果存放于pbmc[["RNA"]]@data。
#標準化前，每個細胞總的表達量
hist(colSums(pbmc$RNA@data),
breaks = 100,
main = "Total expression before normalisation",
xlab = "Sum of expression")
pbmc <- NormalizeData(pbmc, normalization.method = "LogNormalize", scale.factor = 10000)
#標準化后，每個細胞總的表達量
hist(colSums(pbmc$RNA@data),
breaks = 100,
main = "Total expression after normalisation",
xlab = "Sum of expression")
##4.3 變化基因鑒定
鑒定在細胞間表達高度變化的基因，后續(xù)研究需要集中于這部分基因。Seurat內置的FindVariableFeatures()函數，首先計算每一個基因的均值和方差，并且直接模擬其關系。默認返回2000個基因。
pbmc <- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 2000)
# 10個表達變化最為劇烈的基因
top10 <- head(VariableFeatures(pbmc), 10) #head(pbmc$RNA@var.features,10)
# "PPBP" "LYZ" "S100A9" "IGLL5" "GNLY" "FTL" "PF4" "FTH1" "GNG11" "S100A8"
# 畫出表達變化的基因，從而觀察其分布
plot1 <- VariableFeaturePlot(pbmc)
# 畫出表達變化的基因，標記前10個基因
plot2 <- LabelPoints(plot = plot1, points = top10, repel = TRUE)
plot1
plot2
##4.4 數據縮放
線性轉換縮放數據，ScaleData()函數可以實現此功能。
最終每個基因均值為0，方差為1。
結果存放于pbmc[["RNA"]]@scale.data。
all.genes <- rownames(pbmc)
pbmc <- ScaleData(pbmc, features = all.genes)
設置參數features是因為ScaleData默認處理前面鑒定的差異基因。這一步怎么做都不會影響到后續(xù)pca和聚類，但是會影響做熱圖。
移除影響方差的因素
pbmc <- ScaleData(pbmc, vars.to.regress = "percent.mt")
#5 線性降維分析
##5.1 PCA
對縮放后的數據進行PCA分析，默認使用前面鑒定表達變化大的基因。使用features參數可以重新定義數據集。
pbmc <- RunPCA(pbmc, features = VariableFeatures(object = pbmc))
VizDimReduction, DimPlot, 和 DimHeatmap可以從基因或細胞角度可視化pca結果
#查看對每個主成分影響比較大的基因集
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",split.by = 'ident')
DimHeatmap繪制基于單個主成分的熱圖，細胞和基因的排序都是基于他們的主成分分數。對于數據異質性的探索是很有幫助的，可以幫助用戶選擇用于下游分析的主成分維度。
DimHeatmap(pbmc, dims = 1, cells = 500, balanced = TRUE)
展示多個主成分
DimHeatmap(pbmc, dims = 1:15, cells = 500, balanced = TRUE)
##5.1 數據維度
為了避免單個基因影響，Seurat聚類細胞時使用pca結果。首先需要確定的是使用多少個主成分用于后續(xù)分析。常用有兩種方法，一種是基于零分布的統(tǒng)計檢驗方法，這種方法耗時且可能不會返回明確結果。另一種是主成分分析常用的啟發(fā)式評估。
- JackStraw()
在JackStraw()函數中, 使用基于零分布的置換檢驗方法。隨機抽取一部分基因(默認1%)然后進行pca分析得到pca分數，將這部分基因的pca分數與先前計算的pca分數進行比較得到顯著性p-Value，。根據主成分（pc）所包含基因的p-value進行判斷選擇主成分。最終的結果是每個基因與每個主成分的關聯的p-Value。保留下來的主成分是那些富集小的p-Value基因的主成分。
處理大數據時會花費大量時間；ElbowPlot()內置了一些其它的方法可以減少運行時間。
pbmc <- JackStraw(pbmc, num.replicate = 100)
pbmc <- ScoreJackStraw(pbmc, dims = 1:20)
JackStrawPlot()函數提供可視化方法，用于比較每一個主成分的p-value的分布，虛線是均勻分布；顯著的主成分富集有小p-Value基因，實線位于虛線左上方。下圖表明保留10個pca主成分用于后續(xù)分析是比較合理的。
JackStrawPlot(pbmc, dims = 1:15)
- ElbowPlot
ElbowPlot(pbmc)
啟發(fā)式評估方法生成一個Elbow plot圖。在圖中展示了每個主成分對數據方差的解釋情況（百分比表示），并進行排序。根據自己需要選擇主成分，圖中發(fā)現第9個主成分是一個拐點，后續(xù)的主成分(PC)變化都不大了。
注意：鑒別數據的真實維度不是件容易的事情；除了上面兩種方法，Serat官當文檔還建議將主成分(數據異質性的相關來源有關)與GSEA分析相結合。Dendritic cell 和 NK aficionados可能識別的基因與主成分 12 和 13相關，定義了罕見的免疫亞群 (i.e. MZB1 is a marker for plasmacytoid DCs)。如果不是事先知道的情況下，很難發(fā)現這些問題。
Serat官當文檔因此鼓勵用戶使用不同數量的PC（10、15，甚至50）重復下游分析。其實也將觀察到的，結果通常沒有顯著差異。因此，在選擇此參數時，可以盡量選大一點的維度，維度太小的話對結果會產生不好的影響。
#6 細胞聚類
Seurat v3應用基于圖形的聚類方法，例如KNN方法。具有相似基因表達模式的細胞之間繪制邊緣，然后將他們劃分為一個內聯群體。
在PhenoGraph中，首先基于pca維度中（先前計算的pca數據）計算歐式距離（the euclidean distance），然后根據兩個細胞在局部的重合情況（Jaccard 相似系數）優(yōu)化兩個細胞之間的邊緣權值。此步驟內置于FindNeighbors函數，輸入時先前確定的pc數據。
為了聚類細胞，接下來應用模塊化優(yōu)化技術迭代將細胞聚集在一起。（the Louvain algorithm (default) or SLM [SLM, Blondel et al., Journal of Statistical Mechanics]），FindClusters函數實現這一功能，其中需要注意resolution參數，該參數設置下游聚類分析的“granularity”，更大的resolution會導致跟多的細胞類群。3000左右的細胞量，設置resolution為0.4-1.2是比較合適的。細胞數據集越大，需要更大的resolution參數, 會獲得更多的細胞聚類。
查看細胞屬于那個類群可以使用函數Idents。
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: 96033
Running Louvain algorithm...
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Maximum modularity in 10 random starts: 0.8720
Number of communities: 9
Elapsed time: 0 seconds
#查看細胞屬于那個類群
head(Idents(pbmc), 5)
AAACATACAACCAC-1 AAACATTGAGCTAC-1 AAACATTGATCAGC-1 AAACCGTGCTTCCG-1 AAACCGTGTATGCG-1
0 3 2 5 6
Levels: 0 1 2 3 4 5 6 7 8
#7 非線性降維分析
Seurat提供了一些非線性降維度分析的方法，如tSNE和UMAP，以可視化和探索這些數據集；這一步使用的數據建議與聚類分析使用的pc維度一致。
# install UMAP： reticulate::py_install(packages ='umap-learn')
pbmc <- RunUMAP(pbmc, dims = 1:10)
#畫圖展示
DimPlot(pbmc, reduction = "umap")
#添加細胞類群xiba的標簽
DimPlot(pbmc, reduction = "umap",label = TRUE)
LabelClusters(DimPlot(pbmc, reduction = "umap"),id = 'ident')
此時可以保存數據，方便下次直接導入數據修改圖形。
saveRDS(pbmc, file = "../output/pbmc_tutorial.rds")
#8 尋找差異表達基因 (cluster biomarkers)
Seurat可以通過差異表達分析尋找不同細胞類群的標記基因。FindMarkers函數可以進行此操作，但是默認尋找單個類群(參數ident.1)與其他所有類群陽性和陰性標記基因。FindAllMarkers函數會自動尋找每個類群和其他每個類群之間的標記基因。
min.pct參數:設定在兩個細胞群中任何一個被檢測到的百分比，通過此設定不檢測很少表達基因來縮短程序運行時間。默認0.1
thresh.test參數：設定在兩個細胞群中基因差異表達量?？梢栽O置為0 ，程序運行時間會更長。
max.cells.per.ident參數：每個類群細胞抽樣設置；也可以縮短程序運行時間。
# find all markers of cluster 1
cluster1.markers <- FindMarkers(pbmc, ident.1 = 1, min.pct = 0.25)
head(cluster1.markers, n = 5)
p_val avg_logFC pct.1 pct.2 p_val_adj
IL32 1.894810e-92 0.8373872 0.948 0.464 2.598542e-88
LTB 7.953303e-89 0.8921170 0.981 0.642 1.090716e-84
CD3D 1.655937e-70 0.6436286 0.919 0.431 2.270951e-66
IL7R 3.688893e-68 0.8147082 0.747 0.325 5.058947e-64
LDHB 2.292819e-67 0.6253110 0.950 0.613 3.144372e-63
# 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_logFC pct.1 pct.2 p_val_adj
FCGR3A 7.583625e-209 2.963144 0.975 0.037 1.040018e-204
IFITM3 2.500844e-199 2.698187 0.975 0.046 3.429657e-195
CFD 1.763722e-195 2.362381 0.938 0.037 2.418768e-191
CD68 4.612171e-192 2.087366 0.926 0.036 6.325132e-188
RP11-290F20.3 1.846215e-188 1.886288 0.840 0.016 2.531900e-184
# 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_logFC)
p_val avg_logFC pct.1 pct.2 p_val_adj cluster gene
<dbl> <dbl> <dbl> <dbl> <dbl> <fct> <chr>
1 1.96e-107 0.730 0.901 0.594 2.69e-103 0 LDHB
2 1.61e- 82 0.922 0.436 0.11 2.20e- 78 0 CCR7
3 7.95e- 89 0.892 0.981 0.642 1.09e- 84 1 LTB
4 1.85e- 60 0.859 0.422 0.11 2.54e- 56 1 AQP3
5 0. 3.86 0.996 0.215 0. 2 S100A9
6 0. 3.80 0.975 0.121 0. 2 S100A8
7 0. 2.99 0.936 0.041 0. 3 CD79A
8 9.48e-271 2.49 0.622 0.022 1.30e-266 3 TCL1A
9 2.96e-189 2.12 0.985 0.24 4.06e-185 4 CCL5
10 2.57e-158 2.05 0.587 0.059 3.52e-154 4 GZMK
11 3.51e-184 2.30 0.975 0.134 4.82e-180 5 FCGR3A
12 2.03e-125 2.14 1 0.315 2.78e-121 5 LST1
13 7.95e-269 3.35 0.961 0.068 1.09e-264 6 GZMB
14 3.13e-191 3.69 0.961 0.131 4.30e-187 6 GNLY
15 1.48e-220 2.68 0.812 0.011 2.03e-216 7 FCER1A
16 1.67e- 21 1.99 1 0.513 2.28e- 17 7 HLA-DPB1
17 7.73e-200 5.02 1 0.01 1.06e-195 8 PF4
18 3.68e-110 5.94 1 0.024 5.05e-106 8 PPBP
Seurat可以通過參數test.use設定檢驗差異表達的方法(詳情見 DE vignett)。
cluster1.markers <- FindMarkers(pbmc, ident.1 = 0, logfc.threshold = 0.25, test.use = "roc", only.pos = TRUE)
head(cluster1.markers, n = 5)
Seurat有多種方法可視化標記基因的方法
- VlnPlot: 基于細胞類群的基因表達概率分布
- FeaturePlot：在tSNE 或 PCA圖中畫出基因表達情況
- 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為指定的細胞和基因花表達熱圖。每個類群默認展示top 20標記基因。
top10 <- pbmc.markers %>% group_by(cluster) %>% top_n(n = 10, wt = avg_logFC)
DoHeatmap(pbmc, features = top10$gene) + NoLegend()
#9 Assigning cell type identity to clusters
根據已知細胞類型與基因標記的對應關系，可以為細胞類群找到對應的細胞類型。
| Cluster ID | Markers | Cell Type |
|---|---|---|
| 0 | IL7R, CCR7 | Naive CD4+ T |
| 1 | IL7R, S100A4 | Memory CD4+ |
| 2 | CD14, LYZ | CD14+ Mono |
| 3 | MS4A1 | B |
| 4 | CD8A | CD8+ T |
| 5 | FCGR3A, MS4A7 | FCGR3A+ Mono |
| 6 | GNLY, NKG7 | NK |
| 7 | FCER1A, CST3 | DC |
| 8 | PPBP | Platelet |
new.cluster.ids <- c("Naive CD4 T", "Memory CD4 T", "CD14+ Mono", "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")
參考：
Seurat - Guided Clustering Tutorial, Compiled: April 17, 2020
Getting started with Seurat
