PCA的实现流程
使用上面创建的data_3数据来进行后续操作。首先生成表达矩阵，包含3个基因在200个样本中的表达情况。
> kable(headTail(data_3),booktabs=T,caption = "Expression 3Gene in 200 samples")
Table: Expression 3Gene in 200 samples
| |Ge_1 |Ge_2 |Ge_3 |Group |
|:----|:-----|:----|:----|:-----|
|a1 |4.77 |5.71 |5.61 |a |
|a2 |4.13 |5.65 |6.38 |a |
|a3 |4.15 |6.38 |6.47 |a |
|a4 |4.04 |5.88 |5.82 |a |
|... |... |... |... |NA |
|b97 |18.93 |6.06 |3.57 |b |
|b98 |18.06 |6.29 |4.37 |b |
|b99 |18.74 |5.78 |4.18 |b |
|b100 |19.52 |5.87 |4.82 |b |
# 对数据进行标准化处理
> data_3_cs <- scale(data_3[,1:3],center = T,scale = T)
> kable(headTail(data_3_cs),booktabs=T,caption = "norm Expression 3 gene in 200 samples")
上面的代码是对数据进行标准化和中心化处理（使数据的差异变化幅度在同一水平），将数据转化为均值为0且标准差为1的新数据集。
Table: norm Expression 3 gene in 200 samples
| |Ge_1 |Ge_2 |Ge_3 |
|:----|:-----|:-----|:-----|
|a1 |-0.89 |-1 |-0.87 |
|a2 |-0.98 |-1.22 |-0.7 |
|a3 |-0.97 |1.41 |-0.68 |
|a4 |-0.99 |-0.37 |-0.82 |
|... |... |... |... |
|b97 |0.99 |0.25 |-1.32 |
|b98 |0.88 |1.08 |-1.14 |
|b99 |0.97 |-0.73 |-1.18 |
|b100 |1.07 |-0.44 |-1.04 |
> data_3_cs_cov <- cov(data_3_cs)
> kable(data_3_cs_cov,booktabs=T,
+ caption = "cov for 3 gene in 200 samples")
上面的代码生成协方差矩阵，计算3个基因在200个样本中表达数据的协方差。
Table: cov for 3 gene in 200 samples
| | Ge_1| Ge_2| Ge_3|
|:----|----------:|----------:|----------:|
|Ge_1 | 1.0000000| -0.0808226| -0.1181946|
|Ge_2 | -0.0808226| 1.0000000| -0.0106916|
|Ge_3 | -0.1181946| -0.0106916| 1.0000000|
> data_3_cs_cov_e <- eigen(data_3_cs_cov)
#求解特征值和特征向量
> data_3_cs_cov_e$values #特征值
> [1] 1.1383477 1.0099558 0.8516964
> data_3_cs_cov_e$vectors #特征向量
> [,1] [,2] [,3]
> [1,] 0.7189945 0.02734216 -0.6944778
> [2,] -0.3748044 -0.82622441 -0.4205650
> [3,] -0.5852936 0.56267720 -0.5838028
上面的代码得到特征值和特征变量，下面的代码用于产生新矩阵。
> pc_select <- 3
> label <- paste0("PC",c(1:pc_select))
> data_3_n <- data_3_cs %*% data_3_cs_cov_e$vectors[,1:pc_select] #%*%表示矩阵相乘
> colnames(data_3_n) <- label
> kable(headTail(data_3_n),booktabs=T,
+ caption = "PCA gene matrix for 3 gene in 200 samples")
Table: PCA gene matrix for 3 gene in 200 samples
| |PC1 |PC2 |PC3 |
|:----|:-----|:-----|:-----|
|a1 |0.24 |0.31 |1.55 |
|a2 |0.16 |0.59 |1.6 |
|a3 |-0.83 |-1.57 |0.48 |
|a4 |-0.09 |-0.18 |1.32 |
|... |... |... |... |
|b97 |1.39 |-0.92 |-0.02 |
|b98 |0.89 |-1.51 |-0.4 |
|b99 |1.66 |-0.03 |0.33 |
|b100 |1.54 |-0.19 |0.05 |接下来，比较两种方式对样本的聚类差异情况，设置工作区同时输出两个图，并使用scatterplot3d进行绘图。
colorl <- c("#E38F92","#97B6E1")
colors <- colorl[as.numeric(data_3$Group)]
par(mfrow=c(1,2)) #图片输出区为一行两图的布局
scatterplot3d(data_3[,1:3],color = colors,
angle=45,pch=16,main="before data")
# 生成图例legend("top",legend = levels(data_3$Group),col=colorl,pch=16,xpd=T,horiz = T)
scatterplot3d(data_3_n,color=colors,angle = 45,pch=16,
main="after data")
通过对比上图，可以发现两种数据处理方式形成的样品分组情况不同，在处理后数据右图中，样本的分散程度更大，笔者的理解是其变化特征显示的更广泛，相比左图能够读取更多信息，处理后效果更好（可能是因为此时变量间非线性相关）。
利用prcomp进行PCA分析
pca_data_3 <- prcomp(data_3[,1:3],center=T,scale=T)
str(pca_data_3)
上面的代码对data_3数据进行处理，得到新数据，接着查看一下pca_data_3的数据信息摘要。
List of 5
$ sdev : num [1:3] 1.067 1.005 0.923
$ rotation: num [1:3, 1:3] -0.719 0.3748 0.5853 0.0273 -0.8262 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : chr [1:3] "Ge_1" "Ge_2" "Ge_3"
.. ..$ : chr [1:3] "PC1" "PC2" "PC3"
$ center : Named num [1:3] 11.47 5.99 9.55
..- attr(*, "names")= chr [1:3] "Ge_1" "Ge_2" "Ge_3"
$ scale : Named num [1:3] 7.52 0.277 4.548
..- attr(*, "names")= chr [1:3] "Ge_1" "Ge_2" "Ge_3"
$ x : num [1:200, 1:3] -0.2399 -0.1632 0.833 0.0905 0.3406 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : chr [1:200] "a1" "a2" "a3" "a4" ...
.. ..$ : chr [1:3] "PC1" "PC2" "PC3"
- attr(*, "class")= chr "prcomp"
生成新的数据包含五个变量，按照之前的方法对其进行处理。
data_pca_n <- pca_data_3$x
kable(headTail(data_pca_n),booktabs=T,
caption = "PCA gene matrix")得到prcomp方式的基因表达矩阵，此时存在三个主成分（PC1、2、3）。
Table: PCA gene matrix
| |PC1 |PC2 |PC3 |
|:----|:-----|:-----|:-----|
|a1 |-0.24 |0.31 |1.55 |
|a2 |-0.16 |0.59 |1.6 |
|a3 |0.83 |-1.57 |0.48 |
|a4 |0.09 |-0.18 |1.32 |
|... |... |... |... |
|b97 |-1.39 |-0.92 |-0.02 |
|b98 |-0.89 |-1.51 |-0.4 |
|b99 |-1.66 |-0.03 |0.33 |
|b100 |-1.54 |-0.19 |0.05 |
# 查看特征向量
> pca_data_3$rotation
PC1 PC2 PC3
Ge_1 -0.7189945 0.02734216 -0.6944778
Ge_2 0.3748044 -0.82622441 -0.4205650
Ge_3 0.5852936 0.56267720 -0.5838028接下来，比较两种方式实现PCA分析的结果差异，左图是手动方式，右图是利用prcomp方式，可以看出两种结果具有差异性。
scatterplot3d(data_3_n,color=colors,angle=45,pch=16,
main="PCA by steps")
scatterplot3d(data_pca_n,color=colors,angle=45,pch=16,
main="PCA by prcomp")
创建PCA计算函数
在R语言中自定义一个函数ct_PCA，用于计算处理PCA数据（参数设置对原始数据进行标准化和中心化）
ct_PCA <- function(data,center=T,scale=T){
data_norm <- scale(data, center=center, scale=scale)
data_norm_cov <- crossprod(as.matrix(data_norm)) / (nrow(data_norm)-1)
data_eigen <- eigen(data_norm_cov)
rotation <- data_eigen$vectors
label <- paste0('PC', c(1:ncol(rotation)))
colnames(rotation) <- label
sdev <- sqrt(data_eigen$values)
data_new <- data_norm %*% rotation
colnames(data_new) <- label
ct_pca <- list('rotation'=rotation, 'x'=data_new, 'sdev'=sdev)
return(ct_pca)
}标准化scale操作是指将数据的差异程度相对化，消除固有差异幅度的影响，从同一衡量标准下判断数据的差异性，接下来，分别演示不经过标准化处理和进行标准化处理的结果。
data_pca_noscale_step <- ct_PCA(data_3[,1:3],center=T,scale = F)
#只中心化，不标准化
data_pca_noscale_step$rotation #查看特征向量
PC1 PC2 PC3
[1,] 0.993858995 -0.110611181 -0.003076602
[2,] -0.002918535 0.001590917 -0.999994476
[3,] -0.110615464 -0.993862483 -0.001258325
data_pca_noscale_pc <- data_pca_noscale_step$x利用刚才生成的四种数据，生成四个不同类型的结果图：
par(mfrow=c(2,2)) #设置输出区为2行2列排版，同时输出4副图
scatterplot3d(data_3[,c(1,3,2)],color=colors,
angle=45,pch=16,main="ori plot")
scatterplot3d(data_pca_noscale_pc,color=colors,
angle=45,pch=16,main="PCA noscale")
scatterplot3d(data_3_cs[,c(1,3,2)],color=colors,
angle=45,pch=16,main="ori plot(scale)")
scatterplot3d(data_3_n,color=colors,
angle=45,pch=16,main="PCA scale")依次生成4副图，可以看出上面两张图（没有scale标准化）的分布比较秘籍，而经过scale处理之后数据的分散程度更高（下面两张图），说明标准化处理后数据的相对变化幅度信息被保留，差异细节更清晰，这也是PCA分析的目的所在。
本文中所有代码已整理打包，下载链接：
https://down.jewin.love/?f=/Rscript/PCA.R
参考资料：http://www.ehbio.com
END
© 素材来源于网络，侵权请联系后台删除
点击查看往期推荐：
ggplot2能干啥？漫画式介绍这款绘图神器
如何绘制热图？ggplot2入门笔记
绘制热图时看不出颜色差异？
