R语言作图之PCA作图和散点图

PCA分析和散点图今天主要跟大家演示一下简单的PCA分析，并且以散点图的形式将结果展示出来。首先在进行PCA分析之前，先跟大家稍微讨论下什么是PCA分析。PCA分析又叫主成分分析，其实从字面上来理解我们可以发现它其实是和样品分组相关的。举个简单的例子，我们观察了某种植物的株高、叶片大小、果实大小等等多种性状，并记录每种性状对应的数值。这时候我们想看看根据这些性状信息看看我们观察的样本是否明显的分组现象。每一种性状相当于一个维度。利用PCA分析可以将结果投影到一个低维的向量空间(具体计算就不详述了)。类似的比如我们多个样本的表达谱数据，每个基因在各个样品的表达情况就可以算作一个维度。如果大家对PCA算法感兴趣的话，可以自行百度，在这里就不进行太多的描述了。毕竟今天主要是教大家怎么利用R进行PCA分析和结果展示。还是第一步，我们先准备好我们用来分析的数据。setwd("C:/Users/gaom/Desktop")#打开文件所在路径,并将文件所在目录作为工作目录data<-read.table(file="test_data.txt",header=T,sep="\t")#读取数据，并将首行作为列名一dim(data)#[1]299913head(data)#ID_REFT01T02T03T04T05T06#11007_s_at10.19858611.80567610.86795311.76366012.07223212.108312#21053_at9.5940748.7131089.2470969.4332659.0923299.005518#3117_at8.5817638.6036808.8044258.6617008.6349798.606976#4121_at12.02231512.65532912.62733412.79139012.96176112.885307#51255_g_at7.2285697.2146007.2371317.2934177.2767997.268233#61294_at8.8284879.3802779.2979898.8589858.9957729.126825#T07T08T09T10T11T12#110.64686810.85274410.67589811.13766310.79673711.102408#29.0876819.0272088.9652838.9583099.2750108.940965#38.6258388.5772448.6467518.6258438.6251648.522129##413.40204413.24012613.08888313.23409913.38290313.472223##57.1974407.2626627.2897967.2322497.2023647.306229##69.0023859.0035619.0062789.0067219.0181839.164313上述数据为从GEO数据库随意找的基因表达。其中第一列为基因探针号，后续几列则为T01到T12的12个样品对应的表达量数据，每三个样品为一组。因为数据是拼凑的，所以这里不关注探针具体信息了。准备好数据之后我们就开始进行PCA计算了。其实代码非常简单。pca<-prcomp(t(data[,-1]),scale=T)head(pca$x)#PC1PC2PC3PC4PC5PC6#T01-43.457435-44.9500318.3055713.210563-7.428048114.818150#T0242.067255-19.142248-25.57404121.120294-5.793099014.702922#T03-2.123455-21.512488-11.19247417.58300615.2149034-34.730308#T048.166077-4.77481422.837578-11.3641288.4021038-6.921738#T0518.214073-5.83680718.522768-10.941626-0.6183613-5.548845#T0627.219529-5.51932826.649872-11.054961-4.14804135.097715#PC7PC8PC9PC10PC11PC12#T01-1.9663429.2181269-1.520882-1.0608353.0484982.731227e-13#T025.8321978.97930189.3861871.6687611.7054742.674666e-13#T03-5.168168-9.7483411-11.5703202.618203-4.2214562.738955e-13#T0427.7829867.58290079.726761-3.391763-21.9004852.730871e-13#T057.039535-8.9173716-2.239005-17.51443329.7009062.736544e-13#T06-30.026232-0.8253129-5.20703712.349414-8.9006762.681674e-13summary(pca)#Importanceofcomponents:#PC1PC2PC3PC4PC5PC6#Standarddeviation21.998021.799218.593216.6751816.134615.16897#ProportionofVariance0.16140.15850.11530.092720.08680.07672

#CumulativeProportion0.16140.31980.43510.527800.61460.69133#PC7PC8PC9PC10PC11#Standarddeviation14.4869514.0197813.481413.0911212.8896#ProportionofVariance0.069980.065540.06060.057140.0554#CumulativeProportion0.761310.826850.88750.944601.0000#PC12#Standarddeviation2.859e-13#ProportionofVariance0.000e+00#CumulativeProportion1.000e+00上述数据中，pca$x就是后面我们画pca图要用的数据。而在summary(pca)中我们看到的ProportionofVariance就是各个主成分的方差占所有方差的比加即对应的贡献率。而CumulativeProportion则对应的百分比累积值。从上述结果看这组数据pca结果并不是很好，所以应该肯定会有一些分组的结果不太好。不过我们今天主要是展示结果，就不在意这些细节了。做完上述的计算，下面就进入我们的结果展示阶段。首先用基本画图函数展示。plot(pca$x[,1:2])CMO—on_昌-Io-402040T--402040PC1group<-factor(c(rep("A1",3),rep("A2",3),rep("B1",3),rep("B2",3)))这里我们添加分组信息colour_group<-rainbow(length(unique(group)))#利]用rainbow函数选择颜色colour<-colour_group[as.numeric(factor(group))]#创建颜色向量colour

#[1]"#FF0000FF""#FF0000FF""#FF0000FF""#80FF00FF""#80FF00FF"#[6]"#80FF00FF""#00FFFFFF""#00FFFFFF""#00FFFFFF""#8000FFFF”#[11]"#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(2122,23,24))#添加legendtitle("test")testA1A2△△B2CMO—on_W-o-402040-402040PC1这是我们用基本函数对pca分析结果的展示。除此外我们也可以利用ggplot2包进行相同的图片绘制。示例如下：library(ggplot2)group2<-data.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")p20-*A1

▲A2

■B1B2*250PCA125好了，上面那些基本的结果展示我们已经结束了。下面我们开始把这个图的档次再提高一点。比如，我们画了二维的，现在我们画个三维的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=st_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=inset=-0.5)#设置位置为right后，可以用inset来移到legend位置。"teTRUEtest3DA1A1A1A2A2A2B1B1B1B2B2B2PC1除此之外，我们可以考虑把相同的组进行一个圈定，方便我们更好的观察结果。library(ggfortify)#使用这个包时可能要注意R的版本，我刚开始用较老的版本就用不了这个包。