背景知識：
-
線粒體基因表達(dá)比例：
- 線粒體基因參與細(xì)胞凋亡的啟動(dòng)火窒，因此垂死細(xì)胞會(huì)體現(xiàn)出線粒體基因表達(dá)比例的上升想虎。
- 在我們使用組織塊制作單細(xì)胞懸液的過程中神郊，會(huì)造成部分細(xì)胞細(xì)胞膜的破壞临谱，這時(shí)細(xì)胞質(zhì)中的轉(zhuǎn)錄本會(huì)大量流失跑筝，而線粒體和細(xì)胞核中的轉(zhuǎn)錄本不便碗誉，檢測時(shí)就會(huì)表現(xiàn)出大量比例的線粒體基因表達(dá)召嘶。
- 線粒體基因通常為
MT開頭，我們在下面處理中會(huì)使用正則表達(dá)式以計(jì)算這一系列細(xì)胞的表達(dá)百分比作為質(zhì)控指標(biāo)哮缺。
當(dāng)制備單細(xì)胞懸液時(shí)弄跌，如有兩個(gè)或多個(gè)細(xì)胞未解離的情況下，后續(xù)檢測結(jié)果中會(huì)出現(xiàn)某一細(xì)胞所檢測到的基因數(shù)和UMI數(shù)量的倍增尝苇。
在某些數(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)一的定論的。例如由于樣本來源的差異锤悄，必然會(huì)導(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)
nCount_RNA：每個(gè)細(xì)胞的UMI數(shù)量nFeature_RNA：基因數(shù)percent.mt：線粒體基因百分比在實(shí)際處理中沿腰，往往GEO存儲(chǔ)的都為質(zhì)控好的數(shù)據(jù)或在上游cellranger處理后進(jìn)行了過濾。
# 計(jì)算兩特征之間的關(guān)系
plot1 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "percent.mt")
plot2 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
plot1 + plot2
這里主要是計(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è)樣本都會(huì)存在數(shù)千個(gè)細(xì)胞，而每個(gè)細(xì)胞都會(huì)有最高上萬個(gè)所測基因表達(dá)量企巢。
因此枫慷，為了后續(xù)分析的效率和可行性，我們對數(shù)據(jù)進(jìn)行降維的操作浪规。先鑒定出細(xì)胞中的高變基因（在各個(gè)細(xì)胞中的擾動(dòng)較大）或听，僅僅對細(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
我們在做的時(shí)候可省略可視化步驟
#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ù)
JackStraw()主要是根據(jù)查看PC的P值來進(jìn)行篩選滔驶，最佳PC應(yīng)該為PC的p-value突然開始明顯變大的那個(gè)PC。個(gè)人體會(huì)卿闹，這種方法計(jì)算量巨大揭糕，如果數(shù)據(jù)集比較大計(jì)算時(shí)間很長，而且結(jié)果很難分辨锻霎，不做推薦
肘型圖的結(jié)果也比較難以分辨著角，按我們之前的認(rèn)識應(yīng)該選取拐點(diǎn)出，7或8旋恼，但官方在下游卻設(shè)定為了10吏口。
實(shí)際上該選擇多少并沒有一個(gè)明確的規(guī)定，往往只能通過繼續(xù)向下游分群注釋去做冰更，出現(xiàn)問題了回來調(diào)整或者往下做之后更換幾個(gè)PC再做一次看結(jié)果的重復(fù)性是否良好
這里推薦一個(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')
到這里我們的單細(xì)胞數(shù)據(jù)質(zhì)控奠衔，降維標(biāo)準(zhǔn)流程就結(jié)束了谆刨，后續(xù)就是堆細(xì)胞亞群的注釋塘娶。在后續(xù)的注釋過程中也會(huì)發(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.reibang.com/p/a77b8e7dc76d
問題交流：
Email: xuran@hrbmu.edu.cn
