🤩 scRNA-seq | 吐血整理的单细胞入门教程（质控与过滤）（七）

生信漫卷

发布于 2022-10-31 17:21:19

2.4K0

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))