表达矩阵处理—表达质量的控制

书籍翻译

好的书籍是人类进步的阶梯，但有些人却找不到优秀的阶梯，为此我们开设了书籍翻译这个栏目，作为你学习之路的指路明灯；分享国内外优秀书籍，弘扬分享精神，做一个知识的传播者。

希望大家能有所收获！

书籍翻译历史目录

正

文

7. 清理表达矩阵

表达质量控制（UMI）

7.1.1

简介

一旦基因的表达被定量了，就将其概括为表达矩阵，其中每行对应于基因（或转录物），并且每列对应于单个细胞。通过检查该矩阵，去除在读取QC或mapping QC步骤中未检测到的劣质细胞。在此阶段未能移除低质量细胞可能会增加技术noise，这可能会模糊下游分析中感兴趣的生物信号。

由于目前没有用于执行scRNASeq的标准方法，因此本文将呈现的各种QC测量的预期值可以能在实验之间有显着变化。因此，为了执行QC，我们将寻找相对于数据集的其余部分异常的细胞，而不是与独立的质量标准进行比较。因此，在比较使用不同protocol收集的数据集之间的质量指标时应该小心。

7.1.2

Tung数据集

为了展示细胞QC过程，我们考虑一个数据集，这个数据集是从三个不同的个体产生的诱导多能干细胞（Tung等2017年）在Yoav Gilad在芝加哥大学的实验室中。该实验在Fluidigm  C1平台上进行，使用了方便独特的分子标识符（UMI的）和ERCC spike-ins。数据文件位于tung文件夹中。这些文件是15/03/16制作的原始文件的副本。我们将这些副本用于再现。

library(SingleCellExperiment)
library(scater)
options(stringsAsFactors = FALSE)

加载数据和注释：

molecules <- read.table("tung/molecules.txt", sep = "\t")
anno <- read.table("tung/annotation.txt", sep = "\t", header = TRUE)

检查表达矩阵的一小部分

head(molecules[ , 1:3])##                 NA19098.r1.A01 NA19098.r1.A02 NA19098.r1.A03
## ENSG00000237683             0             0             0
## ENSG00000187634             0             0             0
## ENSG00000188976             3             6             1
## ENSG00000187961             0             0             0
## ENSG00000187583             0             0             0
## ENSG00000187642             0             0             0head(anno)##   individual replicate well     batch     sample_id
## 1   NA19098       r1 A01 NA19098.r1 NA19098.r1.A01
## 2   NA19098       r1 A02 NA19098.r1 NA19098.r1.A02
## 3   NA19098       r1 A03 NA19098.r1 NA19098.r1.A03
## 4   NA19098       r1 A04 NA19098.r1 NA19098.r1.A04
## 5   NA19098       r1 A05 NA19098.r1 NA19098.r1.A05
## 6   NA19098       r1 A06 NA19098.r1 NA19098.r1.A06

数据由3个人和3次重复组成，因此总共有9个批次。

我们使用SingleCellExperiment（SCE）和scater包来进行标准化分析。首先，创建SCE对象：

umi <- SingleCellExperiment(
  assays = list(counts = as.matrix(molecules)),
  colData = anno
)

去除在任何细胞中都不表达的基因：

keep_feature <- rowSums(counts(umi) > 0) > 0
umi <- umi[keep_feature, ]

定义控制特征（基因） - ERCC spike-ins和线粒体基因（由作者提供）：

isSpike(umi, "ERCC") <- grepl("^ERCC-", rownames(umi))
isSpike(umi, "MT") <- rownames(umi) %in%
  c("ENSG00000198899", "ENSG00000198727", "ENSG00000198888",
  "ENSG00000198886", "ENSG00000212907", "ENSG00000198786",
  "ENSG00000198695", "ENSG00000198712", "ENSG00000198804",
  "ENSG00000198763", "ENSG00000228253", "ENSG00000198938",
  "ENSG00000198840")

计算质量指标：

umi <- calculateQCMetrics(
  umi,
  feature_controls = list(
      ERCC = isSpike(umi, "ERCC"),
      MT = isSpike(umi, "MT")
  )
)

7.1.3

细胞QC

7.1.3.1 库大小

接下来，我们考虑每个样本检测到的RNA分子总数（如果我们使用读取计数而不是UMI计数，则这将是reads的总数）。具有少量reads/分子的孔很可能已被破坏或未能捕获细胞，因此应被移除。

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

练习1

1. 我们的过滤器删除了多少个细胞？

2. 期望每个细胞的分子总数应遵循什么分布？

我们的答案

## filter_by_total_counts
## FALSE TRUE
##   46   818

7.1.3.2 检测到的基因

除了确保每个样品的足够测序深度外，我们还希望确保reads在转录组中的分布。因此，我们计算每个样品中检测到的独特基因的总数。

hist(
  umi$total_features,
  breaks = 100
)
abline(v = 7000, col = "red")

从图中我们得出结论，大多数细胞具有7,000-10,000个检测到的基因，这对于高深度scRNA-seq是正常的。然而，这取决于实验方案和测序深度。例如，基于液滴的方法或具有较低测序深度的样品通常每个细胞检测较少的基因。上图中最显着的特征是分布左侧的“重尾”。如果检测率在细胞中相等，则分布应近似于高斯分布。因此，我们去除了分布尾部的那些细胞（少于7,000个检测到的基因）。

练习2

我们的过滤器删除了多少个细胞？

我们的答案

## filter_by_expr_features
## FALSE TRUE
##   116   748

7.1.3.3 ERCCs和MTs

细胞质量的另一个衡量标准是ERCC spike-in RNA和内源RNA 之间的比例。该比率可用于估计捕获细胞中RNA的总量。具有高水平的spike-in RNA的细胞具有低起始量的RNA，可能是由于细胞死亡或受到应激，这可能导致RNA降解。

plotPhenoData(
  umi,
  aes_string(
      x = "total_features",
      y = "pct_counts_MT",
      colour = "batch"
  )
)lotPhenoData(
  umi,
  aes_string(
      x = "total_features",
      y = "pct_counts_ERCC",
      colour = "batch"
  )
)

上述分析表明，来自NA19098.r2批次的大多数细胞具有非常高的ERCC / Endo比率。实际上，作者已经证明该批次包含较小尺寸的细胞。

练习3

创建用于去除批次NA19098.r2的过滤器和具有高线粒体基因表达的细胞（> 10％的总数在一个细胞内）。

我们的答案

## filter_by_ERCC
## FALSE TRUE
## 96   768## filter_by_MT
## FALSE TRUE
## 31   833

练习4

如果您正在检查包含不同大小细胞的数据集（例如，正常细胞和衰老细胞），您期望在ERCC与计数图中看到什么？

回答

您可能会看到对应于较小细胞（正常）的组，ERCC读数的分数高于对应于较大细胞（衰老）的一个单独的组。

7.1.4

细胞过滤

7.1.4.1 手册

现在我们可以根据之前的分析定义一个细胞过滤器:

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

7.1.4.2 自动的

可用的另一种选择scater是在一组QC指标上进行PCA，然后使用自动异常检测来识别可能存在问题的细胞。

默认情况下，以下度量标准用于基于PCA的离群值检测：

• pct_counts_top_100_features

• total_features

• pct_counts_feature_controls

• n_detected_feature_controls

• log10_counts_endogenous_features

• log10_counts_feature_controlsscater首先创建一个矩阵，其中行代表细胞，列代表不同的QC指标。这里，PCA图提供了按质量度量排序的单元格的2D表示。然后使用来自mvoutlier包的方法检测异常值。

umi <- plotPCA(
  umi,
  size_by = "total_features",
  shape_by = "use",
  pca_data_input = "pdata",
  detect_outliers = TRUE,
  return_SCE = TRUE
)table(umi$outlier)##
## FALSE TRUE
##   819   45

7.1.5

比较过滤器

练习5

比较默认，自动和手动细胞过滤器。从这些过滤器绘制异常细胞的维恩图。

提示：使用limma包中的函数vennCounts和vennDiagram函数来制作维恩图。

回答

##
## Attaching package: 'limma'## The following object is masked from 'package:scater':
##
##     plotMDS## The following object is masked from 'package:BiocGenerics':
##
##     plotMA

7.1.6

基因分析

7.1.6.1 基因表达

除了去除质量差的细胞外，通常最好排除我们怀疑技术假象可能会导致结果偏差的基因。此外，检查基因表达谱可以提供关于如何改进实验程序的见解。

考虑前50个表达基因消耗的reads通常是有益的。

plotQC(umi, type = "highest-expression")

7.1.6.2基因过滤

通常移除表达水平被认为是“不可检测的”的基因是一个好主意。如果至少两个细胞含有超过1个来自该基因的转录物，我们将基因定义为可检测的。如果我们考虑read计数而不是UMI计数，则合理的阈值是要求至少两个细胞中至少五条reads。但是，在这两种情况下，阈值都很大程度上取决于测序深度。重要的是要记住，必须在细胞过滤后过滤基因，因为某些基因可能只能在质量较差的细胞中检测到（注意colData(umi)$use过滤器应用于umi数据集）。

filter_genes <- apply(
  counts(umi[ , colData(umi)$use]),
  1,
  function(x) length(x[x > 1]) >= 2
)
rowData(umi)$use <- filter_genestable(filter_genes)## filter_genes
## FALSE TRUE
## 4660 14066

根据细胞类型，方案和测序深度，其他阈值可能也是适当的。

7.17

保存数据

QCed数据集的维度（不要忘记我们上面定义的基因过滤器）：

dim(umi[rowData(umi)$use, colData(umi)$use])
## [1] 14066 657

让我们创建一个带有log转换计数的附加插槽（我们将在下一章中使用它）并从reducedDim插槽中删除已保存的PCA结果：

assay(umi, "logcounts_raw") <- log2(counts(umi) + 1)
reducedDim(umi) <- NULL

保存数据：

saveRDS(umi, file = "tung/umi.rds")

7.1.8

大练习

使用相同Blischak数据的读取计数执行完全相同的QC分析。使用tung/reads.txt文件加载reads。完成后，请将您的结果与我们的结果进行比较（下一章）。

7.19

sessionInfo()

## R version 3.4.3 (2017-11-30)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Debian GNU/Linux 9 (stretch)
##
## Matrix products: default
## BLAS: /usr/lib/openblas-base/libblas.so.3
## LAPACK: /usr/lib/libopenblasp-r0.2.19.so
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=C
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] parallel stats4 methods stats graphics grDevices utils
## [8] datasets base
##
## other attached packages:
## [1] limma_3.34.9 scater_1.6.3
## [3] ggplot2_2.2.1 SingleCellExperiment_1.0.0
## [5] SummarizedExperiment_1.8.1 DelayedArray_0.4.1
## [7] matrixStats_0.53.1 Biobase_2.38.0
## [9] GenomicRanges_1.30.3 GenomeInfoDb_1.14.0
## [11] IRanges_2.12.0 S4Vectors_0.16.0
## [13] BiocGenerics_0.24.0 knitr_1.20
##
## loaded via a namespace (and not attached):
## [1] backports_1.1.2 plyr_1.8.4 lazyeval_0.2.1
## [4] sp_1.2-7 shinydashboard_0.6.1 splines_3.4.3
## [7] digest_0.6.15 htmltools_0.3.6 viridis_0.5.0
## [10] magrittr_1.5 memoise_1.1.0 cluster_2.0.6
## [13] prettyunits_1.0.2 colorspace_1.3-2 blob_1.1.0
## [16] rrcov_1.4-3 xfun_0.1 dplyr_0.7.4
## [19] RCurl_1.95-4.10 tximport_1.6.0 lme4_1.1-15
## [22] bindr_0.1 zoo_1.8-1 glue_1.2.0
## [25] gtable_0.2.0 zlibbioc_1.24.0 XVector_0.18.0
## [28] MatrixModels_0.4-1 car_2.1-6 kernlab_0.9-25
## [31] prabclus_2.2-6 DEoptimR_1.0-8 SparseM_1.77
## [34] VIM_4.7.0 scales_0.5.0 sgeostat_1.0-27
## [37] mvtnorm_1.0-7 DBI_0.7 GGally_1.3.2
## [40] edgeR_3.20.9 Rcpp_0.12.15 sROC_0.1-2
## [43] viridisLite_0.3.0 xtable_1.8-2 progress_1.1.2
## [46] laeken_0.4.6 bit_1.1-12 mclust_5.4
## [49] vcd_1.4-4 httr_1.3.1 RColorBrewer_1.1-2
## [52] fpc_2.1-11 modeltools_0.2-21 pkgconfig_2.0.1
## [55] reshape_0.8.7 XML_3.98-1.10 flexmix_2.3-14
## [58] nnet_7.3-12 locfit_1.5-9.1 labeling_0.3
## [61] rlang_0.2.0 reshape2_1.4.3 AnnotationDbi_1.40.0
## [64] munsell_0.4.3 tools_3.4.3 RSQLite_2.0
## [67] pls_2.6-0 evaluate_0.10.1 stringr_1.3.0
## [70] cvTools_0.3.2 yaml_2.1.17 bit64_0.9-7
## [73] robustbase_0.92-8 bindrcpp_0.2 nlme_3.1-129
## [76] mime_0.5 quantreg_5.35 biomaRt_2.34.2
## [79] compiler_3.4.3 pbkrtest_0.4-7 beeswarm_0.2.3
## [82] e1071_1.6-8 tibble_1.4.2 robCompositions_2.0.6
## [85] pcaPP_1.9-73 stringi_1.1.6 highr_0.6
## [88] lattice_0.20-34 trimcluster_0.1-2 Matrix_1.2-7.1
## [91] nloptr_1.0.4 pillar_1.2.1 lmtest_0.9-35
## [94] data.table_1.10.4-3 cowplot_0.9.2 bitops_1.0-6
## [97] httpuv_1.3.6.1 R6_2.2.2 bookdown_0.7
## [100] gridExtra_2.3 vipor_0.4.5 boot_1.3-18
## [103] MASS_7.3-45 assertthat_0.2.0 rhdf5_2.22.0
## [106] rprojroot_1.3-2 rjson_0.2.15 GenomeInfoDbData_1.0.0
## [109] diptest_0.75-7 mgcv_1.8-23 grid_3.4.3
## [112] class_7.3-14 minqa_1.2.4 rmarkdown_1.8
## [115] mvoutlier_2.0.9 shiny_1.0.5 ggbeeswarm_0.6.0

 相关推荐

关于表达矩阵更多可视化内容，请参考群主：表达矩阵可视化大全

如果你对单细胞转录组研究感兴趣，但又不知道如何入门，也许你可以关注一下下面的课程

• 全国巡讲第9、10站-武汉和成都（生信技能树爆款入门课）（火热报名中）

单细胞天地欢迎你

单细胞天地

生信技能树