背景知識：

1. 線粒體基因表達(dá)比例：

   • 線粒體基因參與細(xì)胞凋亡的啟動，因此垂死細(xì)胞會體現(xiàn)出線粒體基因表達(dá)比例的上升。
   • 在我們使用組織塊制作單細(xì)胞懸液的過程中，會造成部分細(xì)胞細(xì)胞膜的破壞，這時(shí)細(xì)胞質(zhì)中的轉(zhuǎn)錄本會大量流失，而線粒體和細(xì)胞核中的轉(zhuǎn)錄本不便，檢測時(shí)就會表現(xiàn)出大量比例的線粒體基因表達(dá)。
   • 線粒體基因通常為MT開頭，我們在下面處理中會使用正則表達(dá)式以計(jì)算這一系列細(xì)胞的表達(dá)百分比作為質(zhì)控指標(biāo)。

2. 當(dāng)制備單細(xì)胞懸液時(shí)，如有兩個(gè)或多個(gè)細(xì)胞未解離的情況下，后續(xù)檢測結(jié)果中會出現(xiàn)某一細(xì)胞所檢測到的基因數(shù)和UMI數(shù)量的倍增。

3. 在某些數(shù)據(jù)中，還需根據(jù)檢測目標(biāo)和情況進(jìn)行其他指標(biāo)的篩選，例如紅細(xì)胞的剔除、細(xì)胞周期的選定等。。。

我們對數(shù)進(jìn)行質(zhì)控的過程就是對上述指標(biāo)進(jìn)行閾值的設(shè)定以剔除我們不需要的數(shù)據(jù)，避免對下游分析的影響。但，各項(xiàng)指標(biāo)閾值設(shè)定為多少？需要對那些指標(biāo)進(jìn)行篩選？是沒有一個(gè)統(tǒng)一的定論的。例如由于樣本來源的差異，必然會導(dǎo)致測得的基因數(shù)、UMI總數(shù)、線粒體比例的不同，這時(shí)我們就要大量閱讀所做樣本的相關(guān)文獻(xiàn)，找一個(gè)普遍性的閾值或根據(jù)對數(shù)據(jù)進(jìn)行可視化之后查看樣本分布情況，剔除離群樣本。

我們用pbmc3k數(shù)據(jù)集進(jìn)行代碼流程演示。強(qiáng)烈建議仔細(xì)閱讀seurat官網(wǎng)指南。

library(Seurat)
library(SeuratData)
library(ggplot2)
library(patchwork)
library(dplyr)
load(file = "pbmc.rdata") #Seurat對象讀取
pbmc
pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-") #計(jì)算線粒體基因百分比
#小提琴圖可視化質(zhì)控指標(biāo)
VlnPlot(pbmc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)

image

# 計(jì)算兩特征之間的關(guān)系
plot1 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "percent.mt") 
plot2 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "nFeature_RNA") 
plot1 + plot2

image

這里主要是計(jì)算了兩個(gè)指標(biāo)之間的相關(guān)性，個(gè)人理解count和feature應(yīng)該有較強(qiáng)的共線性，而mt比例應(yīng)無較強(qiáng)相關(guān)性。

#設(shè)定各指標(biāo)閾值使用subset函數(shù)取子集
pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)  #這里使用的標(biāo)準(zhǔn)為nFeature_RNA：200-2500，percent.mt < 5

到這里質(zhì)控的流程就基本結(jié)束了，結(jié)束后我們也可以重復(fù)上面的可視化過程來對質(zhì)控結(jié)果進(jìn)行驗(yàn)證。

Seurat降維標(biāo)準(zhǔn)流程

由于單細(xì)胞測序每個(gè)樣本都會存在數(shù)千個(gè)細(xì)胞，而每個(gè)細(xì)胞都會有最高上萬個(gè)所測基因表達(dá)量。
因此，為了后續(xù)分析的效率和可行性，我們對數(shù)據(jù)進(jìn)行降維的操作。先鑒定出細(xì)胞中的高變基因（在各個(gè)細(xì)胞中的擾動較大），僅僅對細(xì)胞分群維度有貢獻(xiàn)的高變基因進(jìn)行后續(xù)PCA，僅保留貢獻(xiàn)度較大的PCs用于細(xì)胞亞群的鑒定。
主要流程：尋找高變基因——PCA前的scale——PCA——可視化降維t-sne/umpa

個(gè)人理解：
尋找高變基因的意義：減輕PCA的運(yùn)算壓力并有利于PCA的降維。
PCA前的scale：是PCA之前的標(biāo)準(zhǔn)流程，線性轉(zhuǎn)換所有基因表達(dá)使細(xì)胞間的平均表達(dá)值為0。 縮放每個(gè)基因表達(dá)，使細(xì)胞間的方差為1。此步驟使每個(gè)基因在PCA中獲得同等權(quán)重，防止某一高表達(dá)基因占主導(dǎo)地位導(dǎo)致PC鑒定的偏差

pbmc <- NormalizeData(pbmc, normalization.method = "LogNormalize", scale.factor = 1e4)  #對數(shù)據(jù)進(jìn)行標(biāo)準(zhǔn)化
pbmc <- FindVariableFeatures(pbmc, selection.method = 'vst', nfeatures = 2000) #尋找高變基因
# 找出前10高可變基因用于后續(xù)可視化
top10 <- head(VariableFeatures(pbmc), 10)
# 高變基因可視化
plot1 <- VariableFeaturePlot(pbmc)
plot2 <- LabelPoints(plot = plot1, points = top10, repel = TRUE)
plot1 + plot2

image

#PCA前的scale和PCA
pbmc <- ScaleData(pbmc, vars.to.regress = "percent.mt")
pbmc <- RunPCA(pbmc, features = VariableFeatures(object = pbmc)) #默認(rèn)最大PC數(shù)為50，可查閱函數(shù)help自行修改參數(shù)

PCA流程結(jié)束。我們需要挑選合適的PCs進(jìn)行后續(xù)的細(xì)胞亞群鑒定
這里官方給出了兩種方法，JackStraw()和ElbowPlot()

pbmc <- JackStraw(pbmc, num.replicate = 100)#JackStraw方法
pbmc <- ScoreJackStraw(pbmc, dims = 1:20)
JackStrawPlot(pbmc, dims = 1:15) #可視化前15個(gè)PC
ElbowPlot(pbmc)#肘型圖
##默認(rèn)都是計(jì)算50個(gè)PC可根據(jù)需求調(diào)整參數(shù)

image

image

這里推薦一個(gè)技能樹jimmuy寫的一個(gè)方法來確定最佳PC數(shù)量，關(guān)于詳細(xì)的解釋，這里不贅述。我們就根據(jù)官方給定的PC繼續(xù)下游處理

# 細(xì)胞聚類
pbmc <- FindNeighbors(pbmc, dims = 1:10) #這里dim填入選擇的PC數(shù)量
pbmc <- FindClusters(pbmc, resolution = 0.5)  #分群粒度根據(jù)總細(xì)胞數(shù)量選取，一般為0.4-1.2，粒度越大分群越多
head(Idents(pbmc), 5)
table(pbmc$seurat_clusters) 
##UMPA流程并進(jìn)行分群可視化
pbmc <- RunUMAP(pbmc, dims = 1:10)#與上面的PC要一致
DimPlot(pbmc, reduction = 'umap')

image

到這里我們的單細(xì)胞數(shù)據(jù)質(zhì)控，降維標(biāo)準(zhǔn)流程就結(jié)束了，后續(xù)就是堆細(xì)胞亞群的注釋。在后續(xù)的注釋過程中也會發(fā)現(xiàn)各種問題，這時(shí)候就需要對上面選擇的PC、分群粒度進(jìn)行調(diào)整來對亞群進(jìn)行細(xì)分或合并。

參考來源：
https://satijalab.org/seurat/articles/pbmc3k_tutorial.html
http://t.zoukankan.com/emanlee-p-14932294.html
http://www.itdecent.cn/p/a77b8e7dc76d
問題交流：
Email: xuran@hrbmu.edu.cn