ggplot2一张图的诞生（2）——PCA

gzh：BBio

开始之前，说下Rstudio查看帮助文档的方法：

#方法1
#右下角Help菜单栏直接搜索函数或者包

#方法2
?getwd #命令行处键入?+函数，并运行，右下角Help菜单栏会自动弹出帮助文档

读取表达量文件

getwd()  #返回当前工作路径，"C:/Users/jiutun/Documents"
fpkm<-read.table("/Users/jiutun/Documents/fpkm.txt",sep="\t",row.names=1,header=T)

group<-read.table("C:/Users/jiutun/Documents/group.txt",sep="\t",header=T)

image image

分析主成分

pca<-prcomp(fpkm,scale=T)
pcadata<-pca$rotation

pca1<-round(summary(pca)$importance[2,1]*100,1)
pca2<-round(summary(pca)$importance[2,2]*100,1)
pca3<-round(summary(pca)$importance[2,3]*100,1)

xlab<-paste("PCA1(",pca1,"%)",sep="")
ylab<-paste("PCA2(",pca2,"%)",sep="")
zlab<-paste("PCA3(",pca3,"%)",sep="")

绘图

• ggplot2绘制PCA图

library(ggplot2)
library(dplyr)
#合并分组信息
pcadata<-data.frame(sample=rownames(pcadata),pcadata) %>% left_join(group,by="sample")
rownames(pcadata)=pcadata[,1]

p1<-ggplot(pcadata,aes(PC1,PC2,fill=group,color=group))+
    geom_point()+
    geom_text(aes(label=sample),vjust=-0.5,size=3)+
    labs(x=xlab,y=ylab)

image

p1+stat_ellipse(level=0.95,show.legend=F,geom="polygon",alpha=0.2)

image

• scatterplot3d绘制PCA图

library(scatterplot3d)
#继续使用上面的pcadata数据集
p1<-scatterplot3d(pcadata$PC1,pcadata$PC2,pcadata$PC3,
    color=rep(c("red","blue"),each=4),
    xlab=xlab,ylab=ylab,zlab=zlab,
    pch=rep(c(16,17),each=4),
    col.axis="grey50",
    mar=c(3.1,5.1,3.1,5.1),
    col.grid="grey",box=F)

#参数详解
#scatterplot3d(x, y=NULL, z=NULL, color=par("col"), pch=par("pch"),main=NULL, sub=NULL, xlim=NULL, ylim=NULL, zlim=NULL,xlab=NULL, ylab=NULL, zlab=NULL, scale.y=1, angle=40,axis=TRUE, tick.marks=TRUE, label.tick.marks=TRUE,x.ticklabs=NULL, y.ticklabs=NULL, z.ticklabs=NULL,y.margin.add=0, grid=TRUE, box=TRUE, lab=par("lab"),lab.z=mean(lab[1:2]), type="p", highlight.3d=FALSE,mar=c(5,3,4,3)+0.1, bg=par("bg"), col.axis=par("col.axis"),col.grid="grey", col.lab=par("col.lab"), cex.symbols=par("cex"), cex.axis=0.8 * par("cex.axis"),cex.lab=par("cex.lab"), font.axis=par("font.axis"),font.lab=par("font.lab"), lty.axis=par("lty"),lty.grid=par("lty"), lty.hide=NULL, lty.hplot=par("lty"),log="", asp=NA, ...)

#pch：代表点的形状
#main：图标题
#xlab：x轴标题
#angle：x轴和y轴的夹角
#tick.marks：刻度线
#grid：网格线
#box：外部方框 
#lab：x轴和y轴刻度线的个数，并不是严格控制
#type：p代表点，h代表垂直线，l代表线
#cex.symbol：代表形状的缩放
#font.axis：轴上文字的字形，分别1/2/3代替
#mar：代表图形的边距 

image

legend("bottom",legend=levels(pcadata$group),col=c("red","blue"),pch=c(16,17),inset=.75,horiz=T)
text(p1$xyz.convert(pcadata[2:4]),labels=pcadata[,1],pos=4,cex=0.6)

image