Dendrogram.Rmd

# ============================================================
# AA vs Structural Tree Comparison (Tanglegram + Cluster Coloring)
#
# Description:
#   Compares a FoldMason structural guide tree (.nw) with an
#   IQ-TREE amino acid consensus tree (.contree) using a tanglegram.
#   Clusters are cut from the LEFT (structural) dendrogram and used
#   to color tip labels and connecting lines. Branches are neutral gray.
#
# Input:
#   - Structural guide tree:  <prefix>_structure.nw
#   - AA consensus tree:      <prefix>_AA.contree
#
# Output:
#   - <prefix>_AA_vs_Struct_tanglegram.png
#   - <prefix>_AA_vs_Struct_tree_metrics.csv
#
# Dependencies:
#   install.packages(c("ape", "phangorn", "dendextend", "viridisLite"))
#
# Usage:
#   Set `prefix` and `out_dir` below, then source the script.
# ============================================================

library(ape)
library(phangorn)
library(dendextend)
library(viridisLite)

# ----------------------------
# USER PARAMETERS
# ----------------------------
prefix  <- "lysB"                  # shared filename prefix
data_dir <- "data"                 # directory containing input trees
out_dir  <- "results"              # directory for outputs
h_rel    <- 0.25                   # relative height for cluster cutoff (0–1)
random_seed <- 1                   # for reproducible untangling

# ----------------------------
# #1 — Resolve file paths
# ----------------------------
tL_file <- file.path(data_dir, paste0(prefix, "_structure.nw"))
tR_file <- file.path(data_dir, paste0(prefix, "_AA.contree"))

stopifnot(file.exists(tL_file), file.exists(tR_file))

# ----------------------------
# #2 — Read trees
# ----------------------------
tL <- read.tree(tL_file)
tR <- read.tree(tR_file)

cat("Original tip counts:\n")
cat("  LEFT  (structural):", length(tL$tip.label), "\n")
cat("  RIGHT (AA):        ", length(tR$tip.label), "\n")

# ----------------------------
# #3 — Keep only shared tips
# ----------------------------
shared <- intersect(tL$tip.label, tR$tip.label)
tL <- keep.tip(tL, shared)
tR <- keep.tip(tR, shared)

cat("Shared tips:", length(shared), "\n")
stopifnot(setequal(tL$tip.label, tR$tip.label))

# Midpoint root both trees
tL <- midpoint(tL)
tR <- midpoint(tR)

# ----------------------------
# #4 — Resolve polytomies if needed
# ----------------------------
if (!is.binary(tL)) tL <- multi2di(tL)
if (!is.binary(tR)) tR <- multi2di(tR)

tL <- ladderize(tL)
tR <- ladderize(tR)

# ----------------------------
# #5 — Convert phylo → dendrogram via cophenetic distances
#       (avoids ultrametric requirement)
# ----------------------------
dL <- as.dendrogram(hclust(as.dist(cophenetic(tL)), method = "average"))
dR <- as.dendrogram(hclust(as.dist(cophenetic(tR)), method = "average"))

# ----------------------------
# #6 — Define clusters on LEFT tree at relative height
# ----------------------------
maxh <- attr(dL, "height")
clL  <- cutree(dL, h = h_rel * maxh)
k    <- length(unique(clL))
pal  <- setNames(viridis(k), sort(unique(clL)))

stopifnot(all(labels(dL) %in% names(clL)))
stopifnot(all(labels(dR) %in% names(clL)))

cat("Number of clusters (k =", k, ") at h_rel =", h_rel, "\n")

# ----------------------------
# #7 — Color tip labels by LEFT cluster membership
# ----------------------------
dL_col <- set(dL, "labels_col", pal[clL[labels(dL)]])
dR_col <- set(dR, "labels_col", pal[clL[labels(dR)]])

# Neutral branch styling
dL_col <- set(dL_col, "branches_col", "gray20")
dR_col <- set(dR_col, "branches_col", "gray20")
dL_col <- set(dL_col, "branches_lwd", 1)
dR_col <- set(dR_col, "branches_lwd", 1)

# ----------------------------
# #8 — Untangle to reduce line crossings
# ----------------------------
set.seed(random_seed)
unt    <- untangle(dL_col, dR_col, method = "step2side")
dL_col <- unt[[1]]
dR_col <- unt[[2]]

# Color connecting lines by LEFT cluster
line_cols <- pal[clL[labels(dL_col)]]

# ----------------------------
# #9 — Plot tanglegram
# ----------------------------
dir.create(out_dir, showWarnings = FALSE, recursive = TRUE)

png_file <- file.path(out_dir, paste0(prefix, "_AA_vs_Struct_tanglegram.png"))

png(filename = png_file, width = 12, height = 9, units = "in", res = 300)

tanglegram(
  dL_col, dR_col,
  lab.cex                    = 1,
  margin_inner               = 5.3,
  highlight_distinct_edges   = FALSE,
  common_subtrees_color_lines = FALSE,
  lwd                        = 1,
  color_lines                = line_cols
)

dev.off()
cat("Wrote:", png_file, "\n")

# ----------------------------
# #10 — Tree comparison metrics
# ----------------------------

# Normalized Robinson–Foulds distance (topology only)
rf  <- RF.dist(unroot(tL), unroot(tR), normalize = TRUE)

# Entanglement metric
ent <- entanglement(dL_col, dR_col)

cat("Normalized RF distance:", round(rf,  4), "\n")
cat("Entanglement:          ", round(ent, 4), "\n")

csv_file <- file.path(out_dir, paste0(prefix, "_AA_vs_Struct_tree_metrics.csv"))

write.csv(
  data.frame(
    metric = c("Normalized_RF", "Entanglement"),
    value  = c(rf, ent)
  ),
  csv_file,
  row.names = FALSE
)
cat("Wrote:", csv_file, "\n")