单细胞转录组基础分析二：数据质控与标准化
2021-03-10 本文已影响0人
PhageNanoenzyme
本文是参考学习单细胞转录组基础分析二：数据质控与标准化的学习笔记。可能根据学习情况有所改动。
本节及之后的三、四节主要介绍单细胞表达矩阵到细胞类型鉴定的分析步骤，流程图如下：
图片
10X官网下载数据
官方演示数据集的页面如下，Chromium Connect是自动化系统，官方介绍操作造成的批次效应较小。Chromium Next GEM最新版的芯片，由老版的双十字微流控芯片改成了单十字微流控芯片。目前国内的10X试剂基本都是V3版本的了，因此也不要下载V1和V2试剂的数据了。为了保证笔记本电脑能运行，我们下载细胞数比较少的数据，下图红色箭头所示：
图片
本教程的分析都是基于箭头标示的数据。点击之后需要提交个人信息，提交之后就进入数据下载界面 了。
图片
下载红框标示的数据，Seurat分析的输入文件在这里，解压后如下图所示：
图片
数据质控与标准化
上游分析软件Cell Ranger会对数据进行质控，但是在进行下游分析前，我们一般会对数据进行更严格的过滤。
library(Seurat)
运行上面的代码后会在"QC"文件夹下面得到4个文件
图片
打vlnplot_before_qc的文件
图片
上面的小提琴图依次是细胞的基因数量、mRNA分子数量、线粒体基因比例、红细胞基因比例。我们一般根据基因数量和线粒体比例来过滤细胞，细胞的最低基因数一般选择200-500，最大基因数和核糖体比例根据上图来选择，我的建议是minGene=500，maxGene=4000，pctMT=15。
##设置质控标准
运行上面的代码后会在"QC"文件夹下面得到vlnplot_after_qc的文件
图片
可以看到基因数和核糖体比例不正常的细胞都过滤了。以上代码的最后一行是对数据进行标准化，它的作用是让测序数据量不同的细胞的基因表达量具有可比性。计算公式如下：标准化后基因表达量 = log1p（10000基因counts/细胞总counts）*
> ##创建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) #查看样本的细胞数量
SeuratProject
1222
> ##计算质控指标
> #计算细胞中核糖体基因比例
> scRNA[["percent.mt"]] <- PercentageFeatureSet(scRNA, pattern = "^MT-")
> #计算红细胞比例
> 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, #不需要显示点，可以设置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)
> # 运行上面的代码后会在"QC"文件夹下面得到4个文件
> # 图片
> #
> # 打vlnplot_before_qc的文件
> # 图片
> #
> # 上面的小提琴图依次是细胞的基因数量、mRNA分子数量、线粒体基因比例、红细胞基因比例。我们一般根据基因数量和线粒体比例来过滤细胞，细胞的最低基因数一般选择200-500，最大基因数和核糖体比例根据上图来选择，我的建议是minGene=500，maxGene=4000，pctMT=15。
> ##设置质控标准
> print(c("请输入允许基因数和核糖体比例，示例如下：", "minGene=500", "maxGene=4000", "pctMT=20"))
[1] "请输入允许基因数和核糖体比例，示例如下：" "minGene=500"
[3] "maxGene=4000" "pctMT=20"
> minGene=500
> maxGene=4000
> pctMT=15
> ##数据质控
> 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%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
> ##保存中间结果
> saveRDS(scRNA, file="scRNA.rds")
