SAVER 简介

这是一篇讲述如果单细胞测序由于深度不够导致各个基因表达量偏低怎么办的故事：《SAVER: Gene expression recovery for single-cell RNA sequencing》
也就是通过SAVER使表达量升高一些，使基因表达量为 0 的情况减少

主函数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

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

定义式

其中式子中：

1. λ 表示第c个细胞，第g个基因的矫正后的表达量
2. 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是某一个基因在各个细胞中的表达量，那么遍历所有基因就可以得到最终整个矫正后的表达矩阵）