在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
所以在本例中，奇异值分解得到的矩阵的列向量就相当于PCA的主成分结果
。另外，SVD得到的矩阵
相当于
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分解得到的矩阵对角线上的元素是与PCA返回的标准差成比例的，只不过
prcomp()返回的标准差是样本标准差（prcomp()返回样本方差的无偏估计，所以是做了矫正的样本方差。），而
对角线上的元素却是主成分的平方和再开方得到的，并没有除以样本大小：
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
究竟对哪个方向进行中心化取决于你的分析目的，如果你想比较样本距离，那么行是样本且基因列被中心化，但如果你想分析基因，那么行便是基因，样本列被中心化。
阅读原文请戳
网友评论