当我们拿到表达矩阵后,需要对其进行质控(quality control
, QC
)去除质量较差的细胞
,降低噪音,而后再进行数据分析。😘
rm(list = ls())
library(tidyverse)
library(scater)
library(SingleCellExperiment)
library(AnnotationDbi)
library(org.Hs.eg.db)
library(EnsDb.Hsapiens.v86)
这里我们用一下之前介绍的counts
文件和annotation
文件,然后通过SingleCellExperiment
创建SingleCellExperiment
格式的文件。
counts <- read.delim("./molecules.txt", row.names = 1)
annotation <- read.delim("./annotation.txt", stringsAsFactors = T)
counts
annotation
# 注意assays必须是matrix
umi <- SingleCellExperiment(
assays = list(counts = as.matrix(counts)),
colData = annotation)
altExp(umi,"ERCC") <- umi[grep("^ERCC-",rownames(umi)), ]
umi <- umi[grep("^ERCC-",rownames(umi),invert = T), ]
umi
ensdb_genes <- genes(EnsDb.Hsapiens.v86)
MT_names <- ensdb_genes[seqnames(ensdb_genes) == "MT"]$gene_id
is_mito <- rownames(umi) %in% MT_names
table(is_mito)
这里我们用scater
包的perCellQCMetrics
和perFeatureQCMetrics
函数,分别对细胞和特征(基因
)进行分析。🤒
umi_cell <- perCellQCMetrics(umi,subsets=list(Mito=is_mito))
head(umi_cell)[1:4,1:4]
umi_feature <- perFeatureQCMetrics(umi)
head(umi_feature)
我们现在把上面的信息加到SingleCellExperiment
文件中去,即umi
。🥳
umi <- addPerCellQC(umi, subsets=list(Mito=is_mito))
umi <- addPerFeatureQC(umi)
当然你也可以自己设定cutoff
来过滤细胞或者基因。我们先来看一下分布情况。
hist(
umi$total,
breaks = 100
)
abline(v = 25000, col = "red")
hist(
umi_cell$detected,
breaks = 100
)
abline(v = 7000, col = "red")
有时很难找到一个明确的cutoff
,这个时候我们需要引入一个概念叫中位数绝对偏差(median absolute deviations
, MADs
),我们一般认为如果某个变量距离中位数多于3个MAD
,则认为该值是异常值
,需要进行剔除。✌️ 在实际分析中,如果我们发现检测到的基因数量很少,但线粒体基因比例高,一般认为是低质量细胞的标志。🤜🤛
1️⃣ sum
作为过滤条件。
qc.lib2 <- isOutlier(umi_cell$sum,
nmads = 3,
log=TRUE,
type="lower")
attr(qc.lib2, "thresholds")
2️⃣ detected
作为过滤条件。
qc.nexprs2 <- isOutlier(umi_cell$detected,
nmads = 3,
log=TRUE,
type="lower")
attr(qc.nexprs2, "thresholds")
3️⃣ ERCC
作为过滤条件。
qc.spike2 <- isOutlier(umi_cell$altexps_ERCC_percent,
nmads = 3,
type="higher")
attr(qc.spike2, "thresholds")
4️⃣ Mito_percent
作为过滤条件。
qc.mito2 <- isOutlier(umi_cell$subsets_Mito_percent,
nmads = 3,
type="higher")
attr(qc.mito2, "thresholds")
🫶 合并上述所有过滤条件。
discard2 <- qc.lib2 | qc.nexprs2 | qc.spike2 | qc.mito2
DataFrame(LibSize=sum(qc.lib2),
NExprs=sum(qc.nexprs2),
SpikeProp=sum(qc.spike2),
MitoProp=sum(qc.mito2),
Total=sum(discard2))
使用scater
包的quickPerCellQC
函数进行过滤。
reasons <- quickPerCellQC(umi_cell,
sub.fields = c("subsets_Mito_percent", "altexps_ERCC_percent"))
colSums(as.matrix(reasons))
我们把需要过滤的细胞信息加入到SingleCellExperiment
文件中,即umi
。👌
umi$discard <- reasons$discard
这里可以看到,线粒体基因比例高的细胞,一般counts都比较低。👀
plotColData(umi, x="sum", y="subsets_Mito_percent", colour_by="discard")
这里可以看到,线粒体基因比例高的细胞,一般检测到的基因都比较少。
plotColData(umi, x="sum", y="detected", colour_by="discard")
这里可以看到,线粒体基因比例高的细胞,一般检测到的ERCC比例都比较高。😂
plotColData(umi, x="altexps_ERCC_percent", y="subsets_Mito_percent",colour_by="discard")
这里我们将individual
设为分面参数,看一下不同invidual
的过滤细胞情况。🤒
library(scales)
plotColData(umi, x="sum", y="detected",
colour_by="discard", other_fields = "individual") +
facet_wrap(~individual) +
scale_x_continuous(labels = unit_format(unit = "k", scale = 1e-3))
我们再试一下将replicate
设为分面参数,看一下不同replicate
的过滤细胞情况。
plotColData(umi, x="sum", y="detected",
colour_by="discard", other_fields = "replicate") +
facet_wrap(~replicate) +
scale_x_continuous(labels = unit_format(unit = "k", scale = 1e-3))
这里可以看到不同批次间的细胞情况是有一定差异的。还是要做好质控啊!!!😉
最后祝大家早日不卷!~