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标准化。