Engage to Life Energy
单细胞转录组基础分析二：数据质控与标准化
发布日期：2021-04-20浏览：
次
本专题将针对10X Genomics单细胞转录组数据演示各种主流分析，包括基于Seurat的基础分析、以及基于clusterProfiler、Monocle、SingleR等R包的延伸分析。
本节及之后的三、四节主要介绍单细胞表达矩阵到细胞类型鉴定的分析步骤，流程图如下：
10X官网演示数据：
https://support.10xgenomics.com/single-cell-gene-expression/datasets
官方演示数据集的页面如下，Chromium Connect是自动化系统，官方介绍操作造成的批次效应较小。Chromium Next GEM最新版的芯片，由老版的双十字微流控芯片改成了单十字微流控芯片。目前国内的10X试剂基本都是V3版本的了，因此也不要下载V1和V2试剂的数据了。为了保证笔记本电脑能运行，我们下载细胞数比较少的数据，下图红色箭头所示：
本教程的分析都是基于箭头标示的数据。点击之后需要提交个人信息，提交之后就进入数据下载界面了。
下载红框标示的数据，Seurat分析的输入文件在这里，解压后如下图所示：
运行上面的代码后会在"QC"文件夹下面得到4个文件
打vlnplot_before_qc的文件
运行上面的代码后会在"QC"文件夹下面得到vlnplot_after_qc的文件
标准化后基因表达量 = log1p（10000*基因counts/细胞总counts）
上游分析软件Cell Ranger会对数据进行质控，但是在进行下游分析前，我们一般会对数据进行更严格的过滤。
library(Seurat)library(tidyverse)rm(list=ls())dir.create("QC")##创建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) #查看样本的细胞数量
##计算质控指标#计算细胞中核糖体基因比例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)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") ggsave("QC/pearplot_before_qc.pdf", plot = pearplot, width = 12, height = 5) ggsave("QC/pearplot_before_qc.png", plot = pearplot, width = 12, height = 5)
上面的小提琴图依次是细胞的基因数量、mRNA分子数量、线粒体基因比例、红细胞基因比例。我们一般根据基因数量和线粒体比例来过滤细胞，细胞的最低基因数一般选择200-500，最大基因数和核糖体比例根据上图来选择，我的建议是minGene=500，maxGene=4000，pctMT=15。
##设置质控标准print(c("请输入允许基因数和核糖体比例，示例如下：", "minGene=500", "maxGene=4000", "pctMT=20"))minGene=500maxGene=4000pctMT=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)
##保存中间结果saveRDS(scRNA, file="scRNA.rds")
可以看到基因数和核糖体比例不正常的细胞都过滤了。以上代码的最后一行是对数据进行标准化，它的作用是让测序数据量不同的细胞的基因表达量具有可比性。计算公式如下：
