🤩 scRNA-seq | 吐血整理的单细胞入门教程（质控与过滤）（七）
🤩 scRNA-seq | 吐血整理的单细胞入门教程（质控与过滤）（七）
1写在前面
当我们拿到表达矩阵后，需要对其进行质控（quality control, QC）去除质量较差的细胞，降低噪音，而后再进行数据分析。😘
2用到的包
rm(list = ls())
library(tidyverse)
library(scater)
library(SingleCellExperiment)
library(AnnotationDbi)
library(org.Hs.eg.db)
library(EnsDb.Hsapiens.v86)
3示例数据
这里我们用一下之前介绍的counts文件和annotation文件，然后通过SingleCellExperiment创建SingleCellExperiment格式的文件。
3.1 读入数据
counts <- read.delim("./molecules.txt", row.names = 1)
annotation <- read.delim("./annotation.txt", stringsAsFactors = T)
counts
annotation
3.2 创建SingleCellExperiment对象
# 注意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
4线粒体基因的处理
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)
5基础质控
这里我们用scater包的perCellQCMetrics和perFeatureQCMetrics函数，分别对细胞和特征（基因）进行分析。🤒
5.1 细胞信息
umi_cell <- perCellQCMetrics(umi,subsets=list(Mito=is_mito))
head(umi_cell)[1:4,1:4]
5.2 基因信息
umi_feature <- perFeatureQCMetrics(umi)
head(umi_feature)
5.3 整合信息
我们现在把上面的信息加到SingleCellExperiment文件中去，即umi。🥳
umi <- addPerCellQC(umi, subsets=list(Mito=is_mito))
umi <- addPerFeatureQC(umi)
6手动过滤
当然你也可以自己设定cutoff来过滤细胞或者基因。我们先来看一下分布情况。
6.1 counts分布
hist(
umi$total,
breaks = 100
)
abline(v = 25000, col = "red")
6.2 细胞分布
hist(
umi_cell$detected,
breaks = 100
)
abline(v = 7000, col = "red")
6.3 找不到具体的cutoff怎么办？
有时很难找到一个明确的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))
7scater包一键搞定过滤
7.1 过滤
使用scater包的quickPerCellQC函数进行过滤。
reasons <- quickPerCellQC(umi_cell,
sub.fields = c("subsets_Mito_percent", "altexps_ERCC_percent"))
colSums(as.matrix(reasons))
7.2 整合信息
我们把需要过滤的细胞信息加入到SingleCellExperiment文件中，即umi。👌
umi$discard <- reasons$discard
7.3 可视化1
这里可以看到，线粒体基因比例高的细胞，一般counts都比较低。👀
plotColData(umi, x="sum", y="subsets_Mito_percent", colour_by="discard")
7.4 可视化2
这里可以看到，线粒体基因比例高的细胞，一般检测到的基因都比较少。
plotColData(umi, x="sum", y="detected", colour_by="discard")
7.5 可视化3
这里可以看到，线粒体基因比例高的细胞，一般检测到的ERCC比例都比较高。😂
plotColData(umi, x="altexps_ERCC_percent", y="subsets_Mito_percent",colour_by="discard")
7.6 不同batch的区别
这里我们将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))
这里可以看到不同批次间的细胞情况是有一定差异的。还是要做好质控啊！！！😉
最后祝大家早日不卷!~