suppressWarnings(dir.create("~/tmp"))
.libPaths("~/tmp")
library(mmmgee)
library(statVisual)
library(msir)
source("http://viva1109.duckdns.org/RFunctions/mTMAT_funcs.R",encoding = "UTF-8")

mTMAT<-function(Data_Read,Meta_Table,Phenotype,ID,Time,Cov,Tree,Taxonomy_Table,Taxonomy_level="Genus",ncore=1,str="independence",CompositionalBias_correction="Default",mTMAT_type="IM",taxon_filter_zero_prop=0.25,Total_Readcount=NULL){
  #browser()
  names(Taxonomy_Table)[-1]<-c("Kingdom", "Phylum", "Class", "Order", "Family", "Genus" ,"Species")
  if(is.null(Total_Readcount)){
	total_readcount<-apply(Data_Read,2,sum)
  }else{
	total_readcount<-Total_Readcount	
  }
  
  ind_sample<-match(Meta_Table[,1],colnames(Data_Read))
  Meta_Table$total_readcount<-total_readcount[ind_sample]
  otu_df<-Data_Read[,ind_sample]
  otu_mat = as.matrix(otu_df)

  zero.pr = apply(otu_mat, 1, function(x) {length(x[x != 0])/length(x)} )
  df_zero.pr = data.frame(zero.p=sort(zero.pr, decreasing = T))
  ind_otu_remain<-df_zero.pr$zero.p>taxon_filter_zero_prop
  OTU_id_list<-rownames(df_zero.pr)[ind_otu_remain]

  otu_mat_sel = t(otu_mat[OTU_id_list, ])
  otuALLALLcnt<-apply(otu_df,2,sum)
  otu_countall<-apply(otu_mat_sel,1,sum)
# Other_generaCnt<-otuALLALLcnt-otu_countall


  y_vec<-as.numeric(Meta_Table[,Phenotype])
  id_vec<-as.character(Meta_Table[,ID])
  trc_vec<-Meta_Table$total_readcount
  time_var<-Time
  covResult<-makeCovMat(Cov,Meta_Table)

  if(mTMAT_type=="IM"){
	type=16
	ind_type<-1
  }else{
	type=15
	ind_type<-2
  }
  if(!is.rooted(Tree)){
    tree_ori<-root(Tree, 1, r = TRUE)
  }else{
    tree_ori<-Tree
  }
  tree_chosen<-drop.tip(tree_ori,setdiff(tree_ori$tip.label,OTU_id_list))

  Taxo_togo<-Taxonomy_Table[match(OTU_id_list,Taxonomy_Table[,1]),]
  cnt_taxon<-table(Taxo_togo[,Taxonomy_level])
  taxon_anal_list<-names(sort(cnt_taxon,decreasing = T))
	
  outTMATp<-suppressWarnings({mclapply(1:length(taxon_anal_list),function(j,cbd_clinical){
    target_taxon<-taxon_anal_list[j]
    col_ind<-which(Taxo_togo[,Taxonomy_level]==target_taxon)
  
    if(length(col_ind)==1){
      INDI_SINGLE<-T
    }else{
      INDI_SINGLE<-F
    }
    X_P<-matrix(otu_mat_sel[, col_ind],ncol=length(col_ind))
    colnames(X_P)<-colnames(otu_mat_sel)[col_ind]
    
    # tree_chosen
    if(is.rooted(tree_chosen)){
      tree_chosen_1<-tree_chosen
    }else{
      tree_chosen_1<-multi2di(tree_chosen)
    }
   
    tree_chosen2_bf<-drop.tip(tree_chosen_1,setdiff(tree_chosen_1$tip.label,colnames(X_P)))
    
    
    if(!is.rooted(tree_chosen2_bf)){
      tree_chosen2<-root(tree_chosen2_bf, 1, r = TRUE)
    }else{
      tree_chosen2<-tree_chosen2_bf
    }
    
    if(INDI_SINGLE){
      genera_vec<-X_P
      dim(genera_vec)<-NULL
    }else{
      genera_vec<-apply(X_P,1,sum)
    }
    
    Counts_C<-matrix(trc_vec-genera_vec,ncol=1)
    
    indi_tAnal<-list()
    indi_tAnal$subtree_tmp<-list()
    indi_tAnal$subtree_tmp[[1]]<-tree_chosen2
  
  
    
    # dataSet<-simulSet_by_methods_fixed(y=y_vec,dataSet=data.frame(X_P,check.names = F),ttr=trc_vec,type=16,method="TMAT") 
    dataSet<-simulSet_by_methods_fixed(y=y_vec,dataSet=data.frame(X_P,check.names = F),ttr=trc_vec,type=16,method="TMAT") 
    
    # TMAT_pval_bygenus_bf<-TMAT_func(dataSet,indi_tAnal,type=16,total.reads=total.reads,conti=F,cov=NULL,getBeta=T,DetailResult=T)
    
    # output$pval
    if(dim(dataSet[[2]][[1]])[2]==1){
      DMsC<-TMAT_func_dm2(dataSet,indi_tAnal,type=16,total.reads=dataSet$ttr,conti=F,cov=NULL,getBeta=T,DetailResult=F)
      DMsC_2<-DMsC
      DMsC_2$dms<-list(DMsC_2$dms)
      
      DMsC15<-TMAT_func_dm2(dataSet,indi_tAnal,type=15,total.reads=dataSet$ttr,conti=F,cov=NULL,getBeta=T,DetailResult=F)
      DMsC15_2<-DMsC15
      DMsC15_2$dms<-list(DMsC15_2$dms)
      
    }else{
      DMsC<-TMAT_func_dm2(dataSet,indi_tAnal,type=16,total.reads=dataSet$ttr,conti=F,cov=NULL,getBeta=T,DetailResult=F)
      DMsC_2<-DMsC
      
      DMsC15<-TMAT_func_dm2(dataSet,indi_tAnal,type=15,total.reads=dataSet$ttr,conti=F,cov=NULL,getBeta=T,DetailResult=F)
      DMsC15_2<-DMsC15
    }

  
    p_result<-t(sapply(1:length(DMsC_2$dms),function(ind_node){
      
      
      DATA<-data.frame(response=DMsC_2$dms[[ind_node]],phenotype=y_vec,subject=id_vec,covResult)
      
      DATA_sort <- DATA[order(DATA$subject),]
      formula <- formula(paste(c("response~phenotype",names(covResult)),collapse="+"))
      D0<-matrix(c(0,1,rep(0,dim(covResult)[2])),nrow=1) 
      m2<-geem2(formula,id=subject,data=DATA_sort,family=gaussian,corstr=str)
      res1<-mmmgee.test(x=m2,L=list(D0),r=list(c(0)),statistic="score",type="quadratic",biascorr=T)
      c(res1$test$global$p,m2$beta[2])
    }))
    mTMAT_res<-exact_pval(p_result[,1])
  
    p_result15<-t(sapply(1:length(DMsC15_2$dms),function(ind_node){
      
      DATA<-data.frame(response=DMsC_2$dms[[ind_node]],phenotype=y_vec,subject=id_vec,covResult)
  
      DATA_sort <- DATA[order(DATA$subject),]
      formula <- formula(paste(c("response~phenotype",names(covResult)),collapse="+"))
      D0<-matrix(c(0,1,rep(0,dim(covResult)[2])),nrow=1)
      m2<-geem2(formula,id=subject,data=DATA_sort,family=gaussian,corstr=str)
      res1<-mmmgee.test(x=m2,L=list(D0),r=list(c(0)),statistic="score",type="quadratic",biascorr=T)
      c(res1$test$global$p,m2$beta[2])
    }))
    mTMAT_res15<-exact_pval(p_result15[,1])
    beta_TMAT16<-p_result[which.min(p_result[,1]),2]
    beta_TMAT15<-p_result15[which.min(p_result15[,1]),2]
    RES_TMAT16<-list(mTMAT_res,pvals_betas=cbind(p_result[,1],p_result[,2]),output=DMsC,dataSet=dataSet)
    RES_TMAT15<-list(mTMAT_res15,pvals_betas=cbind(p_result15[,1],p_result15[,2]),output=DMsC15,dataSet=dataSet)
    
    RES_pvalues<-c(mTMAT_res,mTMAT_res15)
    RES_beta<-c(beta_TMAT16,beta_TMAT15)

    list(RES_TMAT16,RES_TMAT15,RES_pvalues,RES_beta)
  },cbd_clinical=Meta_Table,mc.cores=ncore)})
  if(class(outTMATp[[1]])=="try-error"){
	#print(outTMATp[[1]])
	cat("[mTMAT failed to converge]\n")
	if(str_detect(outTMATp[[1]][1],"singular")){
		cat("Categorical covariates made it hard to converge. One of the categorical covariates needs to be excluded. \n")
	}
	cat("The error message is returned.\n")
	return(outTMATp[[1]])
  }else{
	  reTable_pvalus<-t(sapply(outTMATp,function(data){
	    data[[3]]
	  }))
	  reTable_betas<-t(sapply(outTMATp,function(data){
	    data[[4]]
	  }))
	  upperlevel<-names(Taxo_togo)[which(names(Taxo_togo)==Taxonomy_level)-1]
	  FDR_result<-apply(reTable_pvalus,2,p.adjust,method="fdr")
	  output_RealTable<-cbind(Taxo_togo[match(taxon_anal_list,Taxo_togo[,Taxonomy_level]),upperlevel],data.frame(taxon_anal_list),reTable_betas[,ind_type],reTable_pvalus[,ind_type],FDR_result[,ind_type])
	  nametogo<-c("mTMAT","mTMAT")[ind_type]
	  names(output_RealTable)<-c(upperlevel,Taxonomy_level,paste0(nametogo,"_beta"),paste0(nametogo,"_pvalue"),paste0(nametogo,"_FDR"))
	  
	  output<-list(resultTable=output_RealTable, output_raw=list(outTMATp=outTMATp,taxonomy_info=Taxo_togo,tree=tree_chosen,time_vec=Meta_Table[time_var],Taxonomy_level=Taxonomy_level,mTMAT_type=mTMAT_type))
	  return(output)
  }
  

  
}


plot_mTMAT<-function(target_taxon=NULL,taxon_index=NULL,output_mTMAT,phenotype_level=c("0","1"),x_label="Time",Cex_node=7,Cex_nodetxt=1.09, ...){
	#### bookmark Input 1 ####
	#browser() 
	
	if(is.null(target_taxon) & is.null(taxon_index)){
		print("target_taxon or taxon_index need to be specified")
	}else{
		if(is.null(target_taxon)){
			ind_direct_taxon<-taxon_index
			#### bookmark Input 2 ####
			target_taxon<-as.character(output_mTMAT$resultTable[ind_direct_taxon,2])
		}else{
			ind_direct_taxon<-which(output_mTMAT$resultTable[,2]==target_taxon)
		}
	}
	

	target_taxon_out<-gsub(":","_",target_taxon)
	# try(dev.off())
	# png(paste0("mTMAT_practice_",target_taxon_out,"_new.png") ,width=1200*scaler,height=700*scaler, type = "cairo")

	#### bookmark Input 3 outTMATpoutput, outTMATdataset####
	tmp_output_raw<-output_mTMAT$output_raw

	if(tmp_output_raw$mTMAT_type=="IM"){
		type<-16
		ind_type<-1
	}else{
		type<-15
		ind_type<-2
	}


	tmp_output<-tmp_output_raw$outTMATp[[ind_direct_taxon]]
	testMainResult<-tmp_output[[ind_type]][[3]]
	datasetTable<-tmp_output[[ind_type]][[4]]
	output_beta<-tmp_output[[ind_type]][[2]]
	Taxonomy_level<-tmp_output_raw$Taxonomy_level
	Taxo_togo<-tmp_output_raw$taxonomy_info
	col_ind<-which(Taxo_togo[,Taxonomy_level]==target_taxon)
	chosen_otuID<-Taxo_togo$X[col_ind]
	tree_chosen<-tmp_output_raw$tree
	time_vec<-tmp_output_raw$time_vec
	names(time_vec)<-"time"

	margins<-c(0.2,1.5,1.5,1.5)
	par(mfrow=c(1,1), mar=margins, oma=c(0, 0, 0, 0))
	layout(matrix(c(1,2,3,4,5,6), 2, 3, byrow = TRUE),widths = c(2,10, 8),heights = c(3,2))

	contrast<-c(1,2)

	#### bookmark Input 4 tree_chosen ####
	if(is.binary(tree_chosen)){
		tree_chosen_1<-tree_chosen
	}else{
		tree_chosen_1<-multi2di(tree_chosen)
	}

	tree_chosen2<-drop.tip(tree_chosen_1,setdiff(tree_chosen_1$tip.label,chosen_otuID))
	set.seed(0)
	tree_chosen2$edge.length[tree_chosen2$edge.length<0.1]<-runif(1,0.07,0.12)

	tiplabels_ready<-paste0(Taxo_togo[match(tree_chosen2$tip.label,Taxo_togo$X),"Species"]," (",tree_chosen2$tip.label,")")
	ind_uncultered<-str_detect(tiplabels_ready,"uncultured")
	tiplabels_ready[ind_uncultered]<- paste0(tiplabels_ready," (",tree_chosen2$tip.label,")")[ind_uncultered]
	tiplabels_ready<-c(paste0("Other ",Taxonomy_level),tiplabels_ready)
	if(length(col_ind)!=1){
		tiplabels_ready<-c(tiplabels_ready,paste0("[",Taxonomy_level,"] ",target_taxon))
	}

	betah<-output_beta[,2]
	# vvaall<-exp(betah)
	vvaall<-exp(betah)
	# vvaallTEXT<-betah

	# colfunc <- colorRampPalette(c("blue","white","red"))
	colfunc <- colorRampPalette(c("#00bec6","white","#ff756b"))

	ind_col<-ceiling(vvaall*199.5-100)
	ind_col[ind_col>200]<-200
	ind_col[ind_col<0]<-1
	col_picked<-colfunc(200)
	col_nodes<-col_picked[ind_col]
	col_nodes2<-col_nodes[c(length(col_nodes),1:(length(col_nodes)-1))]

	x <- rtree(2, tip.label = LETTERS[1:2])
	x$edge.length[1]<-0.02
	x$edge.length[2]<-0
	tree.Groups_toplot<-bind.tree(x,tree_chosen2,where = 1)

	n_tips<-Ntip(tree.Groups_toplot)

	xl <- 1
	yb <- 1
	xr <- 1.5
	yt <- 2
	ncol_picked<-20
	ncol_texts<-10
	col_picked<-colfunc(ncol_picked)

	plot(NA,type="n",xlim=c(1,2),ylim=c(1,2),xaxt="n",yaxt="n",bty="n",xlab="",ylab="")
	mtext(expression(hat(beta[frac("■","●")])), side=3, adj=0, line=-1, cex=1.1, font=2);

	rect(
	xl,
	head(seq(yb,yt,(yt-yb)/ncol_picked),-1),
	xr,
	tail(seq(yb,yt,(yt-yb)/ncol_picked),-1),
	col=col_picked
	)


	ttt_bf<-seq(1,200,length=ncol_texts)

	ttt<-round(log(((ttt_bf-1)+100)/199.5)*10)/10
	mtext(ttt,side=2,at=tail(seq(yb,yt,(yt-yb)/ncol_texts),-1)-0.05,las=2,cex=0.7)
	p_phy<-plot_phylo(tree.Groups_toplot,main=target_taxon,direction="downwards",show.tip.label = FALSE,edge.width=2,font=6,cex.main=1.5, ...)



	p_vals<-output_beta[,1]

	if(!min(p_vals) %in% p_vals[-length(p_vals)]){
		ind_minP_node<- 0
		col_level<-rep(1,n_tips)
		if(as.logical(testMainResult$dir[which.min(p_vals)])){
			col_level[-n_tips]<-2
			col_level[n_tips]<-3
		}else{
			col_level[-n_tips]<-3
			col_level[n_tips]<-2
		}
		# tree.Groups_toplot$node.label
		makedm<-function(ind_node){
			C1<-apply(datasetTable$simData$X_P,1,sum)
			C2<-datasetTable$simData$X_P_comp
			dm<- type_to_dm(C1=C1,C2=C2,type=type,total.reads=datasetTable$ttr)[,1]
			
		}
		nodes_left<-1:length(p_vals)
		nodes_right<-0


	}else{
		ind_minP_node<- which(p_vals[-length(p_vals)]==min(p_vals))

		# title("Title text", adj = 0, line = 0)

		target_plot<-testMainResult$testnodes[[which.min(p_vals)]] 
		col_level<-rep(1,n_tips)
		if(as.logical(testMainResult$dir[which.min(p_vals)])){
			col_level[which(target_plot$group1)]<-2
			col_level[which(target_plot$group2)]<-3
		}else{
			col_level[which(target_plot$group1)]<-3
			col_level[which(target_plot$group2)]<-2
		}
		makedm<-function(ind_node){
			C1<-datasetTable$simData$X_P%*%testMainResult$testnodes[[ind_node]]$group1
			C2<-datasetTable$simData$X_P%*%testMainResult$testnodes[[ind_node]]$group2
			dm<- type_to_dm(C1=C1,C2=C2,type=type,total.reads=datasetTable$ttr)[,1]
			
		}

		nodes_left<-which(target_plot$group1)
		nodes_right<-which(target_plot$group2)

	}


	col_level2<-col_level[c(length(col_level),1:(length(col_level)-1))]
	co_alpha<-0.7
	co_black<-function(co_alpha){rgb(0, 0, 0, alpha=co_alpha)}
	co_red<-function(co_alpha){rgb(1, 0, 0, alpha=co_alpha)}
	co_blue<-function(co_alpha){rgb(0, 0, 1, alpha=co_alpha)}
	coco<-c("white","#ff756b","#00bec6")

	co_alpha<-1
	col_nodes_picked<-sapply(data.frame(col2rgb(col_nodes2)),function(data){
		rgb(data[1],data[2],data[3],alpha=co_alpha,maxColorValue=255)
	})
	mycol<-coco[col_level2]
	pchs<-rep(21,Nnode(tree.Groups_toplot))

	p_vals2_bf<-p.adjust(p_vals[c(length(p_vals),1:(length(p_vals)-1))],method="fdr")
	p_vals2<-round(p_vals2_bf*1000000)/1000000


	tt_targetX<-p_phy[[1]]
	targetX<-tt_targetX[(n_tips+1):length(tt_targetX)]
	tt_targetY<-p_phy[[2]]
	targetY<-tt_targetY[(n_tips+1):length(tt_targetY)]

	text(targetX,targetY-max(tt_targetY)/15,paste0("P-value: ",p_vals2),col="black" ,bg=col_nodes2,cex=1,srt=0)


	node_txt<-paste0("k=",1:Nnode(tree.Groups_toplot)-1)

	Cex_node_seq<-rep(Cex_node,length(node_txt))
	Cex_node_seq[which.min(p_vals2)]<-1.25*Cex_node
	Cex_nodetxt_seq<-rep(Cex_nodetxt,length(node_txt))
	Cex_nodetxt_seq[which.min(p_vals2)]<-1.25*Cex_nodetxt

	nodelabels(pch=pchs, col="black", bg=col_nodes2, cex=Cex_node_seq,srt=0, frame = "none")
	nodelabels(node_txt,col="black" ,bg=col_nodes2,cex=Cex_nodetxt_seq,srt=0, frame = "none")


	if(n_tips==2){
		tmp1<-mycol[1]
		tmp2<-mycol[2]
		mycol<-c(tmp2,tmp1)

		tips_txt<-paste0("m=",c(1,0))
		pch_seq<-rep(15,length(tips_txt))
		pch_seq[2]<-16
		tiplabels(pch=pch_seq, col="black", cex=Cex_node/3,srt=0,frame="none")
		tiplabels(tips_txt,adj=c(0.5,2),pch="", col="black", cex=Cex_nodetxt,srt=0,frame="none")  
		
	}else{
		tips_txt<-paste0("m=",c(1:n_tips-1))
		pch_seq<-rep(15,length(tips_txt))
		pch_seq[nodes_right+1]<-16
		tiplabels(pch=pch_seq, col="black", cex=Cex_node/3,srt=0,frame="none")
		tiplabels(tips_txt,adj=c(0.5,2),pch="", col="black", cex=Cex_nodetxt,srt=0,frame="none")  

	}

	nlabel<-vector("character",Nnode(tree.Groups_toplot))
	df_orinm<-tree_chosen2$tip.label
	dftmp_bf<-data.frame(datasetTable$simData$X_P[,tree_chosen2$tip.label])

	if(length(col_ind)!=1){
		sumVal<-apply(dftmp_bf,1,sum)
		dftmp<-cbind(dftmp_bf,sumVal)

		if(length(df_orinm)>5){
			str_m<-paste0(1,",...,",length(df_orinm))
		}else{
			str_m<-paste0(1:length(df_orinm),collapse=",")
		}

		names(dftmp)<-c(paste0("m=",1:length(df_orinm)),paste0("m=",str_m))
	}else{
		dftmp<-dftmp_bf
		str_m<-""

		names(dftmp)<-paste0("m=",1:length(df_orinm))
	}


	stacked<-stack(dftmp) 

	groups_for_plot_bf<-phenotype_level[datasetTable$y+1]
	groups_for_plot_bf<-factor(groups_for_plot_bf)

	groups_for_plot<-factor(groups_for_plot_bf,levels = phenotype_level[length(phenotype_level):1])
	cbd_stacked<-cbind(stacked,groups_for_plot)
	names(cbd_stacked)<-c("LogCPM","Leaf_nodes","groups")


	# p <- ggplot(data = cbd_stacked, aes(x = Leaf_nodes, y = LogCPM)) + 
	#   geom_boxplot(aes(fill = groups), width = 0.8) + theme_bw()+ theme(legend.position="bottom",plot.title = element_text(hjust = 0.5))+labs(title="Log counts per million by groups",y="Log counts per million",x="Leaf nodes")



	df_dm_tmp_tips<-cbind(cbd_stacked,time=time_vec)
	loess_bf_tip<-data.frame(LogCPM=df_dm_tmp_tips$LogCPM , y_points=df_dm_tmp_tips$LogCPM, Leaf_nodes=df_dm_tmp_tips$Leaf_nodes, groups=df_dm_tmp_tips$groups,time=df_dm_tmp_tips$time)

	spdf_tip<-split(loess_bf_tip,paste0(as.character(loess_bf_tip$Leaf_nodes),"_",loess_bf_tip$groups))
	sp_les_tip<-lapply(spdf_tip,function(data){
	loess_output<-loess.sd(x=data$time, y =data$LogCPM, nsigma = 1,span=0.75)
	cbind(data.frame(x=loess_output$x,y=loess_output$y,upper=loess_output$upper,lower=loess_output$lower),data)
	})

	outToplot_tips<-do.call("rbind",sp_les_tip)


	if(length(unique(outToplot_tips$x))>3){

		# png("lala0406_tips.png",height=1024,width=1024)

		p1.0 <- ggplot(outToplot_tips, aes(x=x, y=y,color=groups,group=groups,fill=groups)) +geom_line()
		p1.2 <- p1.0 + geom_ribbon(aes(ymax=upper, ymin=lower,color=groups,fill=groups), alpha=1/5)
		p1.3 <- p1.2 +  geom_point(aes(y=y_points))
		p1.4 <- p1.3 + theme_classic()+ theme(legend.position="bottom",plot.title = element_text(hjust = 0.5))+labs(title="logCPM for each leaf node",y="logCPM",x=x_label)
		p<-p1.4+facet_wrap(~Leaf_nodes, scales="free_y", nrow=1)

		# p
		# dev.off()
		# system("scp lala0406_tips.png ng:~")

	}else{

	}


	vp <- viewport(height = unit(0.4,"npc"), width=unit(0.6, "npc"), 
		 just = c("left","top"),   y = 0.4, x = 0)
	print(p, vp = vp)

	

	fontsize=min((21-n_tips)/2,5)
	if(length(col_ind)!=1){
		text<-paste(paste0("m=",c((1:n_tips-1),str_m),": ",tiplabels_ready),collapse="\n")  
	}else{
		text<-paste(paste0("m=",c((1:n_tips-1)),": ",tiplabels_ready),collapse="\n")
	}

	sp <- ggplot(NULL, aes(0, 1, label = text))+geom_point(pch="")
	      p2<-sp + geom_text(hjust=0,size=fontsize) + theme_bw()+ xlim(c(0, 1))+theme(axis.line=element_blank(),
	  axis.text.x=element_blank(),
	  axis.text.y=element_blank(),
	  axis.ticks=element_blank(),
	  axis.title.x=element_blank(),
	  axis.title.y=element_blank(),
	  legend.position="none",
	  panel.background=element_blank(), 
	  panel.border=element_blank(), 
	  panel.grid.major=element_blank(),
	  panel.grid.minor=element_blank(),
	  plot.background=element_blank())
	      


	vp <- viewport(height = unit(0.4,"npc"), width=unit(0.6, "npc"), 
		 just = c("left","top"),   y = 0.4, x = 0.6)
	print(p2, vp = vp)

	AllNodesPlot<-F
	if(AllNodesPlot){
		temp_dm<-lapply(2:(length(p_vals)),makedm)
	}else{
		temp_dm<-lapply(ind_minP_node,makedm)
	}


	df_dm<-data.frame(do.call("cbind",temp_dm))



	if(dim(df_dm)[2]!=1){
		stacked_dm<-stack(df_dm) 
		names(stacked_dm)<-c("LogCPM","Nodes")
		names(stacked_dm)[1]<-"LogCPM"
		df_dm<-cbind(stacked_dm,data.frame(groups=groups_for_plot))
		p3 <- ggplot(data = df_dm, aes(x = Nodes, y = LogCPM)) + 
		geom_boxplot(aes(fill = groups), width = 0.8) + theme_bw()+ theme(legend.position="bottom",plot.title = element_text(hjust = 0.5))+labs(title=paste("Log ratio of logCPMs of the each nodes"),y="Log ratio of log CPM")
	}else{
		stacked_dm<-df_dm
		names(stacked_dm)[1]<-"LogCPM"
		df_dm<-cbind(stacked_dm,data.frame(groups=groups_for_plot))

		df_dm_tmp<-cbind(stacked_dm,data.frame(groups=groups_for_plot,time=time_vec))
		loess_bf<-data.frame(LogCPM=df_dm_tmp$LogCPM , y_points=df_dm_tmp$LogCPM,groups=df_dm_tmp$groups,time=df_dm_tmp$time)

		spdf<-split(loess_bf,loess_bf$groups)
		sp_les<-lapply(spdf,function(data){
		loess_output<-loess.sd(x=data$time, y =data$LogCPM, nsigma = 1,span=0.75)
		cbind(data.frame(x=loess_output$x,y=loess_output$y,upper=loess_output$upper,lower=loess_output$lower),data)
		})

		outToplot<-do.call("rbind",sp_les)
		# loess_output<-loess.sd(x=time, y =df_dm_tmp$LogCPM, nsigma = 1)



		if(length(unique(outToplot$x))>3){
		p1.0 <- ggplot(outToplot, aes(x=x, y=y,color=groups,group=groups,fill=groups)) +geom_line()
		p1.2 <- p1.0 + geom_ribbon(aes(ymax=upper, ymin=lower,color=groups,fill=groups), alpha=1/5)
		p1.3 <- p1.2 +  geom_point(aes(y=y_points))
		p3 <- p1.3 + theme_classic()+ theme(legend.position="bottom",plot.title = element_text(hjust = 0.5))+labs(title=substitute(log~(frac(d,l))~d2, list(d2=paste0(" \nfor the most significant node k=",ind_minP_node),d=paste0("logCPM of m=",paste0(nodes_left,collapse = ",")),l=paste0("logCPM of m=",paste0(nodes_right,collapse = ","))) ),y="Log ratio of logCPM",x=x_label)

		}else{

		}


		# p3<-ggplot(data=df_dm, aes(x = groups, y = LogCPM, fill = groups)) +
		#   geom_boxplot() + theme_bw()+ theme(legend.position="bottom",plot.title = element_text(hjust = 0.5))+labs(title=paste0("Log ratio of log CPMs of test node k=",ind_minP_node),y="Log ratio of log CPM")

	}

	vp <- viewport(height = unit(0.6,"npc"), width=unit(0.4, "npc"), 
		 just = c("left","top"),   y = 1, x = 0.6)
	print(p3, vp = vp)
	
}


plot_phylo<-function (x, type = "phylogram", use.edge.length = TRUE, node.pos = NULL, 
                      show.tip.label = TRUE, show.node.label = FALSE, edge.color = "black", 
                      edge.width = 1, edge.lty = 1, font = 3, cex = par("cex"), 
                      adj = NULL, srt = 0, no.margin = FALSE, root.edge = FALSE, 
                      label.offset = 0, underscore = FALSE, x.lim = NULL, y.lim = NULL, 
                      direction = "rightwards", lab4ut = NULL, tip.color = "black", 
                      plot = TRUE, rotate.tree = 0, open.angle = 0, node.depth = 1, 
                      align.tip.label = FALSE, ...) 
{
  # browser()
  Ntip <- length(x$tip.label)
  if (Ntip < 2) {
    warning("found less than 2 tips in the tree")
    return(NULL)
  }
  .nodeHeight <- function(edge, Nedge, yy) .C(node_height, 
                                              as.integer(edge[, 1]), as.integer(edge[, 2]), as.integer(Nedge), 
                                              as.double(yy))[[4]]
  .nodeDepth <- function(Ntip, Nnode, edge, Nedge, node.depth) .C(node_depth, 
                                                                  as.integer(Ntip), as.integer(edge[, 1]), as.integer(edge[, 
                                                                                                                           2]), as.integer(Nedge), double(Ntip + Nnode), as.integer(node.depth))[[5]]
  .nodeDepthEdgelength <- function(Ntip, Nnode, edge, Nedge, 
                                   edge.length) .C(node_depth_edgelength, as.integer(edge[, 
                                                                                          1]), as.integer(edge[, 2]), as.integer(Nedge), as.double(edge.length), 
                                                   double(Ntip + Nnode))[[5]]
  Nedge <- dim(x$edge)[1]
  Nnode <- x$Nnode
  if (any(x$edge < 1) || any(x$edge > Ntip + Nnode)) 
    stop("tree badly conformed; cannot plot. Check the edge matrix.")
  ROOT <- Ntip + 1
  type <- match.arg(type, c("phylogram", "cladogram", "fan", 
                            "unrooted", "radial"))
  direction <- match.arg(direction, c("rightwards", "leftwards", 
                                      "upwards", "downwards"))
  if (is.null(x$edge.length)) {
    use.edge.length <- FALSE
  }
  else {
    if (use.edge.length && type != "radial") {
      tmp <- sum(is.na(x$edge.length))
      if (tmp) {
        warning(paste(tmp, "branch length(s) NA(s): branch lengths ignored in the plot"))
        use.edge.length <- FALSE
      }
    }
  }
  if (is.numeric(align.tip.label)) {
    align.tip.label.lty <- align.tip.label
    align.tip.label <- TRUE
  }
  else {
    if (align.tip.label) 
      align.tip.label.lty <- 3
  }
  if (align.tip.label) {
    if (type %in% c("unrooted", "radial") || !use.edge.length || 
        is.ultrametric(x)) 
      align.tip.label <- FALSE
  }
  if (type %in% c("unrooted", "radial") || !use.edge.length || 
      is.null(x$root.edge) || !x$root.edge) 
    root.edge <- FALSE
  phyloORclado <- type %in% c("phylogram", "cladogram")
  horizontal <- direction %in% c("rightwards", "leftwards")
  xe <- x$edge
  if (phyloORclado) {
    phyOrder <- attr(x, "order")
    if (is.null(phyOrder) || phyOrder != "cladewise") {
      x <- reorder(x)
      if (!identical(x$edge, xe)) {
        ereorder <- match(x$edge[, 2], xe[, 2])
        if (length(edge.color) > 1) {
          edge.color <- rep(edge.color, length.out = Nedge)
          edge.color <- edge.color[ereorder]
        }
        if (length(edge.width) > 1) {
          edge.width <- rep(edge.width, length.out = Nedge)
          edge.width <- edge.width[ereorder]
        }
        if (length(edge.lty) > 1) {
          edge.lty <- rep(edge.lty, length.out = Nedge)
          edge.lty <- edge.lty[ereorder]
        }
      }
    }
    yy <- numeric(Ntip + Nnode)
    TIPS <- x$edge[x$edge[, 2] <= Ntip, 2]
    yy[TIPS] <- 1:Ntip
  }
  z <- reorder(x, order = "postorder")
  if (phyloORclado) {
    if (is.null(node.pos)) 
      node.pos <- if (type == "cladogram" && !use.edge.length) 
        2
    else 1
    if (node.pos == 1) 
      yy <- .nodeHeight(z$edge, Nedge, yy)
    else {
      ans <- .C(node_height_clado, as.integer(Ntip), as.integer(z$edge[, 
                                                                       1]), as.integer(z$edge[, 2]), as.integer(Nedge), 
                double(Ntip + Nnode), as.double(yy))
      xx <- ans[[5]] - 1
      yy <- ans[[6]]
    }
    if (!use.edge.length) {
      if (node.pos != 2) 
        xx <- .nodeDepth(Ntip, Nnode, z$edge, Nedge, 
                         node.depth) - 1
      xx <- max(xx) - xx
    }
    else {
      xx <- .nodeDepthEdgelength(Ntip, Nnode, z$edge, Nedge, 
                                 z$edge.length)
    }
  }
  else {
    twopi <- 2 * pi
    rotate.tree <- twopi * rotate.tree/360
    if (type != "unrooted") {
      TIPS <- x$edge[which(x$edge[, 2] <= Ntip), 2]
      xx <- seq(0, twopi * (1 - 1/Ntip) - twopi * open.angle/360, 
                length.out = Ntip)
      theta <- double(Ntip)
      theta[TIPS] <- xx
      theta <- c(theta, numeric(Nnode))
    }
    switch(type, fan = {
      theta <- .nodeHeight(z$edge, Nedge, theta)
      if (use.edge.length) {
        r <- .nodeDepthEdgelength(Ntip, Nnode, z$edge, 
                                  Nedge, z$edge.length)
      } else {
        r <- .nodeDepth(Ntip, Nnode, z$edge, Nedge, node.depth)
        r <- 1/r
      }
      theta <- theta + rotate.tree
      if (root.edge) r <- r + x$root.edge
      xx <- r * cos(theta)
      yy <- r * sin(theta)
    }, unrooted = {
      nb.sp <- .nodeDepth(Ntip, Nnode, z$edge, Nedge, node.depth)
      XY <- if (use.edge.length) unrooted.xy(Ntip, Nnode, 
                                             z$edge, z$edge.length, nb.sp, rotate.tree) else unrooted.xy(Ntip, 
                                                                                                         Nnode, z$edge, rep(1, Nedge), nb.sp, rotate.tree)
      xx <- XY$M[, 1] - min(XY$M[, 1])
      yy <- XY$M[, 2] - min(XY$M[, 2])
    }, radial = {
      r <- .nodeDepth(Ntip, Nnode, z$edge, Nedge, node.depth)
      r[r == 1] <- 0
      r <- 1 - r/Ntip
      theta <- .nodeHeight(z$edge, Nedge, theta) + rotate.tree
      xx <- r * cos(theta)
      yy <- r * sin(theta)
    })
  }
  if (phyloORclado) {
    if (!horizontal) {
      tmp <- yy
      yy <- xx
      xx <- tmp - min(tmp) + 1
    }
    if (root.edge) {
      if (direction == "rightwards") 
        xx <- xx + x$root.edge
      if (direction == "upwards") 
        yy <- yy + x$root.edge
    }
  }
  if (no.margin) 
    par(mai = rep(0, 4))
  if (show.tip.label) 
    nchar.tip.label <- nchar(x$tip.label)
  max.yy <- max(yy)
  getLimit <- function(x, lab, sin, cex) {
    s <- strwidth(lab, "inches", cex = cex)
    if (any(s > sin)) 
      return(1.5 * max(x))
    Limit <- 0
    while (any(x > Limit)) {
      i <- which.max(x)
      alp <- x[i]/(sin - s[i])
      Limit <- x[i] + alp * s[i]
      x <- x + alp * s
    }
    Limit
  }
  if (is.null(x.lim)) {
    if (phyloORclado) {
      if (horizontal) {
        xx.tips <- xx[1:Ntip]
        if (show.tip.label) {
          pin1 <- par("pin")[1]
          tmp <- getLimit(xx.tips, x$tip.label, pin1, 
                          cex)
          tmp <- tmp + label.offset
        }
        else tmp <- max(xx.tips)
        x.lim <- c(0, tmp)
      }
      else x.lim <- c(1, Ntip)
    }
    else switch(type, fan = {
      if (show.tip.label) {
        offset <- max(nchar.tip.label * 0.018 * max.yy * 
                        cex)
        x.lim <- range(xx) + c(-offset, offset)
      } else x.lim <- range(xx)
    }, unrooted = {
      if (show.tip.label) {
        offset <- max(nchar.tip.label * 0.018 * max.yy * 
                        cex)
        x.lim <- c(0 - offset, max(xx) + offset)
      } else x.lim <- c(0, max(xx))
    }, radial = {
      if (show.tip.label) {
        offset <- max(nchar.tip.label * 0.03 * cex)
        x.lim <- c(-1 - offset, 1 + offset)
      } else x.lim <- c(-1, 1)
    })
  }
  else if (length(x.lim) == 1) {
    x.lim <- c(0, x.lim)
    if (phyloORclado && !horizontal) 
      x.lim[1] <- 1
    if (type %in% c("fan", "unrooted") && show.tip.label) 
      x.lim[1] <- -max(nchar.tip.label * 0.018 * max.yy * 
                         cex)
    if (type == "radial") 
      x.lim[1] <- if (show.tip.label) 
        -1 - max(nchar.tip.label * 0.03 * cex)
    else -1
  }
  if (phyloORclado && direction == "leftwards") 
    xx <- x.lim[2] - xx
  if (is.null(y.lim)) {
    if (phyloORclado) {
      if (horizontal) 
        y.lim <- c(1, Ntip)
      else {
        pin2 <- par("pin")[2]
        yy.tips <- yy[1:Ntip]
        if (show.tip.label) {
          tmp <- getLimit(yy.tips, x$tip.label, pin2, 
                          cex)
          tmp <- tmp + label.offset
        }
        else tmp <- max(yy.tips)
        y.lim <- c(0, tmp)
      }
    }
    else switch(type, fan = {
      if (show.tip.label) {
        offset <- max(nchar.tip.label * 0.018 * max.yy * 
                        cex)
        y.lim <- c(min(yy) - offset, max.yy + offset)
      } else y.lim <- c(min(yy), max.yy)
    }, unrooted = {
      if (show.tip.label) {
        offset <- max(nchar.tip.label * 0.018 * max.yy * 
                        cex)
        y.lim <- c(0 - offset, max.yy + offset)
      } else y.lim <- c(0, max.yy)
    }, radial = {
      if (show.tip.label) {
        offset <- max(nchar.tip.label * 0.03 * cex)
        y.lim <- c(-1 - offset, 1 + offset)
      } else y.lim <- c(-1, 1)
    })
  }
  else if (length(y.lim) == 1) {
    y.lim <- c(0, y.lim)
    if (phyloORclado && horizontal) 
      y.lim[1] <- 1
    if (type %in% c("fan", "unrooted") && show.tip.label) 
      y.lim[1] <- -max(nchar.tip.label * 0.018 * max.yy * 
                         cex)
    if (type == "radial") 
      y.lim[1] <- if (show.tip.label) 
        -1 - max(nchar.tip.label * 0.018 * max.yy * cex)
    else -1
  }
  if (phyloORclado && direction == "downwards") 
    yy <- y.lim[2] - yy
  if (phyloORclado && root.edge) {
    if (direction == "leftwards") 
      x.lim[2] <- x.lim[2] + x$root.edge
    if (direction == "downwards") 
      y.lim[2] <- y.lim[2] + x$root.edge
  }
  asp <- if (type %in% c("fan", "radial", "unrooted")) 
    1
  else NA
  plot.default(0, type = "n", xlim = x.lim, ylim = y.lim, xlab = "", 
               ylab = "", axes = FALSE, asp = asp, ...)
  if (plot) {
    if (is.null(adj)) 
      adj <- if (phyloORclado && direction == "leftwards") 
        1
    else 0
    if (phyloORclado && show.tip.label) {
      MAXSTRING <- max(strwidth(x$tip.label, cex = cex))
      loy <- 0
      if (direction == "rightwards") {
        lox <- label.offset + MAXSTRING * 1.05 * adj
      }
      if (direction == "leftwards") {
        lox <- -label.offset - MAXSTRING * 1.05 * (1 - 
                                                     adj)
      }
      if (!horizontal) {
        psr <- par("usr")
        MAXSTRING <- MAXSTRING * 1.09 * (psr[4] - psr[3])/(psr[2] - 
                                                             psr[1])
        loy <- label.offset + MAXSTRING * 1.05 * adj
        lox <- 0
        srt <- 90 + srt
        if (direction == "downwards") {
          loy <- -loy
          srt <- 180 + srt
        }
      }
    }
    if (type == "phylogram") {
      # if(dim(x$edge)[1]==2){
      #   tmp1<-x$edge[1,]
      #   tmp2<-x$edge[2,]
      #   x$edge[1,]<-tmp2
      #   x$edge[2,]<-tmp1
      # }
      phylogram.plot(x$edge, Ntip, Nnode, xx, yy, horizontal, 
                     edge.color, edge.width, edge.lty)
    }
    else {
      if (type == "fan") {
        ereorder <- match(z$edge[, 2], x$edge[, 2])
        if (length(edge.color) > 1) {
          edge.color <- rep(edge.color, length.out = Nedge)
          edge.color <- edge.color[ereorder]
        }
        if (length(edge.width) > 1) {
          edge.width <- rep(edge.width, length.out = Nedge)
          edge.width <- edge.width[ereorder]
        }
        if (length(edge.lty) > 1) {
          edge.lty <- rep(edge.lty, length.out = Nedge)
          edge.lty <- edge.lty[ereorder]
        }
        circular.plot(z$edge, Ntip, Nnode, xx, yy, theta, 
                      r, edge.color, edge.width, edge.lty)
      }
      else cladogram.plot(x$edge, xx, yy, edge.color, edge.width, 
                          edge.lty)
    }
    if (root.edge) {
      rootcol <- if (length(edge.color) == 1) 
        edge.color
      else "black"
      rootw <- if (length(edge.width) == 1) 
        edge.width
      else 1
      rootlty <- if (length(edge.lty) == 1) 
        edge.lty
      else 1
      if (type == "fan") {
        tmp <- polar2rect(x$root.edge, theta[ROOT])
        segments(0, 0, tmp$x, tmp$y, col = rootcol, lwd = rootw, 
                 lty = rootlty)
      }
      else {
        switch(direction, rightwards = segments(0, yy[ROOT], 
                                                x$root.edge, yy[ROOT], col = rootcol, lwd = rootw, 
                                                lty = rootlty), leftwards = segments(xx[ROOT], 
                                                                                     yy[ROOT], xx[ROOT] + x$root.edge, yy[ROOT], 
                                                                                     col = rootcol, lwd = rootw, lty = rootlty), 
               upwards = segments(xx[ROOT], 0, xx[ROOT], x$root.edge, 
                                  col = rootcol, lwd = rootw, lty = rootlty), 
               downwards = segments(xx[ROOT], yy[ROOT], xx[ROOT], 
                                    yy[ROOT] + x$root.edge, col = rootcol, lwd = rootw, 
                                    lty = rootlty))
      }
    }
    if (show.tip.label) {
      if (is.expression(x$tip.label)) 
        underscore <- TRUE
      if (!underscore) 
        x$tip.label <- gsub("_", " ", x$tip.label)
      if (phyloORclado) {
        if (align.tip.label) {
          xx.tmp <- switch(direction, rightwards = max(xx[1:Ntip]), 
                           leftwards = min(xx[1:Ntip]), upwards = xx[1:Ntip], 
                           downwards = xx[1:Ntip])
          yy.tmp <- switch(direction, rightwards = yy[1:Ntip], 
                           leftwards = yy[1:Ntip], upwards = max(yy[1:Ntip]), 
                           downwards = min(yy[1:Ntip]))
          segments(xx[1:Ntip], yy[1:Ntip], xx.tmp, yy.tmp, 
                   lty = align.tip.label.lty)
        }
        else {
          xx.tmp <- xx[1:Ntip]
          yy.tmp <- yy[1:Ntip]
        }
        text(xx.tmp + lox, yy.tmp + loy, x$tip.label, 
             adj = adj, font = font, srt = srt, cex = cex, 
             col = tip.color)
      }
      else {
        angle <- if (type == "unrooted") 
          XY$axe
        else atan2(yy[1:Ntip], xx[1:Ntip])
        lab4ut <- if (is.null(lab4ut)) {
          if (type == "unrooted") 
            "horizontal"
          else "axial"
        }
        else match.arg(lab4ut, c("horizontal", "axial"))
        xx.tips <- xx[1:Ntip]
        yy.tips <- yy[1:Ntip]
        if (label.offset) {
          xx.tips <- xx.tips + label.offset * cos(angle)
          yy.tips <- yy.tips + label.offset * sin(angle)
        }
        if (lab4ut == "horizontal") {
          y.adj <- x.adj <- numeric(Ntip)
          sel <- abs(angle) > 0.75 * pi
          x.adj[sel] <- -strwidth(x$tip.label)[sel] * 
            1.05
          sel <- abs(angle) > pi/4 & abs(angle) < 0.75 * 
            pi
          x.adj[sel] <- -strwidth(x$tip.label)[sel] * 
            (2 * abs(angle)[sel]/pi - 0.5)
          sel <- angle > pi/4 & angle < 0.75 * pi
          y.adj[sel] <- strheight(x$tip.label)[sel]/2
          sel <- angle < -pi/4 & angle > -0.75 * pi
          y.adj[sel] <- -strheight(x$tip.label)[sel] * 
            0.75
          text(xx.tips + x.adj * cex, yy.tips + y.adj * 
                 cex, x$tip.label, adj = c(adj, 0), font = font, 
               srt = srt, cex = cex, col = tip.color)
        }
        else {
          if (align.tip.label) {
            POL <- rect2polar(xx.tips, yy.tips)
            POL$r[] <- max(POL$r)
            REC <- polar2rect(POL$r, POL$angle)
            xx.tips <- REC$x
            yy.tips <- REC$y
            segments(xx[1:Ntip], yy[1:Ntip], xx.tips, 
                     yy.tips, lty = align.tip.label.lty)
          }
          if (type == "unrooted") {
            adj <- abs(angle) > pi/2
            angle <- angle * 180/pi
            angle[adj] <- angle[adj] - 180
            adj <- as.numeric(adj)
          }
          else {
            s <- xx.tips < 0
            angle <- angle * 180/pi
            angle[s] <- angle[s] + 180
            adj <- as.numeric(s)
          }
          font <- rep(font, length.out = Ntip)
          tip.color <- rep(tip.color, length.out = Ntip)
          cex <- rep(cex, length.out = Ntip)
          for (i in 1:Ntip) text(xx.tips[i], yy.tips[i], 
                                 x$tip.label[i], font = font[i], cex = cex[i], 
                                 srt = angle[i], adj = adj[i], col = tip.color[i])
        }
      }
    }
    if (show.node.label) 
      text(xx[ROOT:length(xx)] + label.offset, yy[ROOT:length(yy)], 
           x$node.label, adj = adj, font = font, srt = srt, 
           cex = cex)
  }
  L <- list(type = type, use.edge.length = use.edge.length, 
            node.pos = node.pos, node.depth = node.depth, show.tip.label = show.tip.label, 
            show.node.label = show.node.label, font = font, cex = cex, 
            adj = adj, srt = srt, no.margin = no.margin, label.offset = label.offset, 
            x.lim = x.lim, y.lim = y.lim, direction = direction, 
            tip.color = tip.color, Ntip = Ntip, Nnode = Nnode, root.time = x$root.time, 
            align.tip.label = align.tip.label)
  assign("last_plot.phylo", c(L, list(edge = xe, xx = xx, yy = yy)), 
         envir = .PlotPhyloEnv)
  invisible(L)
  return(list(xx,yy))
}


makeCovMat<-function(covariates,data){
      cov_mat<-data[covariates]
      indList_num<-which(sapply(cov_mat,is.numeric))
      indList_factor<-which(!sapply(cov_mat,is.numeric))
      covCateList<-lapply(indList_factor,function(ind_factor){
        fastDummies::dummy_cols(cov_mat[ind_factor])[,-1]
      })
      covCate<-do.call("cbind",covCateList)
      if(is.null(covCate)){
        covResult<-cov_mat[indList_num]
      }else{
        covResult<-cbind(cov_mat[indList_num],covCate)  
      }
      return(covResult)
}