查看原文
其他

16是分析之主成分分析(PCA)实战部分

文涛 微生信生物 2022-05-08

下面我们来做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)

这里我们使用的ggbiplotobs.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.5OTU

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)

#挑选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一样画一个箭头,因为我没有找到简单的方法完成这件事,总会这样的,如果大家有好的办法,可以留言给我最好了,我将不胜感激。


最后还是那句话:学习永无止境,分享永不停歇!

您可能也对以下帖子感兴趣

文章有问题?点此查看未经处理的缓存