plot_trees.rmd

---
title: "Plot phylogenetic trees for target sequence capture data"
---

This document describes the process of plotting phylogenetic trees for sequence capture output
The trees need to be in Newick format and should be in the working directory

Results from IQ-TREE and ASTRAL will have multiple node labels (e.g. concordance factors)
To get the concordance trees into Newick properly, use the script `concord_to_newick.py`
Similarly, use `astral_parse.py` to separate node labels for ASTRAL (posterior probabilities and quartet support)

The outgroup file (`outgroup.txt`) should have a single tip sample ID per line for the outgroups
The samples file (`tree_samples.txt`) should have sample IDs and display labels tab-separated, one line per sample
```{r}
suppressMessages(library(ape))
suppressMessages(library(phytools))
concat_tree <- "concat.treefile"
concat_gcf <- "concord_newick_gcf.tre"
concat_scf <- "concord_newick_scf.tre"
astral_p <- "astral_p.tre"
astral_q <- "astral_q.tre"
astral_t <- "astral_poly.tre"
astral_tc <- "astral_poly_collapsed.tre"
outgroup_file <- "outgroup.txt"
samples_file <- "tree_samples.txt"
```

# Read in data
Read in the samples
```{r}
if (file.exists(samples_file)) {
	samples_present <- TRUE
	sample_table <- read.table(samples_file, sep = "\t")
	colnames(sample_table) <- c("ID", "label")
} else {
	samples_present <- FALSE
}
```

Read in the outgroups
```{r}
if (file.exists(outgroup_file)) {
	outgroup_present <- TRUE
	outgroups <- read.table(outgroup_file)[, 1]
} else {
	outgroup_present <- FALSE
}
```

Read in the trees
```{r}
tree_list <- vector("list", 7)
index <- 1

if (file.exists(concat_tree)) {
	tree_list[[index]] <- read.tree(concat_tree)
	concat_present <- TRUE
	concat_index <- index
	index <- index + 1
} else {
	concat_present <- FALSE
}

if (file.exists(concat_gcf)) {
	tree_list[[index]] <- read.tree(concat_gcf)
	concat_gcf_present <- TRUE
	concat_gcf_index <- index
	index <- index + 1
} else {
	concat_gcf_present <- FALSE
}

if (file.exists(concat_scf)) {
	tree_list[[index]] <- read.tree(concat_scf)
	concat_scf_present <- TRUE
	concat_scf_index <- index
	index <- index + 1
} else {
	concat_scf_present <- FALSE
}

if (file.exists(astral_p)) {
	tree_list[[index]] <- read.tree(astral_p)
	astral_p_present <- TRUE
	astral_p_index <- index
	index <- index + 1
} else {
	astral_p_present <- FALSE
}

if (file.exists(astral_q)) {
	tree_list[[index]] <- read.tree(astral_q)
	astral_q_present <- TRUE
	astral_q_index <- index
	index <- index + 1
} else {
	astral_q_present <- FALSE
}

if (file.exists(astral_t)) {
	tree_list[[index]] <- read.tree(astral_t)
	astral_t_present <- TRUE
	astral_t_index <- index
	index <- index + 1
} else {
	astral_t_present <- FALSE
}

if (file.exists(astral_tc)) {
	tree_list[[index]] <- read.tree(astral_tc)
	astral_tc_present <- TRUE
	astral_tc_index <- index
	index <- index + 1
} else {
   astral_tc_present <- FALSE
}
```

If there is an outgroup, root the trees
```{r}
if (outgroup_present) {
	for (index in seq_len(length(tree_list))) {
		tree <- tree_list[[index]]
		if (sum(outgroups %in% tree$tip.label) > 0) {
			these_outgroups <- outgroups[outgroups %in% tree$tip.label]
			if (length(these_outgroups) > 1) {		# root is set at common ancestor of outgroup
				if (is.monophyletic(tree, as.character(these_outgroups))) {
					rootnode <- getMRCA(tree, as.character(these_outgroups))
					position <- 0.5 * tree$edge.length[which(tree$edge[, 2] == rootnode)]
					rooted_tree <- reroot(tree, rootnode, position, edgelabel = TRUE)
					tree_list[[index]] <- rooted_tree
				} else {
					cat("Tree", index, "does not have monophyletic outgroup, so it is not rooted\n")
				}
			} else {	# single tip is outgroup
				tip_number <- which(tree$tip.label == these_outgroups)
				position <- 0.1 * tree$edge.length[which(tree$edge[, 2] == tip_number)]
				rooted_tree <- reroot(tree, node.number = tip_number, position, edgelabel = TRUE)
				tree_list[[index]] <- rooted_tree
			}
		} else {
			cat("Tree", index, "has no outgroups, so it is not rooted\n")
		}
	}
}
```

For later creation of figures, it may be desirable to get the tip order of ladderized trees (before labels are substituted)
This will output text files per tree, with sample IDs per line (order is from base of tree to top of page, i.e. starting from the outgroup)
Based on: https://stackoverflow.com/a/34364914
```{r}
# select names of trees and indices in the list (depending on what files are present)
tree_names <- vector("list")
indices <- vector("list")
if (concat_present) {
	tree_names <- append(tree_names, "concat")
	indices <- append(indices, concat_index)
}
if (astral_p_present) {
	tree_names <- append(tree_names, "astral")
	indices <- append(indices, astral_p_index)
}
if (astral_tc_present) {
	tree_names <- append(tree_names, "astral_poly")
	indices <- append(indices, astral_tc_index)
}
iter <- 1
for (index in indices) {
	lad_tree <- ladderize(tree_list[[index]], right = FALSE)
	# determine which edges are tips and get the order
	is_tip <- lad_tree$edge[, 2] <= length(lad_tree$tip.label)
	ordered_tips <- lad_tree$edge[is_tip, 2]
	# get the tips in order of plotting
	output_ordered_ids <- lad_tree$tip.label[ordered_tips]
	# write to a file
	connection <- file(paste0("tips_", tree_names[[iter]], ".txt"))
	writeLines(output_ordered_ids, connection)
	close(connection)
	iter <- iter + 1
}
```

# Grab support values
Capture concordance factors from the files present
```{r}
# gene concordance factors
if (concat_gcf_present) {
	tree <- tree_list[[concat_gcf_index]]
	mydf <- data.frame(matrix(ncol = 4, nrow = tree$Nnode))
	index <- 1
	for (label in tree$node.label) {
		if (any(label == "", label == "Root")) {
			values <- c("", "", "", "")
		} else {
			values <- strsplit(label, split = "\\/")[[1]]
		}
		mydf[index, 1:4] <- values
		index <- index + 1
	}
	gcf <- as.data.frame(lapply(mydf, as.numeric))
}

# site concordance factors
if (concat_scf_present) {
	tree <- tree_list[[concat_scf_index]]
	mydf <- data.frame(matrix(ncol = 3, nrow = tree$Nnode))
	index <- 1
	for (label in tree$node.label) {
		if (any(label == "", label == "Root")) {
			values <- c("", "", "")
		} else {
			values <- strsplit(label, split = "\\/")[[1]]
		}
		mydf[index, 1:3] <- values
		index <- index + 1
	}
	scf <- as.data.frame(lapply(mydf, as.numeric))
}
```

Capture ASTRAL posterior probabilities and quartet support if present
```{r}
# posterior probabilities
if (astral_p_present) {
	tree <- tree_list[[astral_p_index]]
	mydf <- data.frame(matrix(ncol = 1, nrow = tree$Nnode))
	index <- 1
	for (label in tree$node.label) {
		if (any(label == "", label == "Root")) {
			values <- c("")
		} else {
			values <- strsplit(label, split = "\\/")[[1]][1]
		}
		mydf[index, ] <- values
		index <- index + 1
	}
	mydf <- as.data.frame(lapply(mydf, as.numeric))
	pps <- sapply(mydf[, 1], function(x) format(round(x, 2), nsmall = 2))
	pps[pps == "NA"] <- ""
}

# quartet support
if (astral_q_present) {
	tree <- tree_list[[astral_q_index]]
	mydf <- data.frame(matrix(ncol = 3, nrow = tree$Nnode))
	index <- 1
	for (label in tree$node.label) {
		if (any(label == "", label == "Root")) {
			values <- c("", "", "")
		} else {
			values <- strsplit(label, split = "\\/")[[1]]
		}
		mydf[index, 1:3] <- values
		index <- index + 1
	}
	qs <- as.data.frame(lapply(mydf, as.numeric))
}
```

# Plot trees
Determine likely graphics parameters for the pages needed to display the trees
(The page should probably not be smaller than about 8 wide x 11 tall)
```{r}
num_tips <- length(tree_list[[1]]$tip.label)
# for Arial 12 point, each letter is ~ 1/6 inch; use roughly 1.5 times that (1/6 * 3/2 = 0.25)
my_height <- num_tips * 0.25
if (my_height < 11) {
	my_height <- 11
}
# determine width so that it is narrower with more terminals
my_width <- my_height * (0.8 - 0.001 * num_tips)
# roughly set a "cex" value for displaying pie charts based on number of tips (and therefore width)
my_cex <- 0.65 - 0.0025 * num_tips
```

It may be better to set these manually for your data and preference
```{r}
my_width <- 28
my_height <- 42
my_cex <- 0.25
```

If there are sample labels, substitute them
```{r}
if (samples_present) {
	for (index in seq_len(length(tree_list))) {
		tips <- tree_list[[index]]$tip.label
		new_tips <- tips
		for (tind in seq_len(length(tips))) {
			if (tips[tind] %in% sample_table$ID) {
				new_tips[tind] <- sample_table$label[sample_table$ID == tips[tind]]
			}
		}
		tree_list[[index]]$tip.label <- new_tips
	}
}
```

## Concatenation
Set the name to call the file and trees (run one of these two)
```{r}
if (concat_present) {
	tree_name <- "concatenation_gcf"
	tree_index <- concat_index
	if (concat_gcf_present) {
		concord <- gcf
		piecols <- c("white", "#c1c1c1", "#5b5a5a", "black")
	} else {
		cat("No concordance present, so don't attempt to plot pie charts below\n")
	}
} else {
	cat("The concatenation tree is needed to plot from this section\n")
}
```
*OR*
```{r}
if (concat_present) {
	tree_name <- "concatenation_scf"
	tree_index <- concat_index
	if (concat_scf_present) {
		concord <- scf
		piecols <- c("white", "grey", "black")
	} else {
		cat("No concordance present, so don't attempt to plot pie charts below\n")
	}
} else {
	cat("The concatenation tree is needed to plot from this section\n")
}
```

Start creating the graphic
```{r}
pdf(paste0(tree_name, ".pdf"), width = my_width, height = my_height, family = "ArialMT")
# set margins
par(mar = c(0.5, 0.5, 1, 0.5))
par(oma = c(0, 0, 3, 0))
```

Plot with edge lengths
```{r}
tree <- tree_list[[tree_index]]
plot.phylo(ladderize(tree, right = FALSE),
	no.margin = FALSE,
	font = 1,
	edge.width = 2,
	label.offset = max(nodeHeights(tree)) / 200,
	main = tree_name)
add.scale.bar(x = mean(par("usr")[1:2]),
	y = par("usr")[3] + 1,
	font = 1,
	lwd = 2)
```
*OR* plot in the reverse orientation
```{r}
tree <- tree_list[[tree_index]]
plot.phylo(ladderize(tree, right = FALSE),
	direction = "leftwards",
	no.margin = FALSE,
	font = 1,
	edge.width = 2,
	label.offset = max(nodeHeights(tree)) / 200,
	main = tree_name)
add.scale.bar(x = mean(par("usr")[1:2]),
	y = par("usr")[3] + 1,
	font = 1,
	lwd = 2)
```

Plot node labels (UFbootstrap) on edges
```{r}
drawSupportOnEdges(tree$node.label, adj = c(0.5, -0.5), frame = "none")
```
*OR* plot node labels on nodes
```{r}
nodelabels(tree$node.label, adj = c(-0.05, 0.5), frame = "none")
```

Plot concordance as pie charts on nodes
```{r}
nodelabels(pie = concord, piecol = piecols,	cex = my_cex)
```

Stop creating the graphic (can now run again choosing the other concordance if desired)
```{r}
invisible(dev.off())
```

## ASTRAL
Set the name to call the file and trees
```{r}
if (astral_p_present) {
	tree_name <- "astral"
	tree_index <- astral_p_index
	if (astral_q_present) {
		concord <- qs
		piecols <- c("white", "grey", "black")
	} else {
		cat("No quartet concordance present, so don't attempt to plot pie charts below\n")
	}
} else {
	cat("The ASTRAL tree is needed to plot from this section\n")
}
```

Start creating the graphic
```{r}
pdf(paste0(tree_name, ".pdf"), width = my_width, height = my_height, family = "ArialMT")
# set margins
par(mar = c(0.5, 0.5, 1, 0.5))
par(oma = c(0, 0, 3, 0))
```

Plot with edge lengths
```{r}
tree <- tree_list[[tree_index]]
plot.phylo(ladderize(tree, right = FALSE),
	no.margin = FALSE,
	font = 1,
	edge.width = 2,
	label.offset = max(nodeHeights(tree)) / 200,
	x.lim =  max(nodeHeights(tree)) * 1.3,
	main = tree_name)
add.scale.bar(x = mean(par("usr")[1:2]),
	y = par("usr")[3] + 1,
	font = 1,
	lwd = 2)
```
*OR* plot in the reverse orientation
```{r}
tree <- tree_list[[tree_index]]
plot.phylo(ladderize(tree, right = FALSE),
	direction = "leftwards",
	no.margin = FALSE,
	font = 1,
	edge.width = 2,
	label.offset = max(nodeHeights(tree)) / 200,
	x.lim =  max(nodeHeights(tree)) * 1.3,
	main = tree_name)
add.scale.bar(x = mean(par("usr")[1:2]),
	y = par("usr")[3] + 1,
	font = 1,
	lwd = 2)
```

Plot node labels (posterior prob) on edges
```{r}
drawSupportOnEdges(pps, adj = c(0.5, -0.5), frame = "none")
```
*OR* plot node labels on nodes
```{r}
nodelabels(pps, adj = c(-0.05, 0.5), frame = "none")
```

Plot quartet support as pie charts on nodes
```{r}
nodelabels(pie = concord, piecol = piecols, adj = 0.49, cex = my_cex)
```

Stop creating the graphic
```{r}
invisible(dev.off())
```

## ASTRAL_polytomy
Set the name to call the file and trees (run one of these two)
```{r}
if (astral_t_present) {
	tree_name <- "astral_poly"
	tree_index <- astral_t_index
} else {
	cat("The ASTRAL polytomy tree is needed to plot from this section\n")
}
```
*OR*
```{r}
if (astral_tc_present) {
	tree_name <- "astral_poly_collapse"
	tree_index <- astral_tc_index
} else {
	cat("The ASTRAL polytomy collapsed tree is needed to plot from this section\n")
}
```

Start creating the graphic
```{r}
pdf(paste0(tree_name, ".pdf"), width = my_width, height = my_height, family = "ArialMT")
# set margins
par(mar = c(0.5, 0.5, 1, 0.5))
par(oma = c(0, 0, 3, 0))
```

Plot
```{r}
tree <- tree_list[[tree_index]]
plot.phylo(ladderize(tree, right = FALSE),
	no.margin = FALSE,
	use.edge.length = FALSE,
	node.depth = 2,
	font = 1,
	edge.width = 2,
	label.offset = 0.2,
	main = tree_name)
```

Plot node labels (test p values) on edges
```{r}
drawSupportOnEdges(tree$node.label, adj = c(0.5, -0.5), frame = "none")
```

Stop creating the graphic
```{r}
invisible(dev.off())
```