Hemberg-lab单细胞转录组数据分析（九）- Scater包单细胞过滤

往期系列

Hemberg-lab单细胞转录组数据分析（一）

Hemberg-lab单细胞转录组数据分析（二）

Hemberg-lab单细胞转录组数据分析（三）

Hemberg-lab单细胞转录组数据分析（四）

Hemberg-lab单细胞转录组数据分析（五）

Hemberg-lab单细胞转录组数据分析（六）

Hemberg-lab单细胞转录组数据分析（七）-导入10X和SmartSeq2数据Tabula Muris

Hemberg-lab单细胞转录组数据分析（八）- Scater包输入导入和存储

收藏|北大生信平台"单细胞分析、染色质分析"视频和PPT分享

该如何自学入门生物信息学

生物信息之程序学习

收藏|你想要的生信学习系列教程-宝典在手，生信无忧

细胞质控

文库大小

查看每个样品(细胞）检测到的总分子数 (UMI count)或总reads数 (reads count)，拥有很少的reads或分子数的样品可能是细胞破损或捕获失败，应该移除。

hist(
    umi$total_counts,
    breaks = 100
)
abline(v = 25000, col = "red")

练习：

1. 我们的过滤移除了多少细胞？

2. 每个细胞中检测到的分子数的分布预期是怎样的？

答案

filter_by_total_counts <- (umi$total_counts > 25000)
table(filter_by_total_counts)
## filter_by_total_counts
## FALSE  TRUE 
##    46   818

检测到的基因数

除了确保每个样品的测序深度，也需要保证测序reads在转录本中分布均衡，而不是集中在少数高表达的基因上。每个样品检测到的基因数也是衡量样品质量好坏的一个标准。

# 原文这个地方有误，可能是版本问题
hist(
    umi$total_features_by_counts,
    breaks = 100
)
abline(v = 7000, col = "red")

从图中可以看出，大部分细胞能检测到7,000-10,000基因，这对高深度scRNA-seq是正常的。当然这个受测序深度和实验方法的影响。比如居于droplet的方法或样品测序深度低时每个细胞检测到的基因数要少一些，表现在图上是，左侧拖尾严重。如果细胞之间的基因检出率相当，应该符合正态分布。因此选择移除分布尾部的细胞 (本例中是检测出的基因数少于7000的细胞)。

练习2: 移除了多少细胞？

答案

## filter_by_expr_features
## FALSE  TRUE 
##   116   748

ERCCs和MTs

另外一个测量细胞质量的方式是比较ERCC spike-in测到的reads数与内源转录本测到的reads数的比例，可以衡量捕获到的内源性RNA的总量。如果spike in的reads数很高，则表示起始内源性RNA总量低，可能是由于细胞死亡或胁迫诱导的RNA降解导致的，也有可能是细胞体积小。

plotColData(    umi,    x = "total_features_by_counts",    y = "pct_counts_MT",    colour = "batch")

plotColData(
    umi,
    x = "total_features_by_counts",
    y = "pct_counts_MT",
    colour = "batch"
)

上图显示来源于NA19098.r2批次的细胞有较高的ERCC/内源RNA比例。作者在文章中证实这一点，说这个批次的细胞体积小。

练习 3：移除NA19098.r2批次的细胞和线粒体基因表达量超过10%的细胞。

答案

filter_by_ERCC <- umi$batch != "NA19098.r2"
table(filter_by_ERCC)
## filter_by_ERCC
## FALSE  TRUE 
##    96   768
filter_by_MT <- umi$pct_counts_MT < 10
table(filter_by_MT)
## filter_by_MT
## FALSE  TRUE 
##    31   833

练习 4: 如果研究的数据集细胞大小不同（正常细胞、衰老细胞），那么ERCC与内源基因被测到的比例会是怎么的分布？

答案：小的细胞 (normal)比大的细胞(senescent,衰老)有更高比例的ERCC reads。

细胞过滤

手动过滤

基于前面的分析定义一个过滤器，不满足任何一个条件的细胞都过滤掉:

umi$use <- (
    # sufficient features (genes)
    filter_by_expr_features &
    # sufficient molecules counted
    filter_by_total_counts &
    # sufficient endogenous RNA
    filter_by_ERCC &
    # remove cells with unusual number of reads in MT genes
    filter_by_MT
)
table(umi$use)
## 
## FALSE  TRUE 
##   207   657

自动过滤

scater提供了一个根据质控数据进行PCA分析进而自动挑出异常细胞的方法。默认，下面这些统计量将用于PCA异常细胞检测的分析：

• pct_counts_top_100_features

• total_features_by_counts

• pct_counts_feature_controls

• n_detected_feature_controls

• log10_counts_endogenous_features

• log10_counts_feature_controls

scater首先生成一个行是细胞，列是细胞中对应的上述质控数据的值，然后使用mvoutlier包筛选质控数据与大部分细胞不同的样品定义为低质量细胞。 package on the QC metrics for all cells. This will identify cells that have substantially different QC metrics from the others, possibly corresponding to low-quality cells. We can visualize any outliers using a principal components plot as shown below:

umi <- runPCA(umi, use_coldata = TRUE, 
              detect_outliers = TRUE)
reducedDimNames(umi)
## [1] "PCA_coldata"

鉴定结果存储于umi变量的$outlier部分，指示细胞是否被判断未异常细胞。自动异常细胞检测是很有意义的，可以作为工厂化大批量模式使用，但特异性的手动检测数据集和根据结果、实验调整过滤是推荐的方式。

table(umi$outlier)

## ## FALSE  TRUE 
##   791    73

绘制PCA结果展示异常细胞分布：

plotReducedDim(umi, use_dimred = "PCA_coldata",
               size_by = "total_features_by_counts", shape_by = "use", 
               colour_by="outlier")

手动过滤和自动过滤比较

练习 5: 绘制Venn图比较自动和手动两个方式检测出的异常细胞

提示: 使用limma包里的vennCounts和vennDiagram函数绘制。生信宝典说，使用高颜值在线绘图工具http://www.ehbio.com/ImageGP更方便。

答案

还有一种方式是使用中位数绝对偏差作为判断样品异常的标准。以测序文库大小为例，假设样品中的Total read count是，所有样品中Total read count的中位数是，那么样品 Total read count的绝对偏差就是。 的样品会被移除 （移除测序深度低的样品）。为了增强过滤的鲁棒性，依据样品测序的文库大小和检测到的基因数目过滤时会先对相应对数值进行对数转换。依据ERCC spike-in基因的比例和线粒体基因的比例过滤时，的样品会被移除 （移除检测到的内源基因少的样品）。