### R code from vignette source 'GrowingTrees.Rnw'

###################################################
### code chunk number 1: GrowingTrees.Rnw:51-53
###################################################
options(continue=" ")
options(width=80)


###################################################
### code chunk number 2: expr1
###################################################
library(DECIPHER)

# specify the path to your sequence file:
fas <- "<<path to FASTA file>>"
# OR find the example sequence file used in this tutorial:
fas <- system.file("extdata", "Streptomyces_ITS_aligned.fas", package="DECIPHER")

seqs <- readDNAStringSet(fas) # use readAAStringSet for amino acid sequences
seqs # the aligned sequences


###################################################
### code chunk number 3: expr2
###################################################
seqs <- unique(seqs) # remove duplicated sequences

ns <- gsub("^.*Streptomyces( subsp\\. | sp\\. | | sp_)([^ ]+).*$", "\\2", names(seqs))
names(seqs) <- ns # name by species (or any other preferred names)
seqs <- seqs[!duplicated(ns)] # remove redundant sequences from the same species
seqs


###################################################
### code chunk number 4: expr3
###################################################
getOption("SweaveHooks")[["fig"]]()
set.seed(123) # set the random number seed

tree <- Treeline(seqs,
	method="ML",
	model="GTR+G4",
	reconstruct=TRUE,
	maxTime=0.01)

set.seed(NULL) # reset seed

plot(tree)


###################################################
### code chunk number 5: expr4
###################################################
getOption("SweaveHooks")[["fig"]]()
new_tree <- MapCharacters(tree, labelEdges=TRUE)
plot(new_tree, edgePar=list(p.col=NA, p.border=NA, t.col="#55CC99", t.cex=0.7))
attr(new_tree[[1]], "change") # state changes on first branch left of (virtual) root


###################################################
### code chunk number 6: expr5
###################################################
#attributes(tree) # view all attributes
attr(tree, "members") # number of leaves below this (root) node
attr(tree, "height") # height of the node (in this case, the midpoint root)
attr(tree, "state") # ancestral state reconstruction (if reconstruct=TRUE)
head(attr(tree, "siteLnLs")) # LnL for every alignment column (site)
attr(tree, "score") # best score (in this case, the -LnL)
attr(tree, "model") # either the specified or automatically select transition model
attr(tree, "parameters") # the free model parameters (or NA if unoptimized)
attr(tree, "midpoint") # center of the edge (for plotting)


###################################################
### code chunk number 7: expr6
###################################################
getOption("SweaveHooks")[["fig"]]()
plot(dendrapply(tree,
	function(x) {
		s <- attr(x, "probability") # choose "probability" (aBayes) or "support"
		if (!is.null(s) && !is.na(s)) {
			s <- formatC(as.numeric(s), digits=2, format="f")
			attr(x, "edgetext") <- paste(s, "\n")
		}
		attr(x, "edgePar") <- list(p.col=NA, p.border=NA, t.col="#CC55AA", t.cex=0.7)
		if (is.leaf(x))
			attr(x, "nodePar") <- list(lab.font=3, pch=NA)
		x
	}),
	horiz=TRUE,
	yaxt='n')
# add a scale bar (placed manually)
arrows(0, 0, 0.4, 0, code=3, angle=90, len=0.05, xpd=TRUE)
text(0.2, 0, "0.4 subs./site", pos=3, xpd=TRUE)


###################################################
### code chunk number 8: expr7
###################################################
getOption("SweaveHooks")[["fig"]]()
reps <- 100 # number of bootstrap replicates

tree1 <- Treeline(seqs, verbose=FALSE, processors=1L)

partitions <- function(x) {
	if (is.leaf(x))
		return(NULL)
	x0 <- paste(sort(unlist(x)), collapse=" ")
	x1 <- partitions(x[[1]])
	x2 <- partitions(x[[2]])
	return(list(x0, x1, x2))
}

pBar <- txtProgressBar()
bootstraps <- vector("list", reps)
for (i in seq_len(reps)) {
	r <- sample(width(seqs)[1], replace=TRUE)
	at <- IRanges(r, width=1)
	seqs2 <- extractAt(seqs, at)
	seqs2 <- lapply(seqs2, unlist)
	seqs2 <- DNAStringSet(seqs2)
	
	temp <- Treeline(seqs2, verbose=FALSE)
	bootstraps[[i]] <- unlist(partitions(temp))
	setTxtProgressBar(pBar, i/reps)
}
close(pBar)

bootstraps <- table(unlist(bootstraps))
original <- unlist(partitions(tree1))
hits <- bootstraps[original]
names(hits) <- original
w <- which(is.na(hits))
if (length(w) > 0)
	hits[w] <- 0
hits <- round(hits/reps*100)

labelEdges <- function(x) {
	if (is.null(attributes(x)$leaf)) {
		part <- paste(sort(unlist(x)), collapse=" ")
		attr(x, "edgetext") <- as.character(hits[part])
	}
	return(x)
}
tree2 <- dendrapply(tree1, labelEdges)
attr(tree2, "edgetext") <- NULL

plot(tree2, edgePar=list(t.cex=0.5), nodePar=list(lab.cex=0.7, pch=NA))


###################################################
### code chunk number 9: expr8
###################################################
WriteDendrogram(tree, file="")


###################################################
### code chunk number 10: expr9
###################################################
params <- attr(tree, "parameters")
cat("[", paste(names(params), params, sep="=", collapse=","), "]", sep="", append=TRUE, file="")


###################################################
### code chunk number 11: sessinfo
###################################################
toLatex(sessionInfo(), locale=FALSE)