2021-07-02
上一篇10x單細胞數(shù)據(jù)分析我們介紹了如何對單細胞進行表達量分析，這一篇我們介紹一下單細胞分析流行的R包Seurat的分析步驟。
在開始分析之前，首先需要安裝最新版的R（https://cran.r-project.org/bin/windows/base/）和Rstudio（https://www.rstudio.com/products/rstudio/download/#download）。
打開Rstudio，需要安裝Seurat包，在控制臺中輸入：
install.packages("Seurat")安裝好以后，先將需要的包進行加載：
library(Seurat) library(dplyr) library(magrittr)
測試數(shù)據(jù)選擇10X Genomics 免費提供的外周血單核細胞 (PBMC) 數(shù)據(jù)集，約有2700個細胞（https://cf.10xgenomics.com/samples/cell/pbmc3k/pbmc3k_filtered_gene_bc_matrices.tar.gz）。下載好以后需要進行解壓。從壓縮包的名字中我們可以看到，其實就是上一篇cellranger分析的結果中的filtered_gene_bc_matrices文件夾。測試數(shù)據(jù)的細胞數(shù)較少，普通的臺式機就可以分析。小編在測試的時候，分析完大約花了10分鐘，只用1G左右的內(nèi)存，所以大家也可以在自己電腦上試著分析一下哦！
1、數(shù)據(jù)讀入
Seurat需要的輸入信息為表達量矩陣，矩陣行為基因，列為細胞。使用Seurat提供的Read10X函數(shù)可以很方便的將10x結果讀入到R矩陣中。使用CreateSeuratObject生成Seurat對象，后續(xù)分析都是在該對象上進行操作。
pbmc.data <- Read10X("C:/Users/my/Desktop/pbmc3k_filtered_gene_bc_matrices/filtered_gene_bc_matrices/hg19/")
pbmc <- CreateSeuratObject(counts = pbmc.data, project = "pbmc3k", min.cells = 3, min.features = 200)這里的min.cells=3表示從表達量矩陣中刪除在少于3個細胞中表達的基因（行），min.features=200表示刪除表達的基因數(shù)小于200的細胞（列）。
2、質(zhì)控和過濾
質(zhì)控的目的是去除掉低質(zhì)量的數(shù)據(jù)，包括破損或死亡的細胞、沒捕獲到細胞的empty droplet和捕獲到2個以上細胞的doublets。一般低質(zhì)量的細胞或者empty droplet通常含有很少的基因，而doublets容易測到更多的基因。另一方面，低質(zhì)量或者死亡細胞會測到更多的線粒體基因表達的RNA。
使用PercentageFeatureSet函數(shù)評估每個細胞中的線粒體表達比例：
pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")
pattern = "^MT-"表示線粒體基因的匹配模式，人的線粒體基因以MT-為前綴，小鼠的以Mt-為前綴。
Seurat對象在生成時，已經(jīng)將每個細胞表達的基因數(shù)和reads數(shù)記錄在meta.data的nFeature_RNA和nCount_RNA中。我們使用VlnPlot函數(shù)將nFeature_RNA、nCount_RNA和percent.mt三個指標進行可視化：
VlnPlot(pbmc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
根據(jù)小提琴圖的數(shù)值分布情況確定篩選的閾值。這里，我們以 200 < nFeature_RNA < 2500 且percent.mt < 5為條件對數(shù)據(jù)進行篩選：
pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)
3、標準化
基因的表達矩陣需要經(jīng)過標準化后才能進行后續(xù)的分析。使用NormalizeData函數(shù)進行標準化，方法選擇LogNormalize：
pbmc <- NormalizeData(pbmc, normalization.method = "LogNormalize", scale.factor = 10000)
4、挑選高變基因
接下來需要挑選出在細胞間高變異的基因，這些基因有助于突出單細胞數(shù)據(jù)中的生物學信息。使用FindVariableFeatures函數(shù)來挑選高變基因，默認的情況下，挑選2000個。
pbmc <- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 2000)
5、數(shù)據(jù)降維（PCA）
在PCA之前，需要講數(shù)據(jù)進行歸一化：
all.genes <- rownames(pbmc) pbmc <- ScaleData(pbmc, features = all.genes)
歸一化后的數(shù)據(jù)進行PCA分析。默認的，只用前面挑選的高變基因進行PCA分析，也可以使用features參數(shù)設置自己選擇的基因進行PCA分析。
pbmc <- RunPCA(pbmc, features = VariableFeatures(object = pbmc))
6、聚類分析
在PCA分析后，就要根據(jù)PCA結果來進行聚類分析了。但是在進行聚類分析之前，需要選擇使用對少個主成分進行計算。每個主成分實際上代表的是相關基因集的信息，因此確定多少個主成分是一個重要的步驟，我們可以根據(jù)ElbowPlot函數(shù)來判斷。
pbmc <- JackStraw(pbmc, num.replicate = 100) pbmc <- ScoreJackStraw(pbmc, dims = 1:20) ElbowPlot(pbmc)
通常我們會將圖中拐點所對應的值作為篩選的閾值，本例中我們選擇前10個PC來進行聚類分析（通常使用更高的PC對結果的影響并不大，但是使用較小的PC會對結果產(chǎn)生不利的影響）。
Seurat基于 PCA 空間中的歐幾里德距離構建一個 KNN 圖，并根據(jù)其局部鄰域中的共享重疊優(yōu)化任意兩個單元之間的邊權重，基于Louvain算法對細胞進行聚類分群。
pbmc <- FindNeighbors(pbmc, dims = 1:10) pbmc <- FindClusters(pbmc, resolution = 0.5)
resolution表示分群的分辨率，更高的值會分出更多的細胞群。
7、非線性降維
Seurat提供幾個非線性降維方法，如UMAP和tSNE，來可視化和探索單細胞數(shù)據(jù)。這些算法是通過學習數(shù)據(jù)的底層流形，來將細胞放到二維空間中進行展示：
pbmc <- RunUMAP(pbmc, dims = 1:10) pbmc <- RunTSNE(pbmc, dims = 1:10)
這一步并不會影響細胞的分群結果，所以也可以在分群前運行?？梢酝ㄟ^DimPlot函數(shù)對非線性降維的結果進行展示：
umap.plot<-DimPlot(pbmc, reduction = "umap", label=T) tsne.plot<-DimPlot(pbmc, reduction = "tsne", label=T) umap.plot+tsne.plot
8、鑒定細胞群間Marker基因
為了鑒定各個細胞群對應的細胞類型，我們需要找出每個細胞群相對高表達的基因。使用Seurat的FindAllMarkers函數(shù)分別計算每個亞群與剩下的所有細胞的差異性，通常只保留高表達的基因作為marker：
pbmc.markers <- FindAllMarkers(pbmc, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)
當然，對于PBMC這種已知細胞類型marker的數(shù)據(jù)，也可以直接使用已知的marker。以Seurat官方提供的PBMC marker為例：
對這些marker進行可視化：
VlnPlot(pbmc, features = c("MS4A1", "CD79A"))
FeaturePlot(pbmc, features = c("MS4A1", "GNLY", "CD3E", "CD14", "FCER1A", "FCGR3A", "LYZ", "PPBP", "CD8A"))
根據(jù)marker基因的展示結果可以對細胞群賦予正式的細胞類型
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()
當然，細胞類型鑒定也有一些軟件可以進行分析，如SingleR、Cell-ID等，后續(xù)的單細胞分析請關注派森諾生物官網(wǎng)
