Seurat的normalization和scaling

Seurat的分析流程有两步, 对数据的normalization和scaling. 两种的作用不同，前者是为了处理每个细胞的总count不同的问题，而后者则是让每个基因的表达量的均值为0，方差为1.

normlization对应的函数是NormalizeData，通过数据进行一些列变换，消除文库大小的影响。 它有三种方法, LogNormalize, CLR, RC

默认方法方式是LogNormalize， 即对于每个细胞，将每个基因的count除以总数，然后乘以一个scale.factor, 之后以自然对数进行转换。为了提高效率，Seurat编写了C++代码用于加速

// [[Rcpp::export]]
Eigen::SparseMatrix<double> LogNorm(Eigen::SparseMatrix<double> data, int scale_factor, bool display_progress = true){
  Progress p(data.outerSize(), display_progress);
  Eigen::VectorXd colSums = data.transpose() * Eigen::VectorXd::Ones(data.rows());
  for (int k=0; k < data.outerSize(); ++k){
    p.increment();
    for (Eigen::SparseMatrix<double>::InnerIterator it(data, k); it; ++it){
      it.valueRef() = log1p(double(it.value()) / colSums[k] * scale_factor);
    }
  }
  return data;
}

也就是说，如果你不追求效率，我们是可以用纯R代码实现。

mat <- matrix(data = rbinom(n = 25, size = 5, prob = 0.2), nrow = 5)
mat_norm <- LogNormalize(data = mat)
# LogNormalize等价于
log1p(t(t(mat) / colSums(mat)) * 10000)

Scale这一步对应的函数是ScaleData， 在它处理之后，使得每个基因在所有样本的均值是0，而方差是1。这让每个基因在下游分析中的具有相同的权重，使得高表达基因不那么显著。

而如果需要移除数据集中不需要的变异来源(unwanted sources of variation), ScaleData需要设置额外的参数vars.to.regress，但是作者更加推荐使用SCTransform。

尽管ScaleData对应的源代码非常的长，但是和数据处理相关的就是如下这几条

# Currently, RegressOutMatrix will do nothing if latent.data = NULL
data.scale <- scale.function(
  mat = object[features[block[1]:block[2]], split.cells[[group]], drop = FALSE],
  scale = do.scale,
  center = do.center,
  scale_max = scale.max,
  display_progress = FALSE
)

其中scale.function是根据数据类型来决定是FastRowScale还是FastSparseRowScale, 而这两个代码Seurat也是通过C++进行提速了。因此，如果不追求效率，还是可以用纯R代码实现。

pbmc_small <- ScaleData(pbmc_small)
GetAssayData(pbmc_small, "scale.data")[c("IGLL5"),1:2]
# 等价于
scaled_data <- t(scale(t(GetAssayData(pbmc_small, "data"))))
scaled_data[c("IGLL5"),1:2]

默认R语言的scale是按列处理，对于行为基因，列为样本的数据，可以直接套用。但是Seurat的输入数据是行为基因，列为样本，那么就需要按行scale。换句话说，如果你需要对行进行scale，那么你可以通过Seurat:::FastRowScale()的方式调用Seurat写的FastRowScale函数。此外Seurat编写了大量的Fast*函数，都可以尝试用在代码中。