本文是參考學(xué)習(xí)單細(xì)胞轉(zhuǎn)錄組基礎(chǔ)分析二：數(shù)據(jù)質(zhì)控與標(biāo)準(zhǔn)化的學(xué)習(xí)筆記?？赡芨鶕?jù)學(xué)習(xí)情況有所改動。
本節(jié)及之后的三、四節(jié)主要介紹單細(xì)胞表達(dá)矩陣到細(xì)胞類型鑒定的分析步驟，流程圖如下：

圖片

10X官網(wǎng)下載數(shù)據(jù)

10X官網(wǎng)演示數(shù)據(jù)：

官方演示數(shù)據(jù)集的頁面如下，Chromium Connect是自動化系統(tǒng)，官方介紹操作造成的批次效應(yīng)較小。Chromium Next GEM最新版的芯片，由老版的雙十字微流控芯片改成了單十字微流控芯片。目前國內(nèi)的10X試劑基本都是V3版本的了，因此也不要下載V1和V2試劑的數(shù)據(jù)了。為了保證筆記本電腦能運(yùn)行，我們下載細(xì)胞數(shù)比較少的數(shù)據(jù)，下圖紅色箭頭所示：

圖片

本教程的分析都是基于箭頭標(biāo)示的數(shù)據(jù)。點擊之后需要提交個人信息，提交之后就進(jìn)入數(shù)據(jù)下載界面 了。

圖片

下載紅框標(biāo)示的數(shù)據(jù)，Seurat分析的輸入文件在這里，解壓后如下圖所示：

圖片

數(shù)據(jù)質(zhì)控與標(biāo)準(zhǔn)化

上游分析軟件Cell Ranger會對數(shù)據(jù)進(jìn)行質(zhì)控，但是在進(jìn)行下游分析前，我們一般會對數(shù)據(jù)進(jìn)行更嚴(yán)格的過濾。

library(Seurat)

運(yùn)行上面的代碼后會在"QC"文件夾下面得到4個文件

圖片

打vlnplot_before_qc的文件

圖片

上面的小提琴圖依次是細(xì)胞的基因數(shù)量、mRNA分子數(shù)量、線粒體基因比例、紅細(xì)胞基因比例。我們一般根據(jù)基因數(shù)量和線粒體比例來過濾細(xì)胞，細(xì)胞的最低基因數(shù)一般選擇200-500，最大基因數(shù)和核糖體比例根據(jù)上圖來選擇，我的建議是minGene=500，maxGene=4000，pctMT=15。

##設(shè)置質(zhì)控標(biāo)準(zhǔn)

運(yùn)行上面的代碼后會在"QC"文件夾下面得到vlnplot_after_qc的文件

圖片

可以看到基因數(shù)和核糖體比例不正常的細(xì)胞都過濾了。以上代碼的最后一行是對數(shù)據(jù)進(jìn)行標(biāo)準(zhǔn)化，它的作用是讓測序數(shù)據(jù)量不同的細(xì)胞的基因表達(dá)量具有可比性。計算公式如下：標(biāo)準(zhǔn)化后基因表達(dá)量 = log1p（10000基因counts/細(xì)胞總counts）*

> ##創(chuàng)建seurat對象
> scRNA.counts <- Read10X(data.dir = "filtered_feature_bc_matrix")
> try({scRNA = CreateSeuratObject(scRNA.counts[['Gene Expression']])},silent=T)
> if(exists('scRNA')){} else {scRNA = CreateSeuratObject(scRNA.counts)}
> table(scRNA@meta.data$orig.ident)         #查看樣本的細(xì)胞數(shù)量

SeuratProject 
         1222 
> ##計算質(zhì)控指標(biāo)
> #計算細(xì)胞中核糖體基因比例
> scRNA[["percent.mt"]] <- PercentageFeatureSet(scRNA, pattern = "^MT-")
> #計算紅細(xì)胞比例
> HB.genes <- c("HBA1","HBA2","HBB","HBD","HBE1","HBG1","HBG2","HBM","HBQ1","HBZ")
> HB_m <- match(HB.genes, rownames(scRNA@assays$RNA))
> HB.genes <- rownames(scRNA@assays$RNA)[HB_m]
> HB.genes <- HB.genes[!is.na(HB.genes)]
> scRNA[["percent.HB"]]<-PercentageFeatureSet(scRNA, features=HB.genes)
> head(scRNA@meta.data)
                      orig.ident nCount_RNA nFeature_RNA percent.mt percent.HB
AAACCCAAGGAGAGTA-1 SeuratProject       8288         2620  10.774614          0
AAACGCTTCAGCCCAG-1 SeuratProject       5512         1808   7.964441          0
AAAGAACAGACGACTG-1 SeuratProject       4283         1562   6.187252          0
AAAGAACCAATGGCAG-1 SeuratProject       2754         1225   5.991285          0
AAAGAACGTCTGCAAT-1 SeuratProject       6592         1831   6.614078          0
AAAGGATAGTAGACAT-1 SeuratProject       8845         2048   7.959299          0
> col.num <- length(levels(scRNA@active.ident))
> violin <- VlnPlot(scRNA,
+                   features = c("nFeature_RNA", "nCount_RNA", "percent.mt","percent.HB"), 
+                   cols =rainbow(col.num), 
+                   pt.size = 0.01, #不需要顯示點，可以設(shè)置pt.size = 0
+                   ncol = 4) + 
+   theme(axis.title.x=element_blank(), axis.text.x=element_blank(), axis.ticks.x=element_blank())
> ggsave("QC/vlnplot_before_qc.pdf", plot = violin, width = 12, height = 6)
> ggsave("QC/vlnplot_before_qc.png", plot = violin, width = 12, height = 6)
> plot1 <- FeatureScatter(scRNA, feature1 = "nCount_RNA", feature2 = "percent.mt")
> plot2 <- FeatureScatter(scRNA, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
> plot3 <- FeatureScatter(scRNA, feature1 = "nCount_RNA", feature2 = "percent.HB")
> pearplot <- CombinePlots(plots = list(plot1, plot2, plot3), nrow=1, legend="none")
Warning message:
CombinePlots is being deprecated. Plots should now be combined using the patchwork system. 
> ggsave("QC/pearplot_before_qc.pdf", plot = pearplot, width = 12, height = 5)
> ggsave("QC/pearplot_before_qc.png", plot = pearplot, width = 12, height = 5)
> # 運(yùn)行上面的代碼后會在"QC"文件夾下面得到4個文件
> # 圖片
> # 
> # 打vlnplot_before_qc的文件
> # 圖片
> # 
> # 上面的小提琴圖依次是細(xì)胞的基因數(shù)量、mRNA分子數(shù)量、線粒體基因比例、紅細(xì)胞基因比例。我們一般根據(jù)基因數(shù)量和線粒體比例來過濾細(xì)胞，細(xì)胞的最低基因數(shù)一般選擇200-500，最大基因數(shù)和核糖體比例根據(jù)上圖來選擇，我的建議是minGene=500，maxGene=4000，pctMT=15。
> ##設(shè)置質(zhì)控標(biāo)準(zhǔn)
> print(c("請輸入允許基因數(shù)和核糖體比例，示例如下：", "minGene=500", "maxGene=4000", "pctMT=20"))
[1] "請輸入允許基因數(shù)和核糖體比例，示例如下：" "minGene=500"                             
[3] "maxGene=4000"                             "pctMT=20"                                
> minGene=500
> maxGene=4000
> pctMT=15
> ##數(shù)據(jù)質(zhì)控
> scRNA <- subset(scRNA, subset = nFeature_RNA > minGene & nFeature_RNA < maxGene & percent.mt < pctMT)
> col.num <- length(levels(scRNA@active.ident))
> violin <-VlnPlot(scRNA,
+                  features = c("nFeature_RNA", "nCount_RNA", "percent.mt","percent.HB"), 
+                  cols =rainbow(col.num), 
+                  pt.size = 0.1, 
+                  ncol = 4) + 
+   theme(axis.title.x=element_blank(), axis.text.x=element_blank(), axis.ticks.x=element_blank())
> ggsave("QC/vlnplot_after_qc.pdf", plot = violin, width = 12, height = 6)
> ggsave("QC/vlnplot_after_qc.png", plot = violin, width = 12, height = 6)
> scRNA <- NormalizeData(scRNA, normalization.method = "LogNormalize", scale.factor = 10000)
Performing log-normalization
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
> ##保存中間結(jié)果
> saveRDS(scRNA, file="scRNA.rds")