##Created by Federico Luebert
##Version 7 November 2023

library(ape)

##create topology with 5 tips and read it into a multiPhylo object
txt<-"(((A,B),(C,D)),E);"
top<-list()
top<-lapply(1:1000,read.tree,text=txt)
class(top)<-"multiPhylo"

## add random branch lengths between 0 and 100 to each topology
topol<-lapply(top, function(x) compute.brlen(x, method=runif(8,0,100)))
class(topol)<-"multiPhylo"

#replicate original topology and branch lengths 
topol1<-topol
topol2<-topol


##change toplogy for reploicated ones, 1=shallow and 2=deep 
for(i in 1:1000){
    topol1[[i]]$tip.label<-c("C","B","A","D","E")
    topol2[[i]]$tip.label<-c("E","B","C","D","A")
}    

##make ultrametric trees of each topology
topolultra<-lapply(1:1000,function(x) chronos(topol[[x]],lambda=0.01))
topol1ultra<-lapply(1:1000,function(x) chronos(topol1[[x]],lambda=0.01))
topol2ultra<-lapply(1:1000,function(x) chronos(topol2[[x]],lambda=0.01))

#install.packages("picante")
library(picante)

##create community data matrix
comm<-matrix(nrow=5,ncol=5)
colnames(comm)<-c("A","B","C","D","E")
rownames(comm)<-1:5
comm[1,]<-c(1,1,0,0,0)
comm[2,]<-c(1,0,1,1,1)
comm[3,]<-c(1,1,0,1,1)
comm[4,]<-c(0,1,1,1,0)
comm[5,]<-c(0,0,1,1,1)

##claculate PD for phylograms
pd<-lapply(1:1000, function(x) pd(comm,topol[[x]]))
pd1<-lapply(1:1000, function(x) pd(comm,topol1[[x]]))
pd2<-lapply(1:1000, function(x) pd(comm,topol2[[x]]))

##claculate PD for ultrametric trees
pdu<-lapply(1:1000, function(x) pd(comm,topolultra[[x]]))
pd1u<-lapply(1:1000, function(x) pd(comm,topol1ultra[[x]]))
pd2u<-lapply(1:1000, function(x) pd(comm,topol2ultra[[x]]))

##differences in max and min ranking values for phylograms and chronograms

tab55<-matrix(nrow=6,ncol=5)#matrix for storing differences in in and max PD ranking for five sites in each six versions of the tree across 1000 simulated trees with five tips  
colnames(tab55)<-c(1:5)
rownames(tab55)<-c("phylograms original","chronogram original","phylograms shallow","chronogram shallow","phylograms deep","chronogram deep")

#phylograms original
tab55[1,]<-abs(apply(sapply(1:1000, function(x) rank(pd[[x]]$PD)),1, min)-apply(sapply(1:1000, function(x) rank(pd[[x]]$PD)),1, max))

#chronogram original
tab55[2,]<-abs(apply(sapply(1:1000, function(x) rank(pdu[[x]]$PD)),1, min)-apply(sapply(1:1000, function(x) rank(pdu[[x]]$PD)),1, max))

#phylograms shallow
tab55[3,]<-abs(apply(sapply(1:1000, function(x) rank(pd1[[x]]$PD)),1, min)-apply(sapply(1:1000, function(x) rank(pd1[[x]]$PD)),1, max))

#chronogram shallow
tab55[4,]<-abs(apply(sapply(1:1000, function(x) rank(pd1u[[x]]$PD)),1, min)-apply(sapply(1:1000, function(x) rank(pd1u[[x]]$PD)),1, max))

#phylograms deep
tab55[5,]<-abs(apply(sapply(1:1000, function(x) rank(pd2[[x]]$PD)),1, min)-apply(sapply(1:1000, function(x) rank(pd2[[x]]$PD)),1, max))

#chronograms deep
tab55[6,]<-abs(apply(sapply(1:1000, function(x) rank(pd2u[[x]]$PD)),1, min)-apply(sapply(1:1000, function(x) rank(pd2u[[x]]$PD)),1, max))

##10 tips
##create topology with 10 tips and read it into a multiPhylo object
txt10<-"((((A,B),(C,D)),E),(F,(G,(H,(I,J)))));"
top10<-list()
top10<-lapply(1:1000,read.tree,text=txt10)
class(top10)<-"multiPhylo"

## add random branch lengths between 0 and 100 to each topology
topol10<-lapply(top10, function(x) compute.brlen(x, method=runif(18,0,100)))
class(topol)<-"multiPhylo"

#replicate original topology and branch lengths 
topol10.1<-topol10
topol10.2<-topol10


##change toplogy for reploicated ones, 1=shallow (in this case I also swapt J and H) and 2=deep (in this case I also swapt D and I) 
for(i in 1:1000){
    topol10.1[[i]]$tip.label<-c("C","B","A","D","E","F","G","J","I","H")
    topol10.2[[i]]$tip.label<-c("E","B","C","I","A","F","G","H","D","J")
}    

##make ultrametric trees of each topology
topol10ultra<-lapply(1:1000,function(x) chronos(topol10[[x]],lambda=0.01))
topol10.1ultra<-lapply(1:1000,function(x) chronos(topol10.1[[x]],lambda=0.01))
topol10.2ultra<-lapply(1:1000,function(x) chronos(topol10.2[[x]],lambda=0.01))

##create community data matrix
comm10<-matrix(nrow=5,ncol=10)
colnames(comm10)<-c("A","B","C","D","E","F","G","H","I","J")
rownames(comm10)<-1:5
comm10[1,]<-c(1,1,0,0,0,0,0,0,1,1)
comm10[2,]<-c(1,0,1,1,1,1,1,1,0,1)
comm10[3,]<-c(1,1,0,1,1,1,1,0,1,1)
comm10[4,]<-c(0,1,1,1,0,0,1,1,1,0)
comm10[5,]<-c(0,0,1,1,1,1,1,1,0,0)

##claculate PD for phylograms
pd10<-lapply(1:1000, function(x) pd(comm10,topol10[[x]]))
pd10.1<-lapply(1:1000, function(x) pd(comm10,topol10.1[[x]]))
pd10.2<-lapply(1:1000, function(x) pd(comm10,topol10.2[[x]]))

##claculate PD for ultrametric trees
pd10u<-lapply(1:1000, function(x) pd(comm10,topol10ultra[[x]]))
pd10.1u<-lapply(1:1000, function(x) pd(comm10,topol10.1ultra[[x]]))
pd10.2u<-lapply(1:1000, function(x) pd(comm10,topol10.2ultra[[x]]))

##differences in max and min ranking values for phylograms and chronograms

tab510<-matrix(nrow=6,ncol=5)#matrix for storing differences in in and max PD ranking for five sites in each six versions of the tree across 1000 simulated trees with 10 tips  
colnames(tab510)<-c(1:5)
rownames(tab510)<-c("phylograms original","chronogram original","phylograms shallow","chronogram shallow","phylograms deep","chronogram deep")

#phylograms original
tab510[1,]<-abs(apply(sapply(1:1000, function(x) rank(pd10[[x]]$PD)),1, min)-apply(sapply(1:1000, function(x) rank(pd10[[x]]$PD)),1, max))

#chronogram original
tab510[2,]<-abs(apply(sapply(1:1000, function(x) rank(pd10u[[x]]$PD)),1, min)-apply(sapply(1:1000, function(x) rank(pd10u[[x]]$PD)),1, max))

#phylograms shallow
tab510[3,]<-abs(apply(sapply(1:1000, function(x) rank(pd10.1[[x]]$PD)),1, min)-apply(sapply(1:1000, function(x) rank(pd10.1[[x]]$PD)),1, max))

#chronogram shallow
tab510[4,]<-abs(apply(sapply(1:1000, function(x) rank(pd10.1u[[x]]$PD)),1, min)-apply(sapply(1:1000, function(x) rank(pd10.1u[[x]]$PD)),1, max))

#phylograms deep
tab510[5,]<-abs(apply(sapply(1:1000, function(x) rank(pd10.2[[x]]$PD)),1, min)-apply(sapply(1:1000, function(x) rank(pd10.2[[x]]$PD)),1, max))

#chronograms deep
tab510[6,]<-abs(apply(sapply(1:1000, function(x) rank(pd10.2u[[x]]$PD)),1, min)-apply(sapply(1:1000, function(x) rank(pd10.2u[[x]]$PD)),1, max))

write.table(tab55,"5tips_simul.txt",quote=F,sep="\t")#simulation with five tips
write.table(tab510,"10tips_simul.txt",quote=F,sep="\t")#simulation with ten  tips
write.table(comm10,"comm10_simul.txt",quote=F,sep="\t")#community table