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

## ---- 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(2,2,2,2)) # 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
# NeighbourNet
plot.networx(Nnet,"2D")
## alternatively,
# plot(Nnet,"2D")


## ------------------------------------------------------------------------
# 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), mar=c(1,1,1,1))
#plotBS(raxml.tree, rotate.tree = 180) 
plot(raxml.tree, "u", rotate.tree = 180) 
edgelabels(raxml.tree$edge[,2],col="blue", frame="none")

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


## ------------------------------------------------------------------------
x <- addConfidences(Nnet,raxml.tree)
# 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
par(mfrow=c(1,1))
out.x <- plot.networx(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")

## ------------------------------------------------------------------------
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")

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

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

obj <- addConfidences(Nnet,cnet)
plot.networx(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")

