16是分析之主成分分析（PCA）实战部分
下面我们来做PCA分析，来做出很好看的图片
之前在插几句：非约束排序，不做统计检验，不能做！
# 安装包ggbiplot包
#install.packages("devtools",repo="http://cran.us.r-project.org")
#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一样画一个箭头，因为我没有找到简单的方法完成这件事，总会这样的，如果大家有好的办法，可以留言给我最好了，我将不胜感激。
最后还是那句话：学习永无止境，分享永不停歇！