1、PCA 分析和散点图gaom今天主要跟大家演示一下简单的 PCA 分析，并且以散点图的形式将结果展示出来。 首先在进 行 PCA 分析之前，先跟大家稍微讨论下什么是 PCA 分析。PCA 分析又叫主成分分析，其实从字面上来理解我们可以发现它其实是和样品分组相关的。举个简单的例子，我们观察了某种植物的株高、叶片大小、果实大小等等多种性状，并记录每种性状对应的数值。这时候我们想看看根据这些性状信息看看我们观察的样本是否明显的分组现象。每一种性状相当于一个维度。利用 PCA 分析可以将结果投影到一个低维的向量空间（具体计算就不详述了）。类似的比如我们多个样本的表达谱数据，每个基因在各个样品的表达情况
2、就可以算作一个维度。如果大家对 PCA 算法感兴趣的话，可以自行百度，在这里就不 进行太多的描述了。毕竟今天主要是教大家怎么利用 R 进行 PCA 分析和结果展示。 还是第一步，我们先准备好我们用来分析的数据。setwd(“C:/Users/gaom/Desktop“)#打开文件所在路径，并将文件所在目 录 作 为 工作目录data-read.table(file = “test_data.txt“,header = T,sep = “t“)#读 取数据，并将首行作 为 列名dim(data)# 1 2999 13head(data)# ID_REF T01 T02 T03 T04 T05 T
3、06# 1 1007_s_at 10.198586 11.805676 10.867953 11.763660 12.072232 12.108312# 2 1053_at 9.594074 8.713108 9.247096 9.433265 9.092329 9.005518# 3 117_at 8.581763 8.603680 8.804425 8.661700 8.634979 8.606976# 4 121_at 12.022315 12.655329 12.627334 12.791390 12.961761 12.885307# 5 1255_g_at 7.228569 7.2
4、14600 7.237131 7.293417 7.276799 7.268233# 6 1294_at 8.828487 9.380277 9.297989 8.858985 8.995772 9.126825# T07 T08 T09 T10 T11 T12# 1 10.646868 10.852744 10.675898 11.137663 10.796737 11.102408# 2 9.087681 9.027208 8.965283 8.958309 9.275010 8.940965# 3 8.625838 8.577244 8.646751 8.625843 8.625164
5、8.522129# 4 13.402044 13.240126 13.088883 13.234099 13.382903 13.472223# 5 7.197440 7.262662 7.289796 7.232249 7.202364 7.306229# 6 9.002385 9.003561 9.006278 9.006721 9.018183 9.164313上述数据为从 GEO 数据库随意找的基因表达。其中第一列为基因探针号，后续几列则为 T01 到 T12 的 12 个样品对应的表达量数据，每三个 样品为一组。因为数据是拼凑的，所以这里不关注探针具体信息了。准备好数据之后我 们 就
6、开始进行 PCA 计算了。其实代码非常简单。pca- prcomp(t(data,-1), scale=T)head(pca$x)# PC1 PC2 PC3 PC4 PC5 PC6# T01 -43.457435 -44.950031 8.305571 3.210563 -7.4280481 14.818150# T02 42.067255 -19.142248 -25.574041 21.120294 -5.7930990 14.702922# T03 -2.123455 -21.512488 -11.192474 17.583006 15.2149034 -34.730308# T04 8
7、.166077 -4.774814 22.837578 -11.364128 8.4021038 -6.921738# T05 18.214073 -5.836807 18.522768 -10.941626 -0.6183613 -5.548845# T06 27.219529 -5.519328 26.649872 -11.054961 -4.1480413 5.097715# PC7 PC8 PC9 PC10 PC11 PC12# T01 -1.966342 9.2181269 -1.520882 -1.060835 3.048498 2.731227e-13# T02 5.832197
8、 8.9793018 9.386187 1.668761 1.705474 2.674666e-13# T03 -5.168168 -9.7483411 -11.570320 2.618203 -4.221456 2.738955e-13# T04 27.782986 7.5829007 9.726761 -3.391763 -21.900485 2.730871e-13# T05 7.039535 -8.9173716 -2.239005 -17.514433 29.700906 2.736544e-13# T06 -30.026232 -0.8253129 -5.207037 12.349
9、414 -8.900676 2.681674e-13summary(pca)# Importance of components:# PC1 PC2 PC3 PC4 PC5 PC6# Standard deviation 21.9980 21.7992 18.5932 16.67518 16.1346 15.16897# Proportion of Variance 0.1614 0.1585 0.1153 0.09272 0.0868 0.07672# Cumulative Proportion 0.1614 0.3198 0.4351 0.52780 0.6146 0.69133# PC7
10、 PC8 PC9 PC10 PC11# Standard deviation 14.48695 14.01978 13.4814 13.09112 12.8896# Proportion of Variance 0.06998 0.06554 0.0606 0.05714 0.0554# Cumulative Proportion 0.76131 0.82685 0.8875 0.94460 1.0000# PC12# Standard deviation 2.859e-13# Proportion of Variance 0.000e+00# Cumulative Proportion 1.
11、000e+00上述数据中，pca$x 就是后面我们画 pca 图要用的数据。而在 summary(pca)中我们看到的 Proportion of Variance 就是各个主成分的方差占所有方差的比值，即对应的贡献率。而 Cumulative Proportion 则对应的百分比累积值。从上述结果看这组数据 pca 结果并不是很好，所以应该肯定会有一些分组的结果不太好。不过我们今天主要是展示结果，就不在意这些细节了。 做完上述的计算，下面就进入我们的结果展示阶段。 首先用基本画图函数展示。plot(pca$x,1:2 )group - factor(c(rep(“A1“, 3), rep(“
12、A2“, 3),rep(“B1“, 3),rep(“B2“, 3)#这 里我 们 添加分 组 信息colour_group-rainbow(length(unique(group)#利用 rainbow 函数 选择颜 色colour-colour_groupas.numeric(factor(group)#创 建 颜 色向量colour# 1 “#FF0000FF“ “#FF0000FF“ “#FF0000FF“ “#80FF00FF“ “#80FF00FF“# 6 “#80FF00FF“ “#00FFFFFF“ “#00FFFFFF“ “#00FFFFFF“ “#8000FFFF“# 11
13、“#8000FFFF“ “#8000FFFF“plot(pca$x,1:2,col = colour ,pch = c(21,22,23,24)group)#在 plot 函数中我 们 把分 组 信息和 颜 色方案添加 进 去legend(“topleft“, legend = levels(group),col = colour_group, pch = c(21,22,23,24)#添加 legendtitle(“test“)这是我们 用基本函数 对 pca 分析结果的展示。除此外我们也可以利用 ggplot2包进行相同的图片绘制。示例如下：library(ggplot2)group2-d
14、ata.frame(group)pca_reuslt-as.data.frame(pca$x)pca_reuslt-cbind(pca_reuslt,group2)p-ggplot(pca_reuslt)+geom_point(aes(x=pca_reuslt,1,y=pca_reuslt,2,color=pca_reuslt$group,shape = pca_reuslt$group),size=5)p-p+theme(legend.title =element_blank()+labs(x=“PCA1“,y=“PCA2“)p好了，上面那些基本的结果展示我们已经结束了。下面我 们开始把这个
15、图的档次再提高一点。比如，我们画了二维的，现在我们画个三维的 PCA 结果吧。library(scatterplot3d)par(mar=c(5.1, 4.1, 4.1, 8.1), xpd=TRUE)scatterplot3d(pca_reuslt,1:3,pch=20 ,color=colour,angle=45,main = “test_3D“,cex.symbols = 2,mar=c(5.1, 4.1, 4.1, 8.1)legend(“right“,legend = group,col = colour,pch = 20,bg=“white“,xpd=TRUE,inset= -0.
16、5)#设 置位置 为 right 后，可以用 inset 来移到 legend 位置。除此之外，我们可以考虑把相同的组进行一个圈定，方便我们更好的观察结果。library(ggfortify)#使用 这 个包 时 可能要注意 R 的版本，我 刚 开始用 较 老的版本就用不了 这 个包。data2-cbind(t(data,-1),group2)autoplot(pca,data = data2,colour = group,label = TRUE,frame = TRUE)这次的数据只有 12 个样品，所以点太少了。如果点比较多的话，还可以在autoplot 函数中加上 frame.type = norm，这样就变成用椭圆形把分组的点圈起来了。 大家也可以通过 ggfortity 包的说明发现更多有趣的使用方法。
