/usr/lib/R/site-library/phangorn/doc/IntertwiningTreesAndNetworks.R is in r-cran-phangorn 2.4.0-1.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 | ## ----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
library(magrittr)
## ------------------------------------------------------------------------
## 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
# NeighbourNet
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)
|