PH525x series - Running PCA and SVD in R

在PCA相关的章节最后，系列教程的作者又专门写了一章“在R中运行PCA和SVD”，使用的还是tissuesGeneExpression包中的基因表达数据。

prcomp函数与svd函数

library(tissuesGeneExpression)
data(tissuesGeneExpression)
library(rafalib)
group <- as.fumeric(tab$Tissue)

首先，一般来说，对样本做PCA分析前都需要转置一下数据（让行为样本）。prcomp()函数可以用来返回主成分以及一些其他的变量：

x <- t(e)
pc <- prcomp(x)
names(pc)
## [1] "sdev"     "rotation" "center"   "scale"    "x"

plot(pc$x[,1],pc$x[,2],col=group,main="PCA",xlab="PC1",ylab="PC2")

201912181038.png

因为prcomp()函数默认对列向量进行中心化，所以这个PCA结果与用中心化的数据（在本例中，列是基因）进行奇异值分解的结果是相同的。

##sweep()在这里的作用是令矩阵x的每一列减去该列的均值
cx <- sweep(x,2,colMeans(x),"-")
sv <- svd(cx)
names(sv)
## [1] "d" "u" "v"

plot(sv$u[,1],sv$u[,2],col=group,main = "SVD",xlab = "U1",ylab = "U2")

201912181047.png

所以在本例中，奇异值分解得到的矩阵 $\mathbf{U}$ 的列向量就相当于PCA的主成分结果 $x$ 。另外，SVD得到的矩阵 $\mathbf{V}$ 相当于prcomp()函数返回结果中的旋转矩阵（rotation）:

sv$v[1:5, 1:5]
##            [,1]      [,2]      [,3]      [,4]       [,5]
## [1,]  0.0046966 -0.013275  0.002087  0.017093  0.0006956
## [2,] -0.0021623 -0.002212  0.001543 -0.003346 -0.0034159
## [3,] -0.0030945  0.005870  0.001686  0.003890  0.0019032
## [4,] -0.0007355 -0.002002 -0.002753  0.001776  0.0192205
## [5,]  0.0010133  0.001215 -0.001561  0.003349 -0.0012380

pc$rotation[1:5, 1:5]
##                  PC1       PC2       PC3       PC4        PC5
## 1007_s_at  0.0046966 -0.013275  0.002087  0.017093  0.0006956
## 1053_at   -0.0021623 -0.002212  0.001543 -0.003346 -0.0034159
## 117_at    -0.0030945  0.005870  0.001686  0.003890  0.0019032
## 121_at    -0.0007355 -0.002002 -0.002753  0.001776  0.0192205
## 1255_g_at  0.0010133  0.001215 -0.001561  0.003349 -0.0012380

还有就是SVD分解得到的矩阵 $\mathbf{D}$ 对角线上的元素是与PCA返回的标准差成比例的，只不过prcomp()返回的标准差是样本标准差（prcomp()返回样本方差的无偏估计，所以是做了 $n/(n-1)$ 矫正的样本方差。），而 $\mathbf{D}$ 对角线上的元素却是主成分的平方和再开方得到的，并没有除以样本大小：

head(sv$d^2)
## [1] 673418 285393 182527 127667 108576  81999
head(pc$sdev^2)
## [1] 3582.0 1518.0  970.9  679.1  577.5  436.2
head(sv$d^2/(ncol(e) - 1))
## [1] 3582.0 1518.0  970.9  679.1  577.5  436.2

令方差除以方差之和，我们可以绘制主成分方差解释度的图：

plot(sv$d^2/sum(sv$d^2),xlim = c(0,15),type = "b",pch=16,
  xlab = "PCs",ylab = "variance explained")

201912181109.png

需要注意的是，在进行奇异值分解前如果没有将数据中心化的话，画出的图就会与PCA画出来的有一些不一样：

svNoCenter <- svd(x)
mypar(1,2)
plot(pc$x[,1],pc$x[,2],col=group,main="PCA",xlab = "PC1",ylab = "PC2")
points(0,0,pch=3,cex=4,lwd=4)
plot(svNoCenter$u[,1],svNoCenter$u[,2],col=group,main="SVD not centered",
     xlab="U1",ylab="U2")

201912181114.png

SVD on (genes vs samples) and (samples vs genes)

若不对数据进行中心化，无论矩阵是否转置，奇异值分解的结果都一样：

mypar(1,2)
sv2 <- svd(t(e))
plot(sv2$u[,1],sv2$u[,2],col=group,main = "samples vs genes( typical PCA )",xlab="U1",ylab="U2")
sv1 <- svd(e)
plot(sv1$v[,1],sv1$v[,2],col=group,main="genes vs samples (SVA)",xlab = "V1",ylab="V2")

201912181130.png

究竟对哪个方向进行中心化取决于你的分析目的，如果你想比较样本距离，那么行是样本且基因列被中心化，但如果你想分析基因，那么行便是基因，样本列被中心化。

阅读原文请戳