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

Four taxa and two areas

Now, we load the data included in the package: tree and distribution.

## Libraries

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

## Data

data(package = "blepd")

## trees

data(tree)

str(tree)

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

initialTree <- tree


## distributions

data(distribution)

str(distribution)

dist4taxa <- distribution

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

distXY <- matrix2XY(dist4taxa)

distXY

## plotting

## the tree

plotTree <-  ggtree(initialTree, ladderize=TRUE,
                    color="black", size=0.5, linetype="dotted") +
             geom_tiplab(size=6, color="black") +
             theme_tree2() +
             ggtitle("Four terminals, equal branch length")


##print(plotTree)



## the distribution


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(),  # Remove major gridlines
                      panel.grid.minor = element_blank(),  # Remove minor gridlines
                      axis.line = element_line(colour = NA_character_),  # Remove axis lines
                      axis.ticks = element_blank(),  # Remove axis ticks
                      legend.position = "none"
                      ) +
               labs(title = "Distributions",
                    y = "",
                    x = "Area")                                 


##print(plotDistrib)

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

We 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)

We 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 the posible null model(s), 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") ){

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

         cat( "\n\t Model=",modelo,
             "\n\t Branchs swapped=",rama,"\n" )

        printBlepd(val)

         cat( "\n\n\n" )

     }

  }

As expected, as all braches are EQUAL notinh happens

Now, let us see if the results could be assigned to the longest branches (or not).

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


?evalBranch

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

print.evalBranchAll(AllTerminals)

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

print.evalBranchAll(AllTerminals)


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

print.evalBranchAll(AllInternals)