背景
标准处理流程:读取数据后对矩阵进行标准的处理流程,包括 QC 过滤,数据标准化以及检测差异表达的基因组。
一、数据质控
数据质控就是分别对行(基因)以及列(细胞)进行过滤,目的是过滤掉一些不符合要求的数据,例如油滴中的空细胞,多细胞,以及死细胞。
质控的参数主要有两个:
1、每个细胞测到的 UMI 总数,过低为空细胞,过高为多细胞;同时检测表达的基因数目,过高或过低都有问题;
2、检测每个细胞的线粒体基因的比例,理论上线粒体基因组与核基因组相比,只占很小一部分,当细胞裂解后线粒体比率会增高,cellranger 过滤中不包含该步骤,seurat 可以根据进行线粒体基因比率进行过滤。
代码语言:javascript复制#计算线粒体基因比例
pbmc[['percent.mt']] <- PercentageFeatureSet(pbmc,pattern = "^MT-")
pbmc[['percent.mt']]
#绘制小提琴图
VlnPlot(pbmc,features = c("nFeature_RNA","nCount_RNA","percent.mt"),ncol=3)
小提琴图
nFeature_RNA 代表每个细胞测到的基因数目,nCount 代表每个细胞测到所有基因的表达量之和,percent.mt 代表测到的线粒体基因的比例。
代码语言:javascript复制#绘制散点图
plot1 <- FeatureScatter(pbmc,feature1 = "nCount_RNA",feature2 = "percent.mt")
plot2 <- FeatureScatter(pbmc,feature1 = "nCount_RNA",feature2 = "nFeature_RNA")
plot1
plot2
plot1 plot2
散点图
去除线粒体基因表达比例过高的细胞,和一些极值细胞。
本案例中选择表达基因在 200-3500 之间,线粒体基因比例小于 5%。其他数据集请使用不同的标准。例如当测序量增加,细胞数增多时,选择相应的参数。
代码语言:javascript复制#过滤条件根据不同测序项目进行调整
dim(pbmc) # 13714 2700
pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 3500 & percent.mt < 5)
dim(pbmc) # 13714 2643 基因数没变,细胞减少了
对过滤后的数据集进行可视化。
代码语言:javascript复制plot1 <- FeatureScatter(pbmc, feature1 ="nCount_RNA", feature2 = "percent.mt") NoLegend()
plot2 <- FeatureScatter(pbmc, feature1 ="nCount_RNA", feature2 = "nFeature_RNA") NoLegend()
# plot1 plot2
CombinePlots(plots = list(plot1, plot2))
QC 之后散点图
二、标准化
由于不同细胞中测序到不同的数据,需要对表达量进行数据标准化,标准化有多种方式,默认使用 LogNormalize 的标准化算法,还有 CLR,RC 等。
A = log( 1 ( UMIA ÷ UMITotal ) × 10000 )
代码语言:javascript复制pbmc <- NormalizeData(pbmc, normalization.method = "LogNormalize", scale.factor = 10000)
pbmc@commands #已经存储了标准化的数据
pbmc@commands$NormalizeData.RNA
三、识别高变(表达差异)基因
这一步的目的是鉴定出细胞与细胞之间表达量相差很大的基因,用于后续鉴定细胞类型,一些基因在部分细胞中高表达,部分细胞中低表达。研究表明,利用这些基因在下游分析中更容易找到关联。
单细胞测序中的表达差异基因与 bulk RNAseq 处理组对照组比较不同,而是一个 cluster 与其他 cluster 之间的差别。
代码语言:javascript复制#默认2000个,通过nfeatures调整
pbmc1 <- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 2000)
#前十高变基因
top10 <- head(VariableFeatures(pbmc1), 10)
top10
#可视化
plot1 <- VariableFeaturePlot(pbmc1) guides(color="none")
plot1
plot2 <- LabelPoints(plot = plot1,points = top10,repel = TRUE,xnudge = 0,ynudge = 0)
plot2
plot1 plot2
#默认2000个,通过nfeatures调整成22000超过所有基因13000多。
pbmc1 <- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 22000)
#前十高变基因
top10 <- head(VariableFeatures(pbmc1), 10)
top10
#可视化
plot1 <- VariableFeaturePlot(pbmc1) guides(color="none")
plot1
plot2 <- LabelPoints(plot = plot1,points = top10,repel = TRUE,xnudge = 0,ynudge = 0)
plot2
plot1 plot2
标记高表达基因
写在最后:有时间我们会努力更新的。大家互动交流可以前去论坛,地址在下面,复制去浏览器即可访问,弥补下公众号没有留言功能的缺憾。原地址暂未启用(bioinfoer.com)。
代码语言:javascript复制sx.voiceclouds.cn
有些板块也可以预设为大家日常趣事的分享等,欢迎大家来提建议。