In Dmirandae/blepd: Branch Length Evaluation in Phylogenetic Diversity Analysis

This is a worked example using a four terminals tree, distributed in two areas. Initially, all branches are equal.
# Four taxa and two areas


## Load libraries

library(blepd)
library(ggtree)
library(ggplot2)



## Call the data

# Call the four terminals example

##data(package = "blepd")

## trees

data(tree)


## Plot the tree

plot(tree,use.edge.length =T)

initialTree <- tree


## Load and plot the distribution

data(distribution)

##str(distribution)

dist4taxa <- distribution

## Change distribution to XY to be able to plot the data, otherwise, skip the next step

distXY <- matrix2XY(dist4taxa)

distXY

## plotting

## the tree, using ggtree

plotTree <-  ggtree(initialTree, ladderize=TRUE,
                    color="black", size=0.5, linetype="dotted")  +
             geom_tiplab(size=4.5, color="black")                +
             theme_tree2()                                       +
             ggtitle("Four terminals, all branches are equal to 1")



print(plotTree)



## The distribution using ggplot


plotDistrib <- ggplot(data=distXY,
                      aes(x= Area, y= Terminal , size =150))    +
               geom_point()                                     +
               theme_bw()                                       +
               theme( axis.text.y = element_blank(),
                      panel.grid.major = element_blank(),  
                      panel.grid.minor = element_blank(),  
                      axis.line = element_line(colour = NA_character_),  
                      axis.ticks = element_blank(),  
                      legend.position = "none" )               +
               labs(title = "Distributions",
                    y = "",
                    x = "Area")                                 


##
print(plotDistrib)

cowplot::plot_grid(plotTree, plotDistrib, ncol=2)



# Check whether names in both objects, initialTree and dist4taxa, are the same.


##all(colnames(dist4taxa) == initialTree$tip.label)

all(colnames(dist4taxa) %in% initialTree$tip.label)


#Report the branch length, and calculate the PD values.


initialTree$edge.length

initialPD <- PDindex(tree=initialTree, distribution = dist4taxa)

initialPD



#As expected, there is a tie between both areas.


#Now, we can see the impact of swapping the branch lengths; we will use the three
#available models to swap: "simpleswap", "allswap", "uniform", and we will swap the
#three classes of branches: "terminals", "internals", "all".


##?swapBL


for( modelo in c("simpleswap","allswap","uniform") ){

        for( rama in c("terminals","internals","all") ){

 cat( "\n\t Branchs swapped=\t",rama,".\t",sep="")
            val <- swapBL(tree=initialTree,
                          distribution = dist4taxa,
                          model = modelo,
                          branch = rama
                         )
        printSwapBL(val, compact=TRUE)

         cat( "\n\n" )

     }

  }



#As expected, as all braches are EQUAL, nothing happens


#Now, let's see if the results could be related to the branch length (or not).



## we could use a single command (*evalBranch*)


##?evalBranch

AllTerminals <- evalBranch(tree=initialTree,
                          distribution = dist4taxa,
                          branchToEval = "terminals", 
                          approach="all")


## The structure of the output

##str(AllTerminals)

## Print it as a friendly output

printEvalBranch(AllTerminals, compact=FALSE)

## Now the internals

AllInternals <- evalBranch(tree=initialTree,
                          distribution = dist4taxa,
                          branchToEval = "internals", 
                          approach="all")


printEvalBranch(AllInternals, compact=FALSE)




#Now with a long branch


initialTree <- tree

initialTree$edge.length[1] <- 5

plotTree <-  ggtree(initialTree, ladderize=TRUE,
                    color="black", size=0.5, linetype="dotted")    +
             geom_tiplab(size=4.5, color="black")                  +
             theme_tree2()                                         +
             ggtitle(paste("Four terminals. Branches 1=",initialTree$edge.length[1],
             ".  2=",initialTree$edge.length[2],".",sep=""))


##
print(plotTree)


 for( modelo in c("simpleswap","allswap","uniform") ){

        for( rama in c("terminals","internals","all") ){

           cat( "\n\t Branchs swapped=\t",rama,".\t",sep="")

            val <- swapBL(tree=initialTree,
                          distribution = dist4taxa,
                          model = modelo,
                          branch = rama
                         )

        printSwapBL(val, compact=TRUE)

         cat( "\n\n" )

     }

  }




AllTerminals <- evalBranch(tree=initialTree,
                          distribution = dist4taxa,
                          branchToEval = "terminals", 
                          approach="all")

printEvalBranch(AllTerminals, compact=FALSE)

## Now the internals

AllInternals <- evalBranch(tree=initialTree,
                          distribution = dist4taxa,
                          branchToEval = "internals", 
                          approach="all")


printEvalBranch(AllInternals, compact=FALSE)




# Now with a longer branch


initialTree <- tree

initialTree$edge.length[1] <- 10


plotTree <-  ggtree(initialTree, ladderize=TRUE,
                    color="black", size=0.5, linetype="dotted")       +
             geom_tiplab(size=4.5, color="black")                     +
             theme_tree2()                                            +
             ggtitle(paste("Four terminals. Branches 1=",initialTree$edge.length[1],
             ".  2=",initialTree$edge.length[2],".",sep=""))


##
print(plotTree)


 for( modelo in c("simpleswap","allswap","uniform") ){

        for( rama in c("terminals","internals","all") ){

 cat( "\n\t Branchs swapped=\t",rama,".\t",sep="")
            val <- swapBL(tree=initialTree,
                          distribution = dist4taxa,
                          model = modelo,
                          branch = rama
                         )
        printSwapBL(val, compact=TRUE)

         cat( "\n\n" )

     }

  }




AllTerminals <- evalBranch(tree=initialTree,
                          distribution = dist4taxa,
                          branchToEval = "terminals", 
                          approach="all")

printEvalBranch(AllTerminals, compact=FALSE)

## Now the internals

AllInternals <- evalBranch(tree=initialTree,
                          distribution = dist4taxa,
                          branchToEval = "internals", 
                          approach="all")


printEvalBranch(AllInternals, compact=FALSE)




# With a short branch



initialTree <- tree

initialTree$edge.length[1] <- 0.1


plotTree <-  ggtree(initialTree, ladderize=TRUE,
                    color="black", size=0.5, linetype="dotted")  +
             geom_tiplab(size=4.5, color="black")                +
             theme_tree2()                                       +
             ggtitle(paste("Four terminals. Branches 1=",initialTree$edge.length[1],
             ".  2=",initialTree$edge.length[2],".",sep=""))


##
print(plotTree)


 for( modelo in c("simpleswap","allswap","uniform") ){

        for( rama in c("terminals","internals","all") ){

 cat( "\n\t Branchs swapped=\t",rama,".\t",sep="")

            val <- swapBL(tree=initialTree,
                          distribution = dist4taxa,
                          model = modelo,
                          branch = rama
                         )

        printSwapBL(val, compact=TRUE)

         cat( "\n\n" )

     }

  }




AllTerminals <- evalBranch(tree=initialTree,
                          distribution = dist4taxa,
                          branchToEval = "terminals", 
                          approach="all")

printEvalBranch(AllTerminals, compact=FALSE)

## Now the internals

AllInternals <- evalBranch(tree=initialTree,
                          distribution = dist4taxa,
                          branchToEval = "internals", 
                          approach="all")


printEvalBranch(AllInternals, compact=FALSE)




#With a long branch in branch 2


initialTree <- tree

initialTree$edge.length[2] <- 5

plotTree <-  ggtree(initialTree, ladderize=TRUE,
                    color="black", size=0.5, linetype="dotted")  +
             geom_tiplab(size=4.5, color="black")                +
             theme_tree2()                                       +
             ggtitle(paste("Four terminals. Branches 1=",initialTree$edge.length[1],
             ".  2=",initialTree$edge.length[2],".",sep=""))


##
print(plotTree)


 for( modelo in c("simpleswap","allswap","uniform") ){

        for( rama in c("terminals","internals","all") ){

 cat( "\n\t Branchs swapped=\t",rama,".\t",sep="")

            val <- swapBL(tree=initialTree,
                          distribution = dist4taxa,
                          model = modelo,
                          branch = rama
                         )

        printSwapBL(val, compact=TRUE)

         cat( "\n\n" )

     }

  }




AllTerminals <- evalBranch(tree=initialTree,
                          distribution = dist4taxa,
                          branchToEval = "terminals", 
                          approach="all")

printEvalBranch(AllTerminals, compact=FALSE)

## Now the internals

AllInternals <- evalBranch(tree=initialTree,
                          distribution = dist4taxa,
                          branchToEval = "internals", 
                          approach="all")


printEvalBranch(AllInternals, compact=FALSE)




# Now with a longer branch2


initialTree <- tree

initialTree$edge.length[2] <- 10


plotTree <-  ggtree(initialTree, ladderize=TRUE,
                    color="black", size=0.5, linetype="dotted")  +
             geom_tiplab(size=4.5, color="black")                +
             theme_tree2()                                       +
             ggtitle(paste("Four terminals. Branches 1=",initialTree$edge.length[1],
             ".  2=",initialTree$edge.length[2],".",sep=""))


##
print(plotTree)


 for( modelo in c("simpleswap","allswap","uniform") ){

        for( rama in c("terminals","internals","all") ){

 cat( "\n\t Branchs swapped=\t",rama,".\t",sep="") 

            val <- swapBL(tree=initialTree,
                          distribution = dist4taxa,
                          model = modelo,
                          branch = rama
                         )

        printSwapBL(val, compact=TRUE)

         cat( "\n\n" )

     }

  }




AllTerminals <- evalBranch(tree=initialTree,
                          distribution = dist4taxa,
                          branchToEval = "terminals", 
                          approach="all")

printEvalBranch(AllTerminals, compact=FALSE)

## Now the internals

AllInternals <- evalBranch(tree=initialTree,
                          distribution = dist4taxa,
                          branchToEval = "internals", 
                          approach="all")


printEvalBranch(AllInternals, compact=FALSE)




#now with a shorter branch



initialTree <- tree

initialTree$edge.length[2] <- 0.1


plotTree <-  ggtree(initialTree, ladderize=TRUE,
                    color="black", size=0.5, linetype="dotted")             +
             geom_tiplab(size=4.5, color="black")                           +
             theme_tree2()                                                  +
             ggtitle(paste("Four terminals. Branches 1=",initialTree$edge.length[1],
             ".  2=",initialTree$edge.length[2],".",sep=""))


##
print(plotTree)


 for( modelo in c("simpleswap","allswap","uniform") ){

        for( rama in c("terminals","internals","all") ){

 cat( "\n\t Branchs swapped=\t",rama,".\t",sep="")

            val <- swapBL(tree=initialTree,
                          distribution = dist4taxa,
                          model = modelo,
                          branch = rama
                         )

        printSwapBL(val, compact=TRUE)

         cat( "\n\n" )

     }

  }




AllTerminals <- evalBranch(tree=initialTree,
                          distribution = dist4taxa,
                          branchToEval = "terminals", 
                          approach="all")

printEvalBranch(AllTerminals, compact=FALSE)

## Now the internals

AllInternals <- evalBranch(tree=initialTree,
                          distribution = dist4taxa,
                          branchToEval = "internals", 
                          approach="all")


printEvalBranch(AllInternals, compact=FALSE)