2021-07-02
上一篇(pian)10x單細胞(bao)數據分析我們介紹(shao)了如何(he)對(dui)單細胞(bao)進(jin)行表達(da)量分析，這一篇(pian)我們介紹(shao)一下單細胞(bao)分析流行的(de)R包Seurat的(de)分析步(bu)驟。
在開始分析之前，首先需要安裝最新版的R（//cran.r-project.org/bin/windows/base/）和Rstudio（//www.rstudio.com/products/rstudio/download/#download）。
打開Rstudio，需要安裝Seurat包(bao)，在控制臺中輸入：
install.packages("Seurat")安裝(zhuang)好(hao)以后，先將需要(yao)的包進行(xing)加(jia)載(zai)：
library(Seurat) library(dplyr) library(magrittr)
測試數據選擇10X Genomics 免費提供的外周血單核細胞 (PBMC) 數據集，約有2700個細胞（//cf.10xgenomics.com/samples/cell/pbmc3k/pbmc3k_filtered_gene_bc_matrices.tar.gz）。下載好以后需要進行解壓。從壓縮包的名字中我們可以看到，其實就是上一篇cellranger分析的結果中的filtered_gene_bc_matrices文件夾。測試數據的細胞數較少，普通的臺式機就可以分析。小編在測試的時候，分析完大約花了10分鐘，只用1G左右的內存，所以大家也可以在自己電腦上試著分析一下哦！
1、數據(ju)讀(du)入
Seurat需要(yao)的輸入信息為(wei)表達量矩(ju)(ju)陣，矩(ju)(ju)陣行為(wei)基因，列為(wei)細胞。使(shi)用Seurat提供的Read10X函數可以很(hen)方(fang)便的將(jiang)10x結(jie)果讀入到(dao)R矩(ju)(ju)陣中。使(shi)用CreateSeuratObject生成Seurat對(dui)象，后(hou)續分析都是在該(gai)對(dui)象上(shang)進行操作。
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)這(zhe)里的min.cells=3表(biao)示從表(biao)達(da)(da)量(liang)矩陣中刪除(chu)在少于3個(ge)細(xi)胞(bao)中表(biao)達(da)(da)的基因（行），min.features=200表(biao)示刪除(chu)表(biao)達(da)(da)的基因數小于200的細(xi)胞(bao)（列(lie)）。
2、質控和過濾
質(zhi)控的目的是去(qu)除掉低(di)質(zhi)量的數據，包括破損(sun)或死(si)亡的細胞(bao)、沒捕獲(huo)到細胞(bao)的empty droplet和捕獲(huo)到2個以上(shang)細胞(bao)的doublets。一般(ban)低(di)質(zhi)量的細胞(bao)或者empty droplet通常含有很少的基因(yin)，而doublets容易測到更多(duo)的基因(yin)。另一方面(mian)，低(di)質(zhi)量或者死(si)亡細胞(bao)會測到更多(duo)的線粒體基因(yin)表達的RNA。
使用PercentageFeatureSet函數評估每個細胞中的(de)線粒體表達比例：
pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")
pattern = "^MT-"表示線粒(li)(li)體基因的匹配(pei)模式，人的線粒(li)(li)體基因以(yi)MT-為(wei)前綴(zhui)(zhui)，小鼠(shu)的以(yi)Mt-為(wei)前綴(zhui)(zhui)。
Seurat對象(xiang)在生(sheng)成(cheng)時(shi)，已經將每個細胞表達的(de)基因(yin)數和reads數記錄(lu)在meta.data的(de)nFeature_RNA和nCount_RNA中。我們使用VlnPlot函數將nFeature_RNA、nCount_RNA和percent.mt三(san)個指(zhi)標進(jin)行可視(shi)化：
VlnPlot(pbmc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
根據小提琴圖的(de)數(shu)值分布情(qing)況確定篩(shai)選的(de)閾值。這里，我們以 200 < nFeature_RNA < 2500 且percent.mt < 5為(wei)條件對數(shu)據進行篩(shai)選：
pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)
3、標準(zhun)化(hua)
基因的(de)表達(da)矩陣需要經(jing)過標準化后才能進(jin)行后續的(de)分析。使用NormalizeData函數進(jin)行標準化，方法選擇(ze)LogNormalize：
pbmc <- NormalizeData(pbmc, normalization.method = "LogNormalize", scale.factor = 10000)
4、挑(tiao)選高變基因
接(jie)下來需要(yao)挑選(xuan)出在細(xi)胞間高(gao)變(bian)異的(de)基因(yin)，這些基因(yin)有(you)助于突出單細(xi)胞數據中的(de)生物學信息(xi)。使用FindVariableFeatures函數來挑選(xuan)高(gao)變(bian)基因(yin)，默(mo)認的(de)情況下，挑選(xuan)2000個。
pbmc <- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 2000)
5、數據降維（PCA）
在PCA之(zhi)前(qian)，需(xu)要講(jiang)數據(ju)進行歸一化(hua)：
all.genes <- rownames(pbmc) pbmc <- ScaleData(pbmc, features = all.genes)
歸一化(hua)后的數據進行PCA分(fen)析(xi)。默認的，只用前面挑選(xuan)的高變基因進行PCA分(fen)析(xi)，也可以(yi)使用features參數設置自己選(xuan)擇的基因進行PCA分(fen)析(xi)。
pbmc <- RunPCA(pbmc, features = VariableFeatures(object = pbmc))
6、聚類分析
在(zai)PCA分(fen)析(xi)后，就要根(gen)據PCA結(jie)果來(lai)進行(xing)(xing)聚(ju)類分(fen)析(xi)了。但(dan)是(shi)(shi)在(zai)進行(xing)(xing)聚(ju)類分(fen)析(xi)之前，需要選擇使用對(dui)少個(ge)主(zhu)成(cheng)分(fen)進行(xing)(xing)計算。每個(ge)主(zhu)成(cheng)分(fen)實際上代表的(de)是(shi)(shi)相關基因集的(de)信息，因此(ci)確定(ding)多少個(ge)主(zhu)成(cheng)分(fen)是(shi)(shi)一個(ge)重要的(de)步驟，我們(men)可以根(gen)據ElbowPlot函數(shu)來(lai)判(pan)斷。
pbmc <- JackStraw(pbmc, num.replicate = 100) pbmc <- ScoreJackStraw(pbmc, dims = 1:20) ElbowPlot(pbmc)
通常我(wo)們會將圖中(zhong)拐點所對應的值作為(wei)篩選的閾(yu)值，本例中(zhong)我(wo)們選擇(ze)前10個PC來進行聚類分析（通常使(shi)用更高(gao)的PC對結果的影響(xiang)并不(bu)大(da)，但是使(shi)用較小的PC會對結果產生不(bu)利(li)的影響(xiang)）。
Seurat基于 PCA 空間中的(de)歐(ou)幾里德距離構建一個(ge) KNN 圖(tu)，并(bing)根(gen)據其局部(bu)鄰(lin)域中的(de)共享重疊優化任意兩(liang)個(ge)單元之(zhi)間的(de)邊權重，基于Louvain算法(fa)對細(xi)胞進(jin)行聚(ju)類分群。
pbmc <- FindNeighbors(pbmc, dims = 1:10) pbmc <- FindClusters(pbmc, resolution = 0.5)
resolution表(biao)示(shi)分(fen)群的(de)分(fen)辨(bian)率，更高的(de)值會分(fen)出更多的(de)細胞群。
7、非(fei)線性降維
Seurat提(ti)供幾個非線性降(jiang)維方法，如(ru)UMAP和tSNE，來可視(shi)化和探索單細(xi)胞數據。這(zhe)些(xie)算法是通過學習數據的底層流形，來將(jiang)細(xi)胞放到二維空間(jian)中進(jin)行展示：
pbmc <- RunUMAP(pbmc, dims = 1:10) pbmc <- RunTSNE(pbmc, dims = 1:10)
這一步并(bing)不會影響細(xi)胞的(de)分(fen)(fen)群結(jie)果，所(suo)以(yi)也可(ke)以(yi)在分(fen)(fen)群前運行(xing)。可(ke)以(yi)通(tong)過(guo)DimPlot函數對非(fei)線性降(jiang)維(wei)的(de)結(jie)果進(jin)行(xing)展示：
umap.plot<-DimPlot(pbmc, reduction = "umap", label=T) tsne.plot<-DimPlot(pbmc, reduction = "tsne", label=T) umap.plot+tsne.plot
8、鑒定細胞群間Marker基因
為(wei)了鑒(jian)定各個細(xi)胞群對應的(de)細(xi)胞類型，我們(men)需要找(zhao)出(chu)每個細(xi)胞群相對高表達(da)的(de)基(ji)(ji)因。使用Seurat的(de)FindAllMarkers函(han)數分別計算每個亞群與(yu)剩下的(de)所(suo)有細(xi)胞的(de)差異性，通常只保留(liu)高表達(da)的(de)基(ji)(ji)因作為(wei)marker：
pbmc.markers <- FindAllMarkers(pbmc, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)
當然，對于PBMC這(zhe)種已(yi)知細胞(bao)類型marker的(de)數據，也可(ke)以直(zhi)接(jie)使用已(yi)知的(de)marker。以Seurat官方提供的(de)PBMC marker為例：
對(dui)這些marker進行可視化：
VlnPlot(pbmc, features = c("MS4A1", "CD79A"))
FeaturePlot(pbmc, features = c("MS4A1", "GNLY", "CD3E", "CD14", "FCER1A", "FCGR3A", "LYZ", "PPBP", "CD8A"))
根(gen)據marker基因的展示結果(guo)可以對細胞群賦予(yu)正式的細胞類型
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等，后續的單細胞分析請關注派森諾生物官網
