inst/doc/IntertwiningTreesAndNetworks.R

In phangorn: Phylogenetic Reconstruction and Analysis

## ----setup, echo=FALSE--------------------------------------------------------
# set global chunk options: images will be bigger
  knitr::opts_chunk$set(fig.width=6, fig.height=6)
#, global.par=TRUE
options(digits = 4)

knitr::knit_hooks$set(small.mar=function(before, options, envir){
   if (before && options$fig.show!='none') par(mar=c(.5,.5,.5,.5))
})

## ----eval=FALSE---------------------------------------------------------------
#  install.packages("phangorn", dependencies=TRUE)
#  # install latest development version needs devtools
#  install.packages("devtools", dependencies=TRUE)
#  library(devtools)
#  install_github("KlausVigo/phangorn")

## -----------------------------------------------------------------------------
library(phangorn)    # load the phangorn library

## -----------------------------------------------------------------------------
## automatically set the correct working directory for the examples below
# setwd(system.file("extdata/trees", package = "phangorn"))
# for this vignette we create a path to the files we want to load 
fdir <- system.file("extdata/trees", package = "phangorn")
## in your case it may look something like this...
# setwd("C:/TreesNetworks/Example Files")

##DNA Matrix, maybe not needed 
woodmouse <- read.phyDat(file.path(fdir, "woodmouse.fasta"),format="fasta") 
## RAxML best-known tree with bipartition support (from previous analysis)
raxml.tree <- read.tree(file.path(fdir,"RAxML_bipartitions.woodmouse"))
## RAxML bootstrap trees (from previous analysis)
raxml.bootstrap <- read.tree(file.path(fdir,"RAxML_bootstrap.woodmouse"))
## MrBayes consensus tree (50% majority rule) (from previous analysis)
mrbayes.tree <- read.nexus(file.path(fdir,"woodmouse.mrbayes.nex.con"))
## MrBayes sample runs 1 and 2 (from previous analysis)
run1 <- read.nexus(file.path(fdir,"woodmouse.mrbayes.nex.run1.t"))
run2 <- read.nexus(file.path(fdir,"woodmouse.mrbayes.nex.run2.t"))
## How many trees are in the MrBayes tree sample?
run1
run2
## Combining the two runs and removing 25% burn-in
mrbayes.trees <- c(run1[251:1001],run2[251:1001])
## NeigbourNet Nexus file generated by SplitsTree (from previous analysis)
Nnet <- read.nexus.networx(file.path(fdir,"woodmouse.nxs"))

## -----------------------------------------------------------------------------
par(mfrow=c(1,2), mar=c(1,1,1,1)) # Setting plot parameters
### Plotting trees with support values:
##  RAxML
plot(raxml.tree)
nodelabels(raxml.tree$node.label, adj = c(1, 0), frame = "none")
##  MrBayes
plot(mrbayes.tree)
nodelabels(mrbayes.tree$node.label, adj = c(1, 0), frame = "none")

par(mfrow=c(1,1)) # Setting plot parameters
# NeighborNet
plot(Nnet,"2D")
## alternatively,
# plot(Nnet,"2D")


## ----fig.width=7, fig.height=4, small.mar=TRUE--------------------------------
# create a vector of labels for the network corresponding to edges in the tree
edge.lab <- createLabel(Nnet, raxml.tree, raxml.tree$edge[,2], "edge")
# could be also 1:27 instead of raxml.tree$edge[,2]

# Show the correspondingly labelled tree and network in R
par(mfrow=c(1,2))  
plot(raxml.tree, "u", rotate.tree = 180, cex=.7) 
edgelabels(raxml.tree$edge[,2],col="blue", frame="none", cex=.7)

# find edges that are in the network but not in the tree
edge.col <- rep("black", nrow(Nnet$edge))
edge.col[ is.na(edge.lab) ] <- "red"
# or a simpler alternative...
edge.col <- createLabel(Nnet, raxml.tree, "black", nomatch="red")

x <- plot(Nnet, edge.label = edge.lab, show.edge.label = T, "2D", 
          edge.color = edge.col, col.edge.label = "blue", cex=.7)
# the above plot function returns an invisible networx object and this object 
# also contains the colors for the edges.


## ----small.mar=TRUE-----------------------------------------------------------
# the scaler argument multiplies the confidence values. This is useful to switch
# confidences values between total, percentage or ratios.   
x <- addConfidences(Nnet,raxml.tree, scaler = .01)
# find splits that are in the network but not in the tree
split.col <- rep("black", length(x$splits))
split.col[ !matchSplits(as.splits(x), as.splits(raxml.tree)) ] <- "red"

# simpler alternative...
split.col2 <- createLabel(x, raxml.tree, label="black", "split", nomatch="red")

# Plotting in R
out.x <- plot(x,"2D",show.edge.label=TRUE, split.color=split.col, 
              col.edge.label = "blue")

## -----------------------------------------------------------------------------
# write.nexus.networx(out.x,"woodmouse.tree.support.nxs")
## or we can also export the splits alone (for usage in software other than SplitsTree)
# write.nexus.splits(as.splits(out.x),"woodmouse.splits.support.nxs")

## ----small.mar=TRUE-----------------------------------------------------------
y <- addConfidences(Nnet, as.splits(raxml.bootstrap))
edge.col <- createLabel(y, raxml.tree, label="black", "edge", nomatch="grey")

y <- plot(y,"2D",show.edge.label=TRUE, edge.color=edge.col)

## Write to SplitsTree for viewing
# write.nexus.networx(y,"NN.with.bs.support.nxs")

## ----small.mar=TRUE-----------------------------------------------------------
cnet <- consensusNet(raxml.bootstrap,prob=0.10)
edge.col <- createLabel(cnet, Nnet, label="black", "edge", nomatch="grey")
cnet <- plot(cnet, "2D", show.edge.label = TRUE, edge.color=edge.col)

edge.col <- createLabel(Nnet, cnet, label="black", "edge", nomatch="grey")
z <- plot(Nnet, "2D", show.edge.label = TRUE, edge.color=edge.col)

obj <- addConfidences(Nnet,cnet)
plot(obj,"2D",show.edge.label=T, edge.color=edge.col, col.edge.label = "blue")

## Write to SplitsTree for viewing
# write.nexus.networx(obj,"Nnet.with.ML.Cnet.Bootstrap.nxs")

## ----fig.width=7, fig.height=6------------------------------------------------
Nnet <- read.nexus.networx(file.path(fdir,"RAxML_distances.Wang.nxs"))
raxml.tree <- read.tree(file.path(fdir,"RAxML_bestTree.Wang.out")) |> unroot()
raxml.bootstrap <- read.tree(file.path(fdir,"RAxML_bootstrap.Wang.out"))

bs_splits <- as.splits(raxml.bootstrap)
tree_splits <- as.splits(raxml.tree) |> unique() |> removeTrivialSplits()
# we overwrite bootstrap values and set the weights 
# to 1e-6 (almost zero), as we plot them on a log scale later on
attr(bs_splits, "weights")[] <- 1e-6
# combine the splits from the bootstrap and neighbor net
# and delete duplicates and add the confidence values
# we get rid of trivial splits
all_splits <- c(Nnet$splits, bs_splits) |> unique() |> removeTrivialSplits() |> 
  addConfidences(bs_splits, scaler=100)

# For easier plotting we create a matrix with the confidences and 
# weights as columns
tab <- data.frame(SplitWeight = attr(all_splits, "weights"), 
                  Bootstrap=attr(all_splits, "confidences"), Tree=FALSE)
# we add a logical variable pto indicate which splits are in the RAxML tree
tab$Tree[matchSplits(tree_splits, all_splits, FALSE)] <- TRUE

tab[is.na(tab[,"Bootstrap"]),"Bootstrap"] <- 0
tab[,"Bootstrap"] <- round(tab[,"Bootstrap"])
rownames(tab) <- apply(as.matrix(all_splits, zero.print = ".", one.print = "|"), 
                       1, paste0, collapse="")
tab[1:10,]

col <- rep("blue", nrow(tab))
col[tab[,"Bootstrap"]==0] <- "green"
col[tab[,"SplitWeight"]==1e-6] <- "red"

pch <- rep(19, nrow(tab))
pch[tab$Tree] <- 17

par(mar=c(5.1, 4.1, 4.1, 8.1), xpd=TRUE)
plot(tab[,"SplitWeight"], tab[,"Bootstrap"], log="x", col=col, pch=pch,  
     xlab="Split weight (log scale)", ylab="Bootstrap (%)")
legend("topright", inset=c(-0.35,0), c("Pattern 1", "Pattern 2", "Pattern 3", 
                       "Pattern in the\nbest tree"), pch=c(19,19,19,17), 
                       col=c("blue", "green", "red", "black"), bty="n")

## -----------------------------------------------------------------------------
YCh <- read.tree(file.path(fdir, "RAxML_bestTree.YCh")) 
mtG <- read.tree(file.path(fdir, "RAxML_bestTree.mtG")) 
ncAI <- read.tree(file.path(fdir, "RAxML_bestTree.AIs")) 
all_data <- read.tree(file.path(fdir, "RAxML_bestTree.3moles")) 
YCh_boot <- read.tree(file.path(fdir, "RAxML_bootstrap.YCh")) 
mtG_boot <- read.tree(file.path(fdir, "RAxML_bootstrap.mtG")) 
ncAI_boot <- read.tree(file.path(fdir, "RAxML_bootstrap.AIs")) 
all_data_boot <- read.tree(file.path(fdir, "RAxML_bootstrap.3moles")) 

## -----------------------------------------------------------------------------
par(mfrow=c(2,2), mar = c(2,2,4,2))
YCh <- plotBS(midpoint(YCh), YCh_boot, "phylogram", p=0, main = "YCh")
mtG <- plotBS(midpoint(mtG), mtG_boot, "phylogram", p=0, main = "mtG")
ncAI <- plotBS(midpoint(ncAI), ncAI_boot, "phylogram", p=0, main = "ncAI")
all_data <- plotBS(midpoint(all_data), all_data_boot, "phylogram", p=0, 
                   main = "All data")

## ----small.mar=TRUE-----------------------------------------------------------
par(mfrow=c(1,1))
cn <- consensusNet(c(YCh, mtG, ncAI))
cn <- addConfidences(cn, YCh_boot) |> addConfidences(mtG_boot, add=TRUE) |> 
  addConfidences(ncAI_boot, add=TRUE) |> addConfidences(all_data_boot, add=TRUE)
plot(cn, "2D", show.edge.label=TRUE)