Lefse分析R语言全套稳定版

Lefse分析R语言全套稳定版写在前面这段时间比较忙，宝贝女儿出生了，所以推文有些慢，但是出来就给大家好东西。比如下面的：lefse分析之前我写过很多了，用R语言重现分析内容，做树图，做柱状图等，但是之前的版本好多人用过后都说不够稳定，因此，这里我重新写一个lefse全套功能，可以解决之前的问题，具体如下：之前LDA判别分析参数设置比较少，这次完善之前出树图往往错误比较多，尤其是做真菌，细菌群落种同时具有古菌等，这是物种分类如果有多个界的情况下没有进行很好的处理，之前依赖的构建树函数换为Y叔的R包MicrobiotaProcess，用其中的函数构建树骨架。对骨架添加阴影，注释等，我改写了函数clade.anno为clade.anno_wt，这些都让树图出来的非常流畅。本次再不同分类等级之间使用Rank分隔符进行分隔，很好的匹配注释数据，后续可以添加更多花样本次更新，可以对不同分类等级的微生物数据都进行lefse分析（这个函数我单独写了下，参见三人成师培训课程，尚未公开分享，如感兴趣，可私聊我）实战library(ggpubr)library(patchwork)library(MicrobiotaProcess)library(ggClusterNet)library(phyloseq)library(ggtree)library(ggtreeExtra)library(tidyverse)data(ps)mytheme1=theme_bw()注意运行下面代码之前，将推文末尾的四个函数运行一下。第一步生成树骨架;p1<-p_base(ps,Top=100)p1第二步：进行LDA判别tablda=LDA_Micro(ps=ps,#phyloseq对象，gCLusterNet提供Top=100,#挑选丰度最高的前100个微生物进行分析p.lvl=0.05,#选择显著性阈值lda.lvl=2.5,#选择LDA判别额阈值seed=11,#设定随机种子，保证分析可重复adjust.p=T#是否矫正P自值)p<-lefse_bar(taxtree=tablda[[2]])p查看数据res=tablda[[2]]head(res)Pvaluesst__ASV_60.00051110f__Bacteria_Rank_Proteobacteria_Rank_Betaproteobacteria_Rank_Burkholderiales_Rank_Unassigned0.00076069g__Bacteria_Rank_Proteobacteria_Rank_Betaproteobacteria_Rank_Burkholderiales_Rank_Unassigned_Rank_Unassigned0.00076069s__Bacteria_Rank_Proteobacteria_Rank_Betaproteobacteria_Rank_Burkholderiales_Rank_Unassigned_Rank_Unassigned_Rank_Unassigned0.00076069c__Bacteria_Rank_Proteobacteria_Rank_Betaproteobacteria0.00391140o__Bacteria_Rank_Proteobacteria_Rank_Betaproteobacteria_Rank_Burkholderiales0.00391140第三步：注释树#注释树p2<-clade.anno_wt(p1,tablda[[1]],alpha=0.3,#设定阴影透明度anno.depth=7#设定级别高于等于7级才将注释添加到图上)p2函数#p1<-p_base(ps,Top=100)###tablda=LDA_Micro(ps=ps,Top=100,p.lvl=0.05,lda.lvl=2,seed=11,adjust.p=T)###注释树##p2<-clade.anno_wt(p1,tablda[[1]],alpha=0.3,anno.depth=7)##ggsave("./cs6.pdf",p2,width=15,height=10)p_base=function(ps,Top=300){alltax=ps%>%filter_OTU_ps(Top)%>%vegan_tax()%>%as.data.frame()alltax$OTU=s(alltax)head(alltax)alltax$Kingdom=paste(alltax$Kingdom,sep="_Rank_")alltax$Phylum=paste(alltax$Kingdom,alltax$Phylum,sep="_Rank_")alltax$Class=paste(alltax$Phylum,alltax$Class,sep="_Rank_")alltax$Order=paste(alltax$Class,alltax$Order,sep="_Rank_")alltax$Family=paste(alltax$Order,alltax$Family,sep="_Rank_")alltax$Genus=paste(alltax$Family,alltax$Genus,sep="_Rank_")alltax$Species=paste(alltax$Genus,alltax$Species,sep="_Rank_")alltax[is.na(alltax)]="Unknown"trda<-convert_to_treedata(alltax)p<-ggtree(trda,layout="circular",size=0.2,xlim=c(30,NA))+geom_point(pch=21,size=3,alpha=1,fill="#FFFFB3")p$data$lab2<-p$data$label%>%strsplit("_Rank_")%>%sapply(function(x)x[length(x)])p$data$lab2=gsub("st__","",p$data$lab2)p$data$nodeSize=1return(p)}LDA_Micro=function(ps=ps,Top=100,p.lvl=0.05,lda.lvl=2,seed=11,adjust.p=F){alltax=ps%>%filter_OTU_ps(Top)%>%vegan_tax()%>%as.data.frame()alltax$OTU=s(alltax)head(alltax)alltax$Kingdom=paste(alltax$Kingdom,sep="_Rank_")alltax$Phylum=paste(alltax$Kingdom,alltax$Phylum,sep="_Rank_")alltax$Class=paste(alltax$Phylum,alltax$Class,sep="_Rank_")alltax$Order=paste(alltax$Class,alltax$Order,sep="_Rank_")alltax$Family=paste(alltax$Order,alltax$Family,sep="_Rank_")alltax$Genus=paste(alltax$Family,alltax$Genus,sep="_Rank_")alltax$Species=paste(alltax$Genus,alltax$Species,sep="_Rank_")otu=ps%>%filter_OTU_ps(Top)%>%vegan_otu()%>%t()%>%as.data.frame()otu_tax=merge(otu,alltax,by="s",all=F)head(otu_tax)rank1<-otu_tax%>%group_by(Kingdom)%>%summarise_if(is.numeric,sum,na.rm=TRUE)colnames(rank1)[1]="id"rank1$id=paste("k__",rank1$id,sep="")rank2<-otu_tax%>%group_by(Phylum)%>%summarise_if(is.numeric,sum,na.rm=TRUE)colnames(rank2)[1]="id"rank2$id=paste("p__",rank2$id,sep="")rank3<-otu_tax%>%group_by(Class)%>%summarise_if(is.numeric,sum,na.rm=TRUE)colnames(rank3)[1]="id"rank3$id=paste("c__",rank3$id,sep="")rank4<-otu_tax%>%group_by(Order)%>%summarise_if(is.numeric,sum,na.rm=TRUE)colnames(rank4)[1]="id"rank4$id=paste("o__",rank4$id,sep="")rank5<-otu_tax%>%group_by(Family)%>%summarise_if(is.numeric,sum,na.rm=TRUE)colnames(rank5)[1]="id"rank5$id=paste("f__",rank5$id,sep="")rank6<-otu_tax%>%group_by(Genus)%>%summarise_if(is.numeric,sum,na.rm=TRUE)colnames(rank6)[1]="id"rank6$id=paste("g__",rank6$id,sep="")rank7<-otu_tax%>%group_by(Species)%>%summarise_if(is.numeric,sum,na.rm=TRUE)colnames(rank7)[1]="id"rank7$id=paste("s__",rank7$id,sep="")rank8<-otu_tax%>%group_by(OTU)%>%summarise_if(is.numeric,sum,na.rm=TRUE)colnames(rank8)[1]="id"rank8$id=paste("st__",rank8$id,sep="")#合并8个分类级all=rbind(rank1,rank2,rank3,rank4,rank5,rank6,rank7,rank8)head(all)#--LDA排序#--------data1=as.data.frame(all)s(data1)=data1$iddata1$id=NULL#-构建phylose对象ps_G_graphlan=phyloseq(otu_table(as.matrix(data1),taxa_are_rows=TRUE),sample_data(ps))ps_G_graphlan#----提取OTU表格otu=as.data.frame((vegan_otu(ps_G_graphlan)))otu[otu==0]<-1map=as.data.frame(sample_data(ps_G_graphlan))#otu=(otu_table)claslbl=map$Group%>%as.factor()set.seed(seed)#KWranksumtestrawpvalues<-apply(otu,2,function(x)kruskal.test(x,claslbl)$p.value);#--得到计算后得到的p值ord.inx<-order(rawpvalues)rawpvalues<-rawpvalues[ord.inx]clapvalues<-p.adjust(rawpvalues,method="fdr")#p.adjustwil_datadf<-as.data.frame(otu[,ord.inx])head(wil_datadf)ldares<-MASS::lda(claslbl~.,data=wil_datadf)#ldaresldamean<-as.data.frame(t(ldares$means))ldameanclass_no<<-length(unique(claslbl))ldamean$max<-apply(ldamean[,1:class_no],1,max);ldamean$min<-apply(ldamean[,1:class_no],1,min);#---计算LDAldamean$LDAscore<-signif(log10(1+abs(ldamean$max-ldamean$min)/2),digits=3);head(ldamean)a=rep("A",length(ldamean$max))for(iin1:length(ldamean$max)){name=colnames(ldamean[,1:class_no])a[i]=name[ldamean[,1:class_no][i,]%in%ldamean$max[i]]}ldamean$class=atem1=s(ldamean)tem1%>%as.character()ldamean$Pvalues<-signif(rawpvalues[match(s(ldamean),names(rawpvalues))],digits=5)ldamean$FDR<-signif(clapvalues,digits=5)resTable<-ldameanrawNms<-rownames(resTable);rownames(resTable)<-gsub("`",'',rawNms);if(adjust.p){de.Num<-sum(clapvalues<=p.lvl&ldamean$LDAscore>=lda.lvl)}else{de.Num<-sum(rawpvalues<=p.lvl&ldamean$LDAscore>=lda.lvl)}if(de.Num==0){current.msg<<-"Nosignificantfeatureswereidentifiedwithgivencriteria.";}else{current.msg<<-paste("Atotalof",de.Num,"significantfeatureswithgivencriteria.")}print(current.msg)#sortbypvalueord.inx<-order(resTable$Pvalues,resTable$LDAscore)resTable<-resTable[ord.inx,,drop=FALSE]resTable<-resTable[,c(ncol(resTable),1:(ncol(resTable)-1))]resTable<-resTable[,c(ncol(resTable),1:(ncol(resTable)-1))]#resTable%>%tail()if(adjust.p){taxtree=resTable[clapvalues<=p.lvl&ldamean$LDAscore>=lda.lvl,]}else{#taxtree=resTable[ldamean$Pvalues<=p.lvl&ldamean$LDAscore>=lda.lvl,]taxtree=resTable[ldamean$Pvalues<=p.lvl,]}#-提取所需要的颜色colour=c('darkgreen','red',"blue","#4DAF4A","#984EA3","#FF7F00","#FFFF33","#A65628","#F781BF")selececol=colour[1:length(levels(as.factor(taxtree$class)))]names(selececol)=levels(as.factor(taxtree$class))A=rep("a",length(s(taxtree)))for(iin1:length(s(taxtree))){A[i]=selececol[taxtree$class[i]]}taxtree$color=A#taxtree<-taxtree[s(taxtree)!="k__Bacteria",]#node_ids<-p0$data#anno<-rep("white",nrow(p1$data))lefse_lists=data.frame(node=s(taxtree),color=A,Group=taxtree$class,stringsAsFactors=FALSE)return(list(lefse_lists,taxtree))}##注释树#p1<-clade.anno_wt(p,lefse_lists,alpha=0.3,anno.depth=7)##ggsave("./cs7.pdf",p1,width=15,height=10)##p1<-clade.anno_wt(p,lefse_lists,alpha=0.3,anno.depth=8)#p1##ggsave("./cs9.pdf",p1,width=10,height=10)#gtree=p1#anno.data=tablda#alpha=0.3#anno.depth=7#anno.x=10#anno.y=40clade.anno_wt<-function(gtree,anno.data,alpha=0.2,anno.depth=5,anno.x=10,anno.y=40){short.labs<-c(letters,paste(letters,1:500,sep=""))get_offset<-function(x){(x*0.2+0.2)^2}get_angle<-function(node){data<-gtree$datasp<-tidytree::offspring(data,node)$nodesp2<-c(sp,node)sp.df<-data[match(sp2,data$node),]mean(range(sp.df$angle))}anno.data<-arrange(anno.data,node)hilight.color<-anno.data$colornode_list<-anno.data$nodenode_ids<-(gtree$data%>%filter(label%in%node_list)%>%arrange(label))$nodeanno<-rep("yellow",nrow(gtree$data))#---添加阴影#-------i=1for(iin1:length(node_ids)){n<-node_ids[i]color<-hilight.color[i]anno[n]<-colormapping<-gtree$data%>%filter(node==n)nodeClass<-as.numeric(mapping$nodeDepth)offset<-get_offset(nodeClass)gtree<-gtree+geom_hilight(node=n,fill=color,alpha=alpha,extend=offset)}#gtree$layers<-rev(gtree$layers)#gtree<-gtree+geom_point2(aes(size=I(nodeSize)),fill=anno,#shape=21)short.labs.anno<-NULLi=1#gtree$data#--添加标签#--------for(iin1:length(node_ids)){n<-node_ids[i]mapping<-gtree$data%>%filter(node==n)nodeClass<-as.numeric(mapping$nodeDepth)if(nodeClass<=anno.depth){lab<-short.labs[1]short.labs<-short.labs[-1]if(is.null(short.labs.anno)){short.labs.anno=data.frame(lab=lab,annot=mapping$lab2,stringsAsFactors=F)}else{short.labs.anno=rbind(short.labs.anno,c(lab,mapping$lab2))}}else{lab<-mapping$lab2}offset<-get_offset(nodeClass)-0.4angle<-get_angle(n)+90gtree<-gtree+geom_cladelabel(node=n,label=lab,angle=angle,fontsize=1+sqrt(nodeClass),offset=offset,barsize=NA,hjust=0.5)}if(!is.null(short.labs.anno)){anno_shapes=sapply(short.labs.anno$lab,utf8ToInt)stable.p<-ggpubr::ggtexttable(short.labs.anno,rows=NULL,theme=ttheme(colna