学习笔记的主要内容是在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
再添加一个基因
刚刚是两个基因，现在再加一个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)
通过上面的演示，已经基本了解PCA的作用了，通过PCA分析能将不同基因在不同样本中的表达量分成几类，接下来，用简单的例子来演示流程。
