|
下面我们来做PCA分析，来做出很好看的图片 之前在插几句：非约束排序，不做统计检验，不能做！ # 安装包ggbiplot包 #install.packages("devtools",repo="http://cran.us.") #library(devtools) #install_github("vqv/ggbiplot") library(ggbiplot) #修改工作目录 setwd("E:/Shared_Folder/HG_usearch_HG") #导入需要的作图文件 design =read.table("map_lxdjhg.txt", header=T, row.names= 1,sep="\t") head(design) # 读取OTU表并进行数据整理 otu_table =read.delim("otu_table8825.txt", row.names= 1, header=T, sep="\t") head(otu_table) # 过滤数据并排序 sub_design =design[rownames(design) %in% colnames(otu_table),] count = otu_table[,rownames(sub_design)] 以上我们导入数据并整理数据，接下来我们开始分析PCA： 函数prcomp默认选项center=T(当然中心化是一定要做的)，但是不做标准化，我们调整参数scale. = TRUE，将方差缩放为单位方差最为合理 # 进行PCA分析 otu.pca <-prcomp(t(count), scale. = TRUE) 这里我们使用的ggbiplot，obs.scale = 1, var.scale = 1,这两个参数表示我们使用标尺出图，ellipse= TRUE我们添加椭圆置信区间默认是0.95置信区间，可以修改吗？对不起，不行，写死了，已经在函数里面； p=ggbiplot(otu.pca,obs.scale = 1, var.scale = 1, groups = sub_design$SampleType,ellipse = TRUE,var.axes = T,circle = TRUE)+ scale_color_discrete(name = '') + theme(legend.direction = 'horizontal',legend.position = 'top') # p 下面这幅图算是将PCA的所有要素展示的非常好了，这里有一个新的东西，之前也没说过，就是： 平衡贡献圆（circle of equilibrium contribution），它的半径表示变量的向量长度对排序的平均贡献率，可以这么讲，如果变量的箭头线长度长于圆的半径，则它对排序空间的贡献大于变量的平均贡献率； 我们使用扩增子16s，未免otu数目过多，什么都看不清楚，自然需要我们将变量减少 # 因此我们首先关注一下高丰度菌的影响 # 转换原始数据为百分比 norm =t(t(count)/colSums(count,na=T)) * 100 # normalization to total 100 head(norm) # 筛选mad值大于0.5的OTU norm2 =norm[apply(norm,1,mean)>1,] dim(norm2) #重新计算PCA otu.pca <-prcomp(t(norm2)) p<-ggbiplot(otu.pca,obs.scale = 1, var.scale = 1, groups = sub_design$SampleType,ellipse = TRUE,var.axes = T,circle = TRUE) p #挑选0.001丰度以上的菌 norm3 =norm[apply(norm,1,mean)>0.1,] otu.pca1 <-prcomp(t(norm3)) p<-ggbiplot(otu.pca1,obs.scale = 1, var.scale = 1, groups = sub_design$SampleType,ellipse = TRUE,var.axes = T,circle = TRUE) p #当然我们仅仅展示样品点坐标也是可以的只需修改var.axes = F p<-ggbiplot(otu.pca1,obs.scale = 1, var.scale = 1, groups = sub_design$SampleType, ellipse= TRUE,var.axes = F) #可以修改主题，我修改了一份自己觉得美观的主题（我没有什么审美），大家有需要可以借鉴 p+theme_bw()+scale_colour_manual(values= mi)+labs(x=paste("PC 1 (", 61.7, "%)", sep=""), y=paste("PC 2 (", 14.5,"%)", sep=""))+ labs(title = "PCA")+ theme(plot.title = element_text(hjust =0.5))+ guides(color=guide_legend(title =NULL),shape=guide_legend(title = NULL))+ geom_hline(yintercept=0) +geom_vline(xintercept=0) ggbiplot是一款PCA分析结果可视化的R包工具，但是毕竟也有很多局限在里面的，比如： 我需要修改样品表示方式，我需要化多个置信区间，可能就不可以实现，所以将PCA分析结果提取出来，使用ggplot作图系统来作图会更加自由一些； 下面我们开始提取主成分分析结果 #现在我们来提取作图坐标 #首先是样品排序轴坐标 predict(otu.pca) #也可以使用下面这条命令提取 yangpin<-otu.pca$x yangpin=as.data.frame(yangpin) yangpin$SampleType=sub_design$SampleType #提取荷载坐标 bianliang<-otu.pca$rotation bianliang=as.data.frame(bianliang) #提取特征根,这里提供的并不是特征值而是标准差，需要求其平方才是特征值 eig=otu.pca$sdev eig=eig*eig #在这里我设定了随机种子，方便两种形式图形比较 set.seed(10) ggplot(data=yangpin,aes(x=yangpin$PC1,y=yangpin$PC2,group=SampleType,color=SampleType))+geom_point()+ stat_ellipse(type = "t", linetype =2)+ labs(x=paste("PC 1 (", format(100 *eig[1] / sum(eig), digits=4), "%)", sep=""), y=paste("PC 2 (", format(100 *eig[2] / sum(eig), digits=4), "%)", sep=""))+ labs(title = "PCA") #########我们将载荷也画上去 p+stat_ellipse(data=yangpin,aes(group=SampleType,color=SampleType),level = 0.95,type ="t", linetype = 2)+ geom_point(data=bianliang,aes(x=PC1,y=PC2),size=2)+ geom_text(data=bianliang,aes(x=PC1,y=PC2,label=spece),color="blue",hjust=1.2,vjust=-1.3,size=3) 这里我们没有像ggbioplot一样画一个箭头，因为我没有找到简单的方法完成这件事，总会这样的，如果大家有好的办法，可以留言给我最好了，我将不胜感激。 最后还是那句话：学习永无止境，分享永不停歇！ |
|
|
