前言

我们在使用DESeq2做下游分析的时候，往往会利用到DESeq2的标准化函数：rlog和vst函数进行标准化，而这两种标准化方式是有一定区别的：

rlog要求

rlog更适用于样本量小于30的数据，相应的vst负责的是样本量大于30的数据

负二项分布模型

对于计数型响应变量，类别（因子类）决策变量，往往利用Poisson回归来拟合模型，但是Poisson回归对离散值的要求很大，因此大多少情况下采用的是负二项回归来拟合模型，而经过负二项分布拟合的data的系数通常采用的是wald test来检验其显著性

那么判断组别间的差异，可以在两个组之间建立线性模型（负二项回归模型）来判断：

以下3个回归系数均是与control组进行比较，那么系数为正数，则代表相对于control组是上调的；为负值，则代表相对于control组是下调的

rlog核心函数

我们先看看整体的rlog函数布局：

  rld <- rlogData(object, intercept, betaPriorVar)
  se <- SummarizedExperiment(
           assays = rld,
           colData = colData(object),
           rowRanges = rowRanges(object),
           metadata = metadata(object))
  dt <- DESeqTransform(se)
  attr(dt,"betaPriorVar") <- attr(rld, "betaPriorVar")
  if (!is.null(attr(rld,"intercept"))) {
    mcols(dt)$rlogIntercept <- attr(rld,"intercept")
  }
  dt

通过源代码我们知道，我们要获得Normalized的data，其核心代码为：rld <- rlogData(object, intercept, betaPriorVar)，而rld返回的即为我们的Normalized table

因此 我们有必要查看一下rlogData的代码：

  lambda <- 1/rep(betaPriorVar,ncol(modelMatrix))
  # except for intercept which we set to wide prior
  if ("Intercept" %in% modelMatrixNames) {
    lambda[which(modelMatrixNames == "Intercept")] <- 1e-6
  }
  
  fit <- fitNbinomGLMs(object=objectNZ, modelMatrix=modelMatrix,
                       lambda=lambda, renameCols=FALSE,
                       alpha_hat=mcols(objectNZ)$dispFit,
                       betaTol=1e-4, useOptim=FALSE,
                       useQR=TRUE)
  normalizedDataNZ <- t(modelMatrix %*% t(fit$betaMatrix))

  normalizedData <- buildMatrixWithZeroRows(normalizedDataNZ, mcols(object)$allZero)

  # add back in the intercept, if finite
  if (!missing(intercept)) {
    normalizedData <- normalizedData + ifelse(infiniteIntercept, 0, intercept)
  }
  rownames(normalizedData) <- rownames(object)
  colnames(normalizedData) <- colnames(object)
  attr(normalizedData,"betaPriorVar") <- betaPriorVar
  if ("Intercept" %in% modelMatrixNames) {
    fittedInterceptNZ <- fit$betaMatrix[,which(modelMatrixNames == "Intercept"),drop=FALSE]
    fittedIntercept <- buildMatrixWithNARows(fittedInterceptNZ, mcols(object)$allZero)
    fittedIntercept[is.na(fittedIntercept)] <- -Inf
    attr(normalizedData,"intercept") <- fittedIntercept
  }
 
  #最后返回normalizedData
  normalizedData

由rlog的源代码我们知道，输入进去的dds对象需要经过负二项分布的拟合（fitNbinomGLMs），这一步主要是利用负二项回归来估计回归系数的值，而最终返回负二项分布回归的系数值

那么，经过拟合后返回fit对象，其中 fit$betaMatrix 指代的是负二项回归的回归系数，由于fitNbinomGLMs()里面部分是由C++写的，所以如果你要跑源代码，需要把对应脚本给加载下：

library(Rcpp)

Rcpp::sourceCpp('DESeq2.cpp')

我们不妨模拟运行下 rlog 里面的核心代码（运行前加载好内部函数）：

dds <- DESeqDataSetFromMatrix(data, DataFrame(condition), 
                              design=~condition)
dds <- DESeq(dds)
object = dds

samplesVector <- as.character(seq_len(ncol(object)))

samplesVector <- factor(samplesVector,levels=unique(samplesVector))

samples <- factor(c("null_level",as.character(samplesVector)),
                  levels=c("null_level",levels(samplesVector)))
modelMatrix <- stats::model.matrix.default(~samples)[-1,]
modelMatrixNames <- colnames(modelMatrix)
modelMatrixNames[modelMatrixNames == "(Intercept)"] <- "Intercept"  


# only continue on the rows with non-zero row sums
objectNZ <- object[!mcols(object)$allZero,]
stopifnot(all(!is.na(mcols(objectNZ)$dispFit)))

logCounts <- log2(counts(objectNZ,normalized=TRUE) + 0.5)
logFoldChangeMatrix <- logCounts - log2(mcols(objectNZ)$baseMean + 0.5)
logFoldChangeVector <- as.numeric(logFoldChangeMatrix)
varlogk <- 1/mcols(objectNZ)$baseMean + mcols(objectNZ)$dispFit
weights <- 1/varlogk   
betaPriorVar <- matchWeightedUpperQuantileForVariance(logFoldChangeVector, rep(weights,ncol(objectNZ)))

stopifnot(length(betaPriorVar)==1)

lambda <- 1/rep(betaPriorVar,ncol(modelMatrix))
# except for intercept which we set to wide prior


fit <- fitNbinomGLMs(object=objectNZ, modelMatrix=modelMatrix,
                     lambda=lambda, renameCols=FALSE,
                     alpha_hat=mcols(objectNZ)$dispFit,
                     betaTol=1e-4, useOptim=FALSE,
                     useQR=TRUE)
normalizedDataNZ <- t(modelMatrix %*% t(fit$betaMatrix))

normalizedData <- buildMatrixWithZeroRows(normalizedDataNZ, mcols(object)$allZero)
# add back in the intercept, if finite
rownames(normalizedData) <- rownames(object)
colnames(normalizedData) <- colnames(object)
attr(normalizedData,"betaPriorVar") <- betaPriorVar

normalizedData

尤其注意的是

 normalizedDataNZ <- t(modelMatrix %*% t(fit$betaMatrix)) 

这句命令，我们不妨打开 t(fit$betaMatrix) 看看：

t(fit$betaMatrix)

而modelMatrix如下图所示：

modelMatrix

normalizedDataNZ <- t(modelMatrix %*% t(fit$betaMatrix))

t(normalizedDataNZ)

根据矩阵的乘法原理，我们举个例子，比方说normalizedDataNZ第一行第一列的数据10.052690，它是fit$betaMatrix的intercept（第一行第一列的值）加上其第二行第一列的值，也就是每一个负二项回归的截距项加上对应的比较组的斜率（该例子的回归系数解释为sample1与control=0做负二项回归的回归系数）

小tip

当然在fitNbinomGLMs的源代码中，

 # here we transform the betaMatrix and betaSE to a log2 scale
  betaMatrix <- log2(exp(1))*betaRes$beta_mat

##计算betaRes的函数:
  betaRes <- fitBetaWrapper(ySEXP = counts(object), xSEXP = modelMatrix,
                            nfSEXP = normalizationFactors,
                            alpha_hatSEXP = alpha_hat,
                            beta_matSEXP = beta_mat,
                            lambdaSEXP = lambdaNatLogScale,
                            weightsSEXP = weights,
                            useWeightsSEXP = useWeights,
                            tolSEXP = betaTol, maxitSEXP = maxit,
                            useQRSEXP=useQR, minmuSEXP=minmu)

  # Note on deviance: the 'deviance' calculated in fitBeta() (C++)
  # is not returned in mcols(object)$deviance. instead, we calculate
  # the log likelihood below and use -2 * logLike.
  # (reason is that we have other ways of estimating beta:
  # above intercept code, and below optim code)

我们可以看到为了得到betaMatrix，软件对拟合出来的负二项回归系数（betaRes$beta_mat）做了一个乘 log2(exp(1)) 的处理，而fitBetaWrapper函数的实质就是利用负二项回归去做拟合（这一部分是C++写的），从而得到负二项回归的回归系数
C++写的这部分主要是先对row count文件做一个log标准化处理，log标准化原理具体可参考学习：StatQuest-Normalization，或者参考Introduction to DGE - ARCHIVED:
然后在对每一个基因对应于每一个sample去拟合负二项回归系数

总结

引入负二项回归来标准化，可以利用测序数据有生物学重复的优势，从一定程度上消除离群值得影响

实质上我们通过rlog标准化的raw count数据，实质上是先将row count文件做一个log标准化处理，然后对每一个基因做一个负二项回归，其响应变量为基因的表达量，决策变量为不同的分组信息。对于分组信息，相当于是sample1，2，3，4，5，6 这6个组别分别做负二项回归，得到不同的回归系数（推测这6组都是和control=0进行比较），可参照t(fit$betaMatrix)表格来看

在同一套坐标体系下：
为什么要强调在同一套坐标体系下，是因为在这样的条件下，各个基因的表达量以及拟合的负二项回归方程可以相互比较

回归系数的大小可以理解为一组分组的表达量相对于另一组分组的增长的幅度，比方说A与B两组进行比较，其负二项回归系数为5，那么从某种程度上理解可以理解为B组比A组高表达5个单位；
如果有截距项，截距项为8，A组的负二项回归系数（β1）为5，，B组的负二项回归系数（β2）为3，那么A组的基因表达量为8+5，B组的基因表达量为8+3，因此回归系数加上截距项即代表基因表达量的高低。
由于我们的决策变量是因子型变量，比方说sample 1和control=0进行建模，那么sample 1的因子为1，control=0的因子为0，假设负二项回归的回归方程为 y = βx + intercept，那么对于control组，其基因表达量为intercept；对于sample 1组，其基因表达量为β + intercept，这也就印证了上一段叙述中我们对基因表达量的讨论

因此，在与相同的control组比较的条件下，回归系数与截距项的和其实就是基因表达量的高低

因此，对于DESeq2这样的寻找差异基因的软件，将其转换为回归系数是不影响结果的