质控过滤以及标准化

生信喵实验柴

发布于 2022-10-25 19:38:43

1.2K0

发布于 2022-10-25 19:38:43

文章被收录于专栏：生信喵实验柴

背景

标准处理流程：读取数据后对矩阵进行标准的处理流程，包括 QC 过滤，数据标准化以及检测差异表达的基因组。

一、数据质控

数据质控就是分别对行（基因）以及列（细胞）进行过滤，目的是过滤掉一些不符合要求的数据，例如油滴中的空细胞，多细胞，以及死细胞。

质控的参数主要有两个：

1、每个细胞测到的 UMI 总数，过低为空细胞，过高为多细胞；同时检测表达的基因数目，过高或过低都有问题；

2、检测每个细胞的线粒体基因的比例，理论上线粒体基因组与核基因组相比，只占很小一部分，当细胞裂解后线粒体比率会增高，cellranger 过滤中不包含该步骤，seurat 可以根据进行线粒体基因比率进行过滤。

#计算线粒体基因比例
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 代表测到的线粒体基因的比例。

#绘制散点图
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%。其他数据集请使用不同的标准。例如当测序量增加，细胞数增多时，选择相应的参数。

#过滤条件根据不同测序项目进行调整
dim(pbmc) # 13714  2700
pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 3500 & percent.mt < 5)
dim(pbmc) # 13714  2643 基因数没变，细胞减少了

对过滤后的数据集进行可视化。

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 )

pbmc <- NormalizeData(pbmc, normalization.method = "LogNormalize", scale.factor = 10000)
pbmc@commands #已经存储了标准化的数据
pbmc@commands$NormalizeData.RNA

三、识别高变（表达差异）基因

这一步的目的是鉴定出细胞与细胞之间表达量相差很大的基因，用于后续鉴定细胞类型，一些基因在部分细胞中高表达，部分细胞中低表达。研究表明，利用这些基因在下游分析中更容易找到关联。

单细胞测序中的表达差异基因与 bulk RNAseq 处理组对照组比较不同，而是一个 cluster 与其他 cluster 之间的差别。

#默认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