R tips：细究FactoMineR的z-score标准化细节

R中的做主成分分析(PCA)有很多函数，如R自带的prcomp、princomp函数以及FactoMineR包中PCA函数，要论分析简单和出图优雅还是FactoMineR的PCA函数（绘图可以搭配factoextra包）。

比如对于R内建数据集iris，可以如下进行PCA分析绘图：

local({
  dat <- iris[,-5]
  group <- iris$Species

  dat %>% 
    FactoMineR::PCA(graph = F) %>%
    factoextra::fviz_pca_ind(habillage = group,
                             addEllipses = TRUE,
                             label = "none")
})

绘图如下：

有关于FactoMineR的PCA函数所进行的降维分析，现在关注它的一个问题：它进行PCA分析时是否对数据进行了标准化处理以及如何设置进行标准化处理？

先直接说结论：FactoMineR的PCA函数默认进行z-score标准化，永远进行均值中心化处理。

在进行PCA、聚类等数据分析时，先进行数据标准化处理往往是必不可少的的步骤，它们的作用在于调整好各个特征（基因）的权重。一般情况下更关注这些特征在不同样本间的变化趋势，而不是绝对表达量。举例来说，对一个转录组表达矩阵进行降维分析，那么我们一般并不希望高表达基因和低表达基因的分析权重有不同：假如CD45整体高表达在1000左右而CD3低表达在100左右，但是一般情况下并不希望CD45的作用要高于CD3的作用，其他亦然，这个时候就可以通过标准化处理，比如z-score处理，将他们调整为统一的数据尺度。

当然如果是极低表达的基因，甚至他们就是一些实验噪声，那么可以在进行PCA分析前将这些基因过滤掉。 如果确实需要将表达量的高低表达水平纳入到分析中，可以进行log处理，在保留数据高低趋势的情况下也尽可能的收缩数据范围。

在FactoMineR中是默认进行z-score处理的，z-score处理就是将特征（基因）减去均值，除以标准差。FactoMineR的PCA函数的帮助文档比较隐晦，只有一个scale.unit参数是用于是否将数据调整为单位方差，看起来和z-score有点关系，其他就不太清晰了，所以我们直接看源码（本文只关注它的z-score的处理过程）：

FactoMineR::PCA
# function (X, scale.unit = TRUE, ncp = 5, ind.sup = NULL, quanti.sup = NULL, 
#           quali.sup = NULL, row.w = NULL, col.w = NULL, graph = TRUE, 
#           axes = c(1, 2)) 
# {
#   moy.ptab <- function(V, poids) {
#     as.vector(crossprod(poids/sum(poids), as.matrix(V)))
#   }
#   ec.tab <- function(V, poids) {
#     ecart.type <- sqrt(as.vector(crossprod(poids/sum(poids), as.matrix(V^2))))
#     ecart.type[ecart.type <= 1e-16] <- 1
#     return(ecart.type)
#   }
#   fct.eta2 <- function(vec, x, weights) {
#     VB <- function(xx) {
#       return(sum((colSums((tt * xx) * weights)^2)/ni))
#     }
#     tt <- tab.disjonctif(vec)
#     ni <- colSums(tt * weights)
#     unlist(lapply(as.data.frame(x), VB))/colSums(x * x * weights)
#   }
#   X <- as.data.frame(X)
#  ...省略...

函数源码的开头就是先定义了三个函数，其中前两个就是用于计算列均值和标准差的，特别是标准差的计算，为了保证z-score标准化不出现NA值，FactoMineR::PCA的处理其实非常聪明（见后述）。

后面紧接着的一条命令就是将输入数据X先转换为数据框，这里说一下X其实就是表达矩阵，但是不同于常规的"列是样本行是基因"的表达矩阵，它其实是"行是样本列是基因"（重要），所以进行PCA分析时，往往需要转置，且可以知道z-score是对基因进行标准化，而不是样本。

先接着往下看，跳过大概小三十行就可以看到如下代码，这里开始就是进行z-score处理的源码了：

# ...省略...
# if (is.null(row.w)) 
#   row.w <- rep(1, nrow(X))
# row.w.init <- row.w
# row.w <- row.w/sum(row.w)

# if (is.null(col.w)) 
#   col.w <- rep(1, ncol(X))

row.w和col.w是对行与列所分配的权重，默认不分配权重，所以就全是1，代码中是使用rep函数创建的一个全是1的向量：rep(1, nrow(X))与rep(1, ncol(X))。另外这里还将row.w转换为了比例值（除以所有权重和），比如有100个样本，那么默认情况下每个样本的权重就是1/100。

减去均值

接下来就是计算每一个基因的均值，然后每个基因各自减去自己的均值。

# centre <- moy.ptab(X, row.w)
# data <- X
# X <- t(t(as.matrix(X)) - centre)

#<><><><>moy.ptab<><><><>
#moy.ptab <- function(V, poids) {
#  as.vector(crossprod(poids/sum(poids), as.matrix(V)))
#}
#<><><><><><><><><><><><>

moy.tab就是FactoMineR的PCA函数开头创建的一个函数，用于计算标准差，此处调用时传入了表达矩阵X和行权重row.w。

而在moy.tab函数内部，是使用矩阵乘法crossprod实现的基因均值计算，另外由于row.w已经是比例值了，其实这个函数的函数体在这里可以简化为：as.vector(crossprod(poids, as.matrix(V)))，因为sum(poids)的值就是1。矩阵乘法代表表达矩阵的每一列都是和这个行权重的线性组合，其结果就是一个均值。

后面的代码就是将原来的表达矩阵减去这个均值向量即可，之所以要转置是因为R中的矩阵默认是进行列方向的自动对齐。

可以发现这个过程是没有参数控制的，所以FactoMineR的PCA函数一定会进行均值中心化处理。

除以标准差

再往下就是将每一个基因的标准差调为1，也就是先计算每一列的标准差，再将每一列除以各自的标准差。

#  ...省略...
# if (scale.unit) {  
#   ecart.type <- ec.tab(X, row.w)
#   X <- t(t(X)/ecart.type)
# }
# else ecart.type <- rep(1, length(centre))
#  ...省略...

#<><><><>ec.tab<><><><>
#ec.tab <- function(V, poids) {
#  ecart.type <- sqrt(as.vector(crossprod(poids/sum(poids), as.matrix(V^2))))
#  ecart.type[ecart.type <= 1e-16] <- 1
#  return(ecart.type)
#}
#<><><><><><><><><><><>

ec.tab函数是FactoMineR的PCA函数开头创建的第二个函数，用于计算标准差，此处调用同样是传入表达矩阵X和行权重row.w。

在ec.tab函数内部，计算标准差的是（1）先计算方差：crossprod(poids/sum(poids), as.matrix(V^2)，同样的sum(poids)的值也是1可以省略，对表达矩阵的每一列的平方值进行线性组合就是方差值：

；（2）对方差开方即可获得标准差。

下面就是FactoMineR处理的比较稳健地方：

它将小于1e-16的值设为1，这是为了解决有一些基因的标准差是0的问题，如果表达量除以0的话，在R中会出现Inf或者NaN值，后续计算就容易出现报错，所以这里将非常小的值直接调整为1了。

那么调整为1有没有道理？ 方差非常小甚至是0的情况下，那么就说明表达数据就是近乎一样的值，而经过中心化以后，这些值其实近乎是0。当都除以1的时候其实还是一群近乎0的值，这种值在聚类也不会起到太多的作用，所以调为1是比较合理的，就是不做任何处理的意思。 像这种基因由于在数据分析中起不到太大作用，其实也是可以直接丢弃的。 另外，R中scale函数也是进行的z-score标准化，如果不注意这个scale函数就会引入Inf或者NaN值，然后就可能是代码莫名报错。

计算好了标准差后，同样的道理需要先将原始表达矩阵转置，将每一列除以各自的标准差即可：X <- t(t(X)/ecart.type)。

另外看源码就可以知道，这个除以标准差的过程是可以控制的，由参数scale.unit控制，默认是TRUE，所以FactoMineR的PCA函数默认进行z-score标准化。