共计 9112 个字符，预计需要花费 23 分钟才能阅读完成。
学习笔记的次要内容是在 R 语言中利用 ggplot2 进行 PCA 剖析和绘图，包含简略剖析与操作流程，比照不同形式失去的后果差别，提供脚本代码供练习.
PCA 剖析的原理
在解决基因差别表白数据时，有时候须要剖析其中因素的影响最大，判断后果的关系，这个时候能够用 PCA 分析法，之前发过一篇 PCA 剖析的简介和数学原理解析，如果有趣味点击这里查看，明天的笔记次要围绕实际操作过程进行分享。笔者学习时参考易汉博的教程，感觉这个教程挺好的，举荐给大家，也能够在学习过程中一起交换。
PCA 剖析示例
创立演示数据
count <- 100 #设置变量个数为 100
Ge_1a <- rnorm(count,4,0.6) #生成 100 个遵从均值为 4 标准差为 0.6 正态分布的数字
Ge_1b <- rnorm(count,19,0.4)
gro_a <- rep('a',count) #生成 100 个 a，代表 a 组
gro_b <- rep('b',count)演示数据为 Ge_1 的表白量（每个基因包含两组类型的值各 100 个，且两个组的表白量有差别），接下来创立依据数据创立矩阵，设置样本的名称标签，增加新列R，并生成一个表格输入基因在 200 个样本中的表白量，每一行为一个样品，每一列为基因的表白值。
c_data <- data.frame(Ge_1=c(Ge_1a,Ge_1b),
Group=c(gro_a,gro_b))
label <- c(paste0(gro_a,1:count),paste0(gro_b,1:count))
row.names(c_data) <- label
c_data$R <- rep(0,count*2)
kable(headTail(c_data),booktabs=TRUE,
caption="Expr Profile For Ge_1 in 200 samples")生成了 200 行 3 列的表格数据，后果如下：
Table: Expr Profile For Ge_1 in 200 samples
| |Ge_1 |Group |R |
|:----|:-----|:-----|:---|
|a1 |4.77 |a |0 |
|a2 |4.13 |a |0 |
|a3 |4.15 |a |0 |
|a4 |4.04 |a |0 |
|... |... |NA |... |
|b97 |18.93 |b |0 |
|b98 |18.06 |b |0 |
|b99 |18.74 |b |0 |
|b100 |19.52 |b |0 |加载 R 包
library(knitr)
library(psych)
library(reshape2)
library(ggplot2)
library(ggbeeswarm)
library(scatterplot3d)
library(useful)
library(ggfortify)须要加载上述 R 包，如果没有请先装置后载入 R 包。
绘制图像
p <- ggplot(c_data,aes(Ge_1,R)) + geom_quasirandom(aes(color=factor(Group))) +theme(legend.position = c(0.5,0.8)) +
theme(legend.title = element_blank()) +
scale_fill_discrete(name="Group") +
theme(axis.line.y=element_blank(),
axis.text.y=element_blank(),
axis.ticks.y=element_blank(),
axis.title.y=element_blank()) +
ylim(-0.5,5) + xlim(0,25)
p 利用 ggplot 函数进行绘图，发现 200 个样本在 Ge1 基因表白量上分为了两类（起因是刚刚生成数据时分成了两个不同类型的组，表白量存在差别）
增加一个基因
刚刚是只有 Ge1 的状况，接下来再创立一个 Ge_2（办法和刚刚相似），看看两个基因时状况会产生什么变动？
创立一组均值为 6 标准差为 0.3 的正态分布随机数据，并设置行名构建矩阵，输入表白矩阵。须要留神的是：Ge_2的表白量保持稳定（a 组和 b 组的表白程度相当），不像 Ge_1 存在表白量差别。
> count <- 100
> Ge_2a <- rnorm(count,6,0.3)
> Ge_2b <- rnorm(count,6,0.3)
> c2_data <- data.frame(Ge_1=c(Ge_1a,Ge_1b),Ge_2=c(Ge_2a,Ge_2b),
+ Group=c(gro_a,gro_b))
> row.names(c2_data) <- label
> kable(headTail(c2_data),booktabs=T,
+ caption="Expression for Ge_1 and Ge_2 in 200 samples")
Table: Expression for Ge_1 and Ge_2 in 200 samples
| |Ge_1 |Ge_2 |Group |
|:----|:-----|:----|:-----|
|a1 |4.77 |5.71 |a |
|a2 |4.13 |5.65 |a |
|a3 |4.15 |6.38 |a |
|a4 |4.04 |5.88 |a |
|... |... |... |NA |
|b97 |18.93 |6.06 |b |
|b98 |18.06 |6.29 |b |
|b99 |18.74 |5.78 |b |
|b100 |19.52 |5.87 |b |利用 ggplot 函数作图，数据为 c2data，此时能显示出Ge1 和Ge2的散布状况，能够看出在 Ge_1(x 轴) 上分成了两类，而 Ge_2 上分类趋势很小（因为 Ge_2 自身就没什么差别分组）
p <- ggplot(c2_data,aes(Ge_1,Ge_2)) +
geom_point(aes(color=factor(Group))) +
theme(legend.position = c(0.5,0.8)) +
theme(legend.title = element_blank()) +
ylim(0,10)+xlim(0,25)
p<img src=”https://jewin.oss-cn-hangzhou.aliyuncs.com/pic/image-20221102202914671.png” style=”zoom:80%;” />
再增加一个基因
刚刚是两个基因，当初再加一个Ge_3，这个基因的表白量差别设置的更大一些，设置该基因分成两个组，而且每个组的表白量也存在两种类型，所以这个基因对样本分类的作用更大。
> count <- 100
> Ge_3a <- c(rnorm(count/2,6,0.4),rnorm(count/2,14,0.3))
> Ge_3b <- c(rnorm(count/2,14,0.3),rnorm(count/2,4,0.4))
> data_3 <- data.frame(Ge_1=c(Ge_1a,Ge_1b),
+ Ge_2=c(Ge_2a,Ge_2b),
+ Ge_3=c(Ge_3a,Ge_3b),
+ Group=c(gro_a,gro_b))
> data_3 <- as.data.frame(data_3)
> data_3$Group <- as.factor(data_3$Group)
> row.names(data_3) <- label
>
> kable(headTail(data_3),booktabs=T,caption = "Expression 3 Genes in 200 samples")
Table: Expression 3 Genes 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 |生成一组色彩变量，用于辨别不同类别。每个数据向底面做垂直投影，能够看出在 x 轴方向（Ge_1）和 z 轴（Ge_3）上投影时在不同地位分成两类，而在 y 轴（Ge_2）上投影位于同一区域，所以能够看出 Ge_2 对样本分类的贡献度最小。
colorl <- c("#E19F90", "#96B4E9")
colors <- colorl[as.numeric(data_3$Group)]
scatterplot3d(data_3[,1:3],color=colors,xlim=c(0,24),
ylim=c(0,24),zlim=c(0,24),type="h",
angle=45,pch=16)
legend("top",legend=levels(data_3$Group),col=colorl,
pch=16,xpd=T,horiz=T)<img src=”https://jewin.oss-cn-hangzhou.aliyuncs.com/pic/image-20221102204238292.png” style=”zoom:50%;” />
通过下面的演示，曾经根本理解 PCA 的作用了，通过 PCA 剖析能将不同基因在不同样本中的表白量分成几类，接下来，用简略的例子来演示流程。
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=/R…
参考资料：http://www.ehbio.com
本文由 mdnice 多平台公布
