10X单细胞（10X空间转录组）数据降噪分析（Imputation）之SAVER-X

hello，大家好，今天给大家分享一个用于10X单细胞数据降噪的分析软件，SAVER-X，其实这个软件的前身是SAVER，关于SAVER的文章在SAVER: gene expression recovery for single-cell RNA sequencing | Nature Methods，2018年6月发表于nature methods，影响因子30分，我们简单回顾一下SVAER，不展开，感兴趣大家可以看一下文章。

scRNA-seq数据的imputation(Imputation of ScRNA-Seq Data)

单细胞RNA-seq数据通常包含许多由于原始RNA扩增失败而导致的缺失值或dropouts。（注解： dropouts是指空beats）。dropout 的比例是跟 protocol 相关，而且与每个细胞测序的reads数密切相关。dropout event 增加了细胞间的变异。导致信号对每个基因的影响，模糊了基因-基因关系的检测。由于在scRNA-seq中可能无法检测到大量真正表达的转录本，因此dropout 可能会很大程度上影响下游的分析。Imputation是一种用替代值替换缺失数据(dropouts)的有效方法。虽然已经提出了一些方法来估算bulk RNA-seq数据，但它们并不直接适用于scRNA-seq数据。最近为scRNA-seq开发了几种imputation方法。包括SAVER，MAGIC，ScImpute，DrImpute，AutoImpute。SAVER是一个基于贝叶斯的模型，用于基于umi的scRNA-seq数据来恢复所有基因的真实表达水平。MAGIC通过构建基于马尔可夫亲和度的基因表达图来进行基因表达的计算。ScImpute的开发者认为，SAVER和MAGIC可能会导致未受dropout影响的基因的表达变化，而ScImpute可以利用其他类似细胞中不太可能受dropout影响的相同基因的信息，在不引入新的偏差的情况下计算dropout值。DrImpute是一种基于集群的方法，可以有效地将dropout中的零从真正的零中分离出来。AutoImpute是一种基于自编码的方法，它通过学习scRNA-seq数据的固有分布来寻找（impute）缺失的值。最近，Zhang等人评估了不同的imputation方法，发现这些方法的性能与其模型假设和可扩展性相关。

简单回顾一下代码和原理

主函数saver

首先我们先读入文件，数据下载地址为：

data.path <- "expression_mRNA_17-Aug-2014.txt"

# Need to remove first 10 rows of metadata
raw.data <- read.table(data.path, header = FALSE, skip = 11, row.names = 1, 
                       check.names = FALSE)
cortex <- as.matrix(raw.data[, -1])

cellnames <- read.table(data.path, skip = 7, nrows = 1, row.names = 1, 
                        stringsAsFactors = FALSE)
colnames(cortex) <- cellnames[-1]

那么我们的数据长这样：

行名代表不同的细胞，列名代表不同的基因

然后参观下主函数saver

  x = cortex
  do.fast = TRUE
  ncores = 1
  size.factor = NULL
  npred = NULL
  pred.cells = NULL
  pred.genes = NULL
  pred.genes.only = FALSE
  null.model = FALSE
  mu = NULL
  estimates.only = FALSE
  
  x <- clean.data(x)
  np <- dim(x)
  ngenes <- as.integer(np[1])
  ncells <- as.integer(np[2])

  message(ngenes, " genes, ", ncells, " cells")

  # assign size factor
  sf.out <- calc.size.factor(x, size.factor, ncells)
  sf <- sf.out[[1]]
  scale.sf <- sf.out[[2]]

  gene.names <- rownames(x)
  cell.names <- colnames(x)

  # Set up parallel backend
  cl.create <- FALSE
  ncores <- getDoParWorkers()
  
  # if prior means are provided

  # assign pred.cells and pred.genes
  pred.cells <- get.pred.cells(pred.cells, ncells)
  pred.genes <- get.pred.genes(x, pred.genes, npred, ngenes)
  npred <- length(pred.genes)


  good.genes <- which(Matrix::rowMeans(x) >= 0.1)
#筛选好的单细胞表达矩阵定义为x.est
  x.est <- t(as.matrix(log(sweep(x[good.genes, ] + 1, 2, sf, "/"))))
  
  out <- saver.fit(x, x.est, do.fast, ncores, sf, scale.sf, pred.genes,
                     pred.cells, null.model, ngenes = nrow(x),
                     ncells = ncol(x), gene.names, cell.names,
                     estimates.only)
  }


  class(out) <- "saver"
  out

我们观察得知，最终的结果需要经过saver.fit函数的计算，我们通过观察saver.fit，发现它的结果需要经过calc.estimate函数的计算，而calc.estimate函数里面的核心函数又是expr.predict
所以当我们目光集中在expr.predict时：

  x = x.est
  y = y[1, ]  ##y[1, ]相当于单细胞表达矩阵中的一个基因在各个细胞中的表达
  lambda.max = NULL
  lambda.min = NULL

  if (is.null(lambda.max)) {
    cv <- tryCatch(
# 拟合泊松分布，并获得最佳参数 λ
      suppressWarnings(glmnet::cv.glmnet(x[pred.cells, ], t(y[pred.cells]),
                                         family="poisson", dfmax = 300,
                                         nfolds = 5)),
      error = function(cond) {
        return(NA)
      }
    )

# newx表示抽样的背景集
    mu <- c(predict(cv, newx = x, s = "lambda.min",
                    type="response"))
    lambda.max <- cv$lambda[1]
    lambda.min <- cv$lambda.min
    min.ind <- which(cv$lambda == cv$lambda.min)
    sd.cv <- (cv$cvm[1] - cv$cvm[min.ind]) / cv$cvsd[min.ind]
  } 

那么我们所输入的数据是啥呢？
首先是x.est，这个其实可以理解为经过筛选以后好的基因的单细胞表达矩阵：

那么y[1,]（其实只是取了其中一行做例子，这里是按行遍历）其实就是其中的一个基因在各个细胞中的表达

这里的cv.glmnet表示的是经过交叉验证来拟合泊松分布，换言之就是每个细胞的基因表达量的分布是满足于泊松分布的（推荐看单细胞组学的单细胞数据特征），其中x[pred.cells, ]，相当于背景集的总体分布；而 t(y[pred.cells])相当于其中的一次抽样（pred.cells相当于选中的细胞），相当于其中的一个样本（矩估计的思想）
第二步就是确定最佳λ的值（λ是泊松分布的参数），然后用该泊松分布模型来预测该基因在各细胞中的表达情况（相当于在泊松分布里面抽样，newx表示抽样的背景集，并在该分布中抽一次样，作为该基因矫正后的表达量），作为该基因矫正后的表达量

那么我们说 y[1,]相等于其中的一个基因，那么这一个基因联合背景集x.est拟合一个泊松分布，再从里面抽样作为抽样预测的表达值mu
遍历所有基因（比方说由 n 个基因），那么就会有 n 个泊松分布，这 n 个泊松分布分别对应于每一个基因，然后从每一个基因各自的分布里面抽样，作为每一个基因抽样预测的表达值mu

如上图所示，每一个基因对应一个泊松分布，并在这个分布中重新抽样，作为每个基因抽样预测的表达值mu

当然，作者对矫正后的基因在各细胞中的表达做了如下定义：

定义式

其中式子中：
λ 表示第c个细胞，第g个基因的矫正后的表达量
Ygc表示抽样完的第c个细胞，第g个基因表达量，即为mu

mu表示抽样预测后的基因表达量，为一列向量值：

而y[1,]表示抽样前的基因表达量

这么矫正以后感觉似乎表达量提高了一些，0值少了不少

之后SAVER又更具拟合的先验分布做出计算（以下代码即为计算定义式的代码）：

##prior.beta，prior.alpha是通过极大似然估计求解的
# mu 表示抽样预测后的基因表达量，为一列向量值
post.alpha <- prior.alpha + mu
# sf 为标准化因子
post.beta <- prior.beta + sf
lambda.hat <- post.alpha/post.beta

那么最后每一个函数都返回它们的lambda.ha，即为最后的结果（lambda.ha是某一个基因在各个细胞中的表达量，那么遍历所有基因就可以得到最终整个矫正后的表达矩阵）,这个大家尝试一下就可以了。

来到我们今天的重点SVAER-X,大家看名字就知道，这个软件是SAVER的改进版,文章在Data denoising with transfer learning in single-cell transcriptomics,2019年8月发表于nature methods，影响因子30分,我们简单分享一下文献，最后看一下代码。

摘要

单细胞转录组测序（single-cell RNA-seq, scRNA-seq）数据有高噪音和稀疏的特点。原文作者展示了跨数据集的迁移学习可显著提高数据的质量。通过将深度自动编码器与贝叶斯模型相结合，原文开发的SAVER-X软件可从不同实验室、不同条件和不同物种的数据中提取可迁移的基因关系，以对新的目标数据集进行降噪。

前言

在scRNA-seq的研究中，技术噪音让细胞状态之间的精确识别模糊化，而且低表达基因不能被准确量化。当测序深度较低或感兴趣的细胞类型较少时，现有的降噪方法往往表现不佳；而且这些方法也忽略了公共的数据集，这些公共数据集可能包含有助于降噪的相关信息。现在逐渐建立了小鼠的细胞图谱，很快就会有人体每个器官的详细图谱。可公开利用的scRNA-seq数据集包含与新产生数据相关的细胞类型和基因特征信息。然而，目前还不清楚如何跨平台、跨物种、跨组织地借用信息。此外，这种迁移学习方法不能引入bias或迫使新数据失去其独特的特征。

结果

原文作者开发了一种新的降噪方SAVER-X，该方法通过表达复原利用外部数据进行单细胞分析；它将贝叶斯分层模型与一个可预处理的深度自动编码器相结合。尽管神经网络已经成为了其他单细胞方法的基础，但现有的工具只对周围的数据起作用。此外，本文中广泛的基准测试和数据强调，除了SAVER-X软件的前身SAVER之外，大多数方法对真实基因表达都产生有偏估计，并引入了特殊的基因相关性。SAVER-X是建立在SAVER软件的核心模型基础上，结合自动编码器后端和两阶段训练机制来利用公共数据资源。

SAVER-X将Y的变异分解为:(i)一个可预测的结构化成分(λ)，解释基因间的共享变异；(ii)不可预测的细胞水平波动，这些波动在具有基因特异性分散的基因之间是独立的；和(iii)技术噪音。SAVER-X估计未观察到的真实基因表达，即λ加上独立的细胞水平波动。λ是通过自动编码器(图1b)获得的，该编码器的权重首先在来自相同组织或相似类型组织的细胞上进行预处理，从公共存储库中提取(“预处理数据”；图1a)。权重随后被更新以适应目标数据。这种两阶段训练方案允许适应性保留可迁移特征。人类和小鼠共享许多核心的细胞类型和必要的生物学通路。

图1 | SAVER-X软件的迁移学习框架概要。（图片来源：Wang et al. Nature Methods，2019）a为自动编码器预处理步骤。b为SAVER-X的工作流程。对于具有UMI计数矩阵的目标数据，SAVER-X在没有选定预处理模型的情况下用自动编码器训练目标数据(项目A)，然后使用交叉验证过滤不可预测的基因(项目B)，并使用经验贝叶斯收缩估计最终降噪值(项目C)。

为了可以跨物种学习，SAVER-X中的自动编码器包含了人和小鼠之间的共同网络。此外，SAVER-X采用基于交叉验证的基因过滤和贝叶斯收缩来保留目标数据集特有的表达模式(图1b)。交叉验证识别出与自动编码器不匹配的基因，自动编码器的预测被它们的目标数据均值所取代。贝叶斯收缩方法通过计算加权预测值的平均值和观察到的数值来得到最终的降噪值。

原文作者首先探索了SAVER-X在免疫系统不同细胞类型上进行迁移学习的好处和局限性。在没有流动分类的情况下，scRNA-seq数据中的组织浸润免疫细胞很少，因此，如果没有外部数据的帮助，降噪很有挑战性。软件作者检测了SAVER-X是否对来自人类细胞图谱(HCA)项目(来自脐带血和骨髓的500，000个免疫细胞)和10X基因组学(200，000个外周血单核细胞)的数据进行了预处理，从而有意义地提高了来自健康和疾病组织的免疫细胞的数据质量。与此同时，软件作者把SAVER-X和现有的降噪方法进行了对比，降噪方法是针对来自九种非重叠免疫细胞类型的一组纯化细胞进行的。

对T细胞亚型的可靠鉴定，对组织免疫环境的定性至关重要，然而，在原始的单细胞转录组测序数据中，T细胞亚型常常被混在一起(图2a)。原文作者通过对每种细胞类型随机选择100个细胞创建了一个测试数据集，发现SAVER-X不仅显著提高了肝癌T细胞亚型之间的分离，而且与其他降噪方法相比，还提高了细胞类型识别的准确性(通过矫正后的兰德指数(ARI)来测量)。包含相似细胞类型的数据集显示了迁移方法的有效性，通过降低目标数据中细胞数量或测序深度，原文作者也观察到了更加明显的效果。在极端情况下，即使是覆盖范围仅为60个唯一分子标识符(UMIs)的细胞(通常被丢弃在当前管道中)，也可以通过迁移学习来拯救，以揭示有用的信息。

图2 |人类免疫细胞的SAVER-X降噪。（ 图片来源：Wang et al. Nature Methods，2019）

为了理解迁移学习的局限性，原文作者评估了降噪精度，与预处理数据集和目标数据集之间细胞类型的相似性关系。迁移学习是否能有效地降噪预处理数据中不存在的细胞状态呢？考虑上面分析的纯化T细胞。当对所有的T细胞亚型进行SAVER-X预处理时，标记基因的聚类和表达量化得到了显著改善(图2b)。然而，即使细胞类型(CD4+调节性T (Treg)细胞)完全在预处理中缺失，SAVER-X也改善了这种“新”细胞类型的鉴定和标记基因定量。而且，为了确定预处理数据中细胞类型的富集是否提高了目标数据中的降噪精度，原文作者对富集了Treg细胞的T细胞进行了SAVER-X预处理，发现会产生任何明显的差异。因此，SAVER-X不需要预处理数据和目标数据之间的细胞类型完美匹配，重要的是，可以提高预处理数据中不存在的新细胞类型的量化。

作为辅助手段，原文作者还使用CITE-seq技术检测关键免疫标记的基因蛋白相关性。在CITE-seq中发现蛋白质丰度和它们的同源基因RNA表达之间的相关性非常低。而且还发现，对于PBMC和脐带血单个核细胞(CBMC)的CITE-seq数据，SAVER-X的降噪表达估计(在HCA和PBMC 10X基因组学上进行了预处理)与它们的蛋白质产物有明显更高的相关性。与其他方法相比，当目标数据集包含100和1000个细胞时，SAVER-X都提高了所有标记的相关性(图2c)。然而，对于具有8，000个细胞的较大数据集，预处理并没有显示出明显的效果。

接下来，原文作者探讨了SAVER-X是否能有效地从健康的免疫细胞中去学习降噪。与非预处理模型相比，对来自健康组织的免疫细胞(HCA和PBMC 10X基因组学)的预处理，使能够更好地表征多受试者中肿瘤浸润免疫细胞的类型(图2d)。同时，迁移学习后，肿瘤相关免疫细胞亚群仍可识别。特别地，SAVER-X保留了升高的免疫球蛋白，这是该疾病特异细胞状态下的标志(图2e)。该细胞群在健康组织中不存在，通过LYZ等标记物进一步验证了该免疫类型的存在。配对肿瘤和正常组织之间的细胞亚型特异基因表达差异，也在所有具有配对组织的受试者中，和对于免疫监测至关重要的两种细胞类型中得以保留(图2f)。

最后，原文作者使用了来自小鼠和人类发育腹侧中脑细胞scRNA-seq数据，进行跨物种迁移学习。对人类数据集中10%的读段（reads）进行降采样，将其降低到每个细胞有中值为452 UMIs，并利用原始数据作为衡量降噪精度的参考。将人类细胞随机分成两组，对一组的reads进行下采样（down-sampled），并保留另一组进行预处理。与未预处理的相比，在匹配的小鼠脑细胞上的SAVER-X预处理，导致人类脑细胞分类的有了显著地改善(图3a)。与单独对人类细胞进行预处理相比，对人和小鼠细胞预处理进一步提高了降噪精度。此外，对腹侧中脑以外区域的细胞进行SAVER-X预处理是有益的，而且对三个人类非UMI数据集以及小鼠细胞进行预处理也是有益的。这些实验证明了跨物种迁移学习在一般和实际环境中的优点，尽管在这些环境中，解剖区域和实验方案可能在预处理数据和目标数据之间有所不同。

然后，原文作者仔细检查了一个对小鼠数据进行预处理的模型，判断是否会对人类特异表达基因的估计产生偏差。他们计算了人类和小鼠之间，细胞类型特异平均表达的对数倍数变化，并鉴定了两个物种之间四种细胞类型的差异表达基因。用处理小鼠细胞的SAVER-X对下采样的人类数据进行降噪，结果保留了对数的倍数变化(图3b)。相反，仅仅依靠自动编码器，而没有交叉验证和收缩，则减少了一些基因的倍数变化。与其他方法不同，SAVER-X还保留了细胞间基因的差异(图3c)。

图3 |小鼠和人类脑细胞之间的迁移学习。（ 图片来源：Wang et al. Nature Methods，2019）

综上所述，结果表明SAVER-X的框架可以利用现有数据，来提高新scRNA-seq数据集的质量。SAVER-X的核心是在一系列研究设计中训练一个深层神经网络，并将该模型应用于新数据，以加强共享的生物模式。迁移学习将scRNA-seq数据分析的方法，从特定研究的统计建模过程转变为跨研究数据集成和信息共享的自动化过程。

我们来看一下代码

Basic Usage

Our input can be either a data file or a dense/sparse matrix. The data file can be '.txt', '.csv' or '.rds' file. The '.rds' file can store either a matrix or a sparse matrix of class 'dgCMatrix'. As a toy example, you may download one of our demo datasets on the web server, shekhar_downsampled.csv, the down-sampled mouse retina data from here. SAVER-X will generate a new folder based on the current time to store intemediate files and final results and return the file name storing the denoised matrix.

SAVER-X without pretraining

library(SAVERX)
file <- saverx("./testdata/shekhar_downsampled.csv")
denoised.data <- readRDS(file)

SAVER-X with a pretrained model

For the demo dataset, we have a pre-trained model for the mouse retina, please download the file, mouse_Retina.hdf5, and you may save it in './mouse_retina.hdf5'

library(SAVERX)
file <- saverx("./testdata/shekhar_downsampled.csv", data.species = "Mouse", 
use.pretrain = T, pretrained.weights.file = "./mouse_retina.hdf5", model.species = "Mouse")
denoised.data <- readRDS(file)

For both cases, object 'denoised.data' is your final denoised matrix. When dealing with large datasets, you can set 'is.large.data = T' to reduce RAM. Set 'clearup.python.session = F' if you plan to run SAVER-X more than once in your R session.

生活很好，有你更好~~~