后臺(tái)有讀者翻到了一年前發(fā)的文獻(xiàn)解讀，請(qǐng)教了一下文章的圖的做法。正好前段時(shí)間剛做過(guò)單細(xì)胞轉(zhuǎn)錄組分析，今天就給大家介紹一下常用工具Seurat的用法。
Seurat 4.0 使用指南
設(shè)置Seurat對(duì)象
示例數(shù)據(jù)
10X Genomics免費(fèi)提供的外周血單核細(xì)胞(PBMC)數(shù)據(jù)集。通過(guò)Illumina NextSeq 500測(cè)序的2700個(gè)單細(xì)胞。示例數(shù)據(jù)下載:https://cf.10xgenomics.com/samples/cell/pbmc3k/pbmc3k_filtered_gene_bc_matrices.tar.gz。(或者后臺(tái)回復(fù)seurat,領(lǐng)取完整代碼及示例數(shù)據(jù))
setwd(".../filtered_gene_bc_matrices/hg19") #設(shè)置工作環(huán)境到數(shù)據(jù)所在文件夾
#安裝和加載所需包
BiocManager::install("Seurat")
BiocManager::install("dplyr")
BiocManager::install("patchwork")
library(dplyr)
library(Seurat)
library(patchwork)
#導(dǎo)入示例數(shù)據(jù)
pbmc.data <- Read10X(data.dir = ".../filtered_gene_bc_matrices/hg19/")#自行填寫數(shù)據(jù)所在文件夾
#創(chuàng)建Seurat對(duì)象
pbmc <- CreateSeuratObject(counts = pbmc.data, project = "pbmc3k", min.cells = 3, min.features = 200)
#過(guò)濾檢測(cè)少于200個(gè)基因的細(xì)胞（min.features = 200）和少于3個(gè)細(xì)胞檢測(cè)出的基因（min.cells = 3）
pbmc
#參數(shù)解釋
CreateSeuratObject(
counts, #未標(biāo)準(zhǔn)化的數(shù)據(jù)，如原始計(jì)數(shù)或TPMs
project = "CreateSeuratObject",#設(shè)置Seurat對(duì)象的項(xiàng)目名稱
assay = "RNA", #與初始輸入數(shù)據(jù)對(duì)應(yīng)的分析名稱
names.field = 1,#對(duì)于每個(gè)cell的初始標(biāo)識(shí)類，從cell的名稱中選擇此字段。例如，如果cell在輸入矩陣 #中被命名為BARCODE_CLUSTER_CELLTYPE，則設(shè)置名稱。字段設(shè)置為3以將初始標(biāo)識(shí)設(shè)置為 #CELLTYPE。
names.delim = "_", #對(duì)于每個(gè)cell的初始標(biāo)識(shí)類，從cell的列名中選擇此分隔符。例如，如果cell命名 #為bar - cluster - celltype，則將此設(shè)置為“-”，以便將cell名稱分離到其組成部分 #中，以選擇相關(guān)字段。
meta.data = NULL, #要添加到Seurat對(duì)象的其他單元級(jí)元數(shù)據(jù)。應(yīng)該是data.frame，其中行是單元格名稱，列 #是附加的元數(shù)據(jù)字段。
...
min.cells #包含至少在這些細(xì)胞檢測(cè)到的features。
min.features #包含至少檢測(cè)到這些features的細(xì)胞
)
> pbmc
An object of class Seurat
13714 features across 2700 samples within 1 assay
Active assay: RNA (13714 features, 0 variable features)
#1個(gè)數(shù)據(jù)集，包含2700個(gè)細(xì)胞，13714個(gè)基因。
數(shù)據(jù)矩陣
# 查看這三個(gè)基因的前三十行矩陣
pbmc.data[c("CD3D", "TCL1A", "MS4A1"), 1:30]
> pbmc.data[c("CD3D", "TCL1A", "MS4A1"), 1:30]
3 x 30 sparse Matrix of class "dgCMatrix"
[[ suppressing 30 column names ‘AAACATACAACCAC-1’, ‘AAACATTGAGCTAC-1’, ‘AAACATTGATCAGC-1’ ... ]]
CD3D 4 . 10 . . 1 2 3 1 . . 2 7 1 . . 1 3 . 2
TCL1A . . . . . . . . 1 . . . . . . . . . . .
MS4A1 . 6 . . . . . . 1 1 1 . . . . . . . . .
CD3D 3 . . . . . 3 4 1 5
TCL1A . 1 . . . . . . . .
MS4A1 36 1 2 . . 2 . . . .
.在矩陣中的值表示0(未檢測(cè)到分子)。由于scRNA-seq矩陣中的大多數(shù)值為0,因此Seurat在任何可能的情況下都使用稀疏矩陣表示。這為Drop-seq/inDrop/10x數(shù)據(jù)節(jié)省了大量?jī)?nèi)存和速度。
標(biāo)準(zhǔn)的預(yù)處理流程
下面的步驟包含了Seurat中的scRNA-seq數(shù)據(jù)的標(biāo)準(zhǔn)預(yù)處理流程。包括QC（質(zhì)控）、數(shù)據(jù)歸一化以及細(xì)胞的選擇和過(guò)濾。
QC和細(xì)胞篩選
常用的質(zhì)控指標(biāo)：
- 每個(gè)細(xì)胞在檢測(cè)到的特異基因數(shù)
- 低質(zhì)量細(xì)胞或空液滴通常只能檢測(cè)到非常少的基因
- 兩個(gè)或多個(gè)細(xì)胞被同時(shí)捕獲通常會(huì)有很高的基因數(shù)
- 每個(gè)細(xì)胞檢測(cè)到的分子總數(shù)（與基因密切相關(guān)）
- 每個(gè)細(xì)胞的線粒體基因比例
- 低質(zhì)量/瀕死細(xì)胞常表現(xiàn)出廣泛的線粒體污染
- 使用
PercentageFeatureSet()函數(shù)計(jì)算線粒體QC指標(biāo) - 使用所有以MT-開頭的基因作為一組線粒體基因
#向pbmc新增一列percent.mt數(shù)據(jù)
pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")
QC指標(biāo)儲(chǔ)存在哪？
每個(gè)細(xì)胞基因數(shù)和總分子數(shù)在建立seurat對(duì)象時(shí)就已經(jīng)自動(dòng)計(jì)算好了。
#展示前5個(gè)細(xì)胞的QC指標(biāo)
head(pbmc@meta.data, 5)
> head(pbmc@meta.data, 5)
orig.ident nCount_RNA
AAACATACAACCAC-1 pbmc3k 2419
AAACATTGAGCTAC-1 pbmc3k 4903
AAACATTGATCAGC-1 pbmc3k 3147
AAACCGTGCTTCCG-1 pbmc3k 2639
AAACCGTGTATGCG-1 pbmc3k 980
nFeature_RNA percent.mt
AAACATACAACCAC-1 779 3.0177759
AAACATTGAGCTAC-1 1352 3.7935958
AAACATTGATCAGC-1 1129 0.8897363
AAACCGTGCTTCCG-1 960 1.7430845
AAACCGTGTATGCG-1 521 1.2244898
#使用小提琴圖可視化QC指標(biāo)
VlnPlot(pbmc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
-
nFeature_RNA代表每個(gè)細(xì)胞測(cè)到的基因數(shù)目。 -
nCount_RNA代表每個(gè)細(xì)胞測(cè)到所有基因的表達(dá)量之和。 -
percent.mt代表測(cè)到的線粒體基因的比例。
#FeatureScatter通常用于可視化 feature-feature 相關(guān)性，
#nCount_RNA 與percent.mt的相關(guān)性
plot1 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "percent.mt")
#nCount_RNA與nFeature_RNA的相關(guān)性
plot2 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
plot1 + plot2 #合并兩圖
過(guò)濾線粒體基因表達(dá)比例過(guò)高的細(xì)胞，和一些極值細(xì)胞（可以根據(jù)小提琴圖判斷，查看兩端離群值）。
pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)
#濾掉 2500 > nFeature_RNA >200 和percent.mt < 5的數(shù)據(jù)
數(shù)據(jù)標(biāo)準(zhǔn)化
默認(rèn)情況下，使用全局縮放歸一化方法“LogNormalize”，用總表達(dá)量對(duì)每個(gè)細(xì)胞的基因表達(dá)式進(jìn)行歸一化，再乘以一個(gè)縮放因子(默認(rèn)為10,000)，然后對(duì)結(jié)果進(jìn)行log轉(zhuǎn)換。標(biāo)準(zhǔn)化的數(shù)值存儲(chǔ)在pbmc[["RNA"]]@data中。
pbmc <- NormalizeData(pbmc, normalization.method = "LogNormalize", scale.factor = 10000)
若所有調(diào)用的參數(shù)都是默認(rèn)值，則可省去。
pbmc <- NormalizeData(pbmc)
鑒定高變基因
接下來(lái)，計(jì)算數(shù)據(jù)集中表現(xiàn)出高細(xì)胞間變異的特征基因(即，它們?cè)谀承┘?xì)胞中高表達(dá)，而在其他細(xì)胞中低表達(dá))。這些基因有助于突出單細(xì)胞數(shù)據(jù)集中的生物信號(hào)。
用FindVariableFeatures()函數(shù)實(shí)現(xiàn)。默認(rèn)情況下，每個(gè)數(shù)據(jù)集返回2000個(gè)features 。這些將用于下游分析，如PCA。
pbmc <- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 2000)
# 查看最高變的10個(gè)基因
top10 <- head(VariableFeatures(pbmc), 10)
# 畫出不帶標(biāo)簽或帶標(biāo)簽基因點(diǎn)圖
plot1 <- VariableFeaturePlot(pbmc)
plot2 <- LabelPoints(plot = plot1, points = top10, repel = TRUE)
plot1 + plot2
數(shù)據(jù)縮放
線性變換(“縮放”)，是在PCA降維之前的一個(gè)標(biāo)準(zhǔn)預(yù)處理步驟。ScaleData()函數(shù)功能:
轉(zhuǎn)換每個(gè)基因的表達(dá)值，使每個(gè)細(xì)胞的平均表達(dá)為0
-
轉(zhuǎn)換每個(gè)基因的表達(dá)值，使細(xì)胞間的方差為1
- 此步驟在下游分析中具有相同的權(quán)重，因此高表達(dá)的基因不會(huì)占主導(dǎo)地位
結(jié)果存儲(chǔ)在
pbmc[["RNA"]]@scale.data中
all.genes <- rownames(pbmc)
pbmc <- ScaleData(pbmc, features = all.genes)
線性降維
接下來(lái)，對(duì)縮放的數(shù)據(jù)執(zhí)行PCA。默認(rèn)情況下，只使用前面確定的變量特性作為輸入，但是如果想選擇不同的子集，可以使用features參數(shù)來(lái)定義。
pbmc <- RunPCA(pbmc, features = VariableFeatures(object = pbmc))
Seurat提供了幾種有用的方法來(lái)可視化細(xì)胞和定義PCA的特性，包括VizDimReduction、DimPlot和DimHeatmap
#查看PCA結(jié)果
print(pbmc[["pca"]], dims = 1:5, nfeatures = 5)
> 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()可以方便地探索數(shù)據(jù)集中異質(zhì)性的主要來(lái)源，并且可以確定哪些PC維度可以用于下一步的下游分析。細(xì)胞和基因根據(jù)PCA分?jǐn)?shù)來(lái)排序。
DimHeatmap(pbmc, dims = 1, cells = 500, balanced = TRUE) #1個(gè)PC 500個(gè)細(xì)胞
DimHeatmap(pbmc, dims = 1:15, cells = 500, balanced = TRUE) #15個(gè)PC
確定數(shù)據(jù)的維度
主成分分析的原理非常簡(jiǎn)單，概括來(lái)說(shuō)就是選擇包含信息量大的維度（features），去除信息量少的“干擾”維度。所以這里會(huì)有個(gè)問(wèn)題——如何知道保留幾個(gè)維度是最佳的呢？我們希望通過(guò)保留盡可能少的維度來(lái)留存盡可能多的信息。Seurat有兩種方法來(lái)確定維度。
- JackStraw
pbmc <- JackStraw(pbmc, num.replicate = 100)
pbmc <- ScoreJackStraw(pbmc, dims = 1:20)
JackStrawPlot(pbmc, dims = 1:15)
可以看出，在10-12個(gè)PC之后，顯著性大幅下降，也就是前10-12個(gè)維度包含了大部分的樣本信息。
- Elbow plot
ElbowPlot(pbmc)
可以看出，PC9-10附近有一個(gè)拐點(diǎn)（“elbow”），這表明大部分真實(shí)信號(hào)是在前10個(gè)pc中捕獲的。
綜合以上方法，選擇10個(gè)主成成分作為參數(shù)用于后續(xù)分析。
細(xì)胞聚類
Seurat使用KNN算法進(jìn)行聚類。
pbmc <- FindNeighbors(pbmc, dims = 1:10)
pbmc <- FindClusters(pbmc, resolution = 0.5)
#dims = 1:10 即選取前10個(gè)主成分來(lái)分類細(xì)胞。
#查看前5個(gè)細(xì)胞的分類ID
head(Idents(pbmc), 5)
非線性降維（UMAP/tSNE）
pbmc <- RunUMAP(pbmc, dims = 1:10)
DimPlot(pbmc, reduction = "umap")
# 顯示在聚類標(biāo)簽
DimPlot(pbmc, reduction = "umap", label = TRUE)
# 使用TSNE聚類
pbmc <- RunTSNE(pbmc, dims = 1:10)
DimPlot(pbmc, reduction = "tsne")
# 顯示在聚類標(biāo)簽
DimPlot(pbmc, reduction = "tsne", label = TRUE)
#保存rds，用于后續(xù)分析
saveRDS(pbmc, file = "../output/pbmc_tutorial.rds")
找差異表達(dá)基因(聚類標(biāo)志cluster biomarkers)
利用 FindMarkers 命令，可以找到找到各個(gè)細(xì)胞類型中與其他類別的差異表達(dá)基因，作為該細(xì)胞類型的生物學(xué)標(biāo)記基因。
dent.1參數(shù)設(shè)置待分析的細(xì)胞類別min.pct參數(shù)，在兩組細(xì)胞中的任何一組中檢測(cè)到的最小百分thresh.test參數(shù)，在兩組細(xì)胞間以一定數(shù)量的差異表達(dá)（平均）max.cells.per.ident參數(shù)，通過(guò)降低每個(gè)類的采樣值，提高計(jì)算速度
# cluster 1的標(biāo)記基因
cluster1.markers <- FindMarkers(pbmc, ident.1 = 1, min.pct = 0.25)
head(cluster1.markers, n = 5)
#找出區(qū)分cluster 5與cluster 0和cluster 3的所有標(biāo)記
cluster5.markers <- FindMarkers(pbmc, ident.1 = 5, ident.2 = c(0, 3), min.pct = 0.25)
head(cluster5.markers, n = 5)
# 找出每個(gè)cluster的標(biāo)記與所有剩余的細(xì)胞相比較，只報(bào)告陽(yáng)性細(xì)胞
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)
可視化
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"))
#每個(gè)聚類前10個(gè)差異基因表達(dá)熱圖(如果小于10，則繪制所有標(biāo)記)
top10 <- pbmc.markers %>% group_by(cluster) %>% top_n(n = 10, wt = avg_log2FC)
DoHeatmap(pbmc, features = top10$gene) + NoLegend()
鑒定細(xì)胞類型
這個(gè)數(shù)據(jù)集的markers與已知細(xì)胞的marker可以輕松配對(duì)。也可以通過(guò)查閱相關(guān)文獻(xiàn)人工注釋，或者利用singleR（挖個(gè)坑，有空再來(lái)填）自動(dòng)注釋。
| Cluster ID | Markers | Cell Type |
|---|---|---|
| 0 | IL7R, CCR7 | Naive CD4+ T |
| 1 | CD14, LYZ | CD14+ Mono |
| 2 | IL7R, S100A4 | Memory CD4+ |
| 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", "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")
往期內(nèi)容：
