EGB - MC03: Introduction to Phylogenetic Analysis Using R

# Calculate the Hamming distance matrix from the aligned sequences
D_hamming <- dist.hamming(phyDat_msa_sample)

# Nearest Neighbour Clustering (NNC)
tree_NNC <- upgma(D_hamming, method = "single")
tree_NNC <- midpoint(tree_NNC)

# Furthest Neighbour (FN)
tree_FN <- upgma(D_hamming, method = "complete")
tree_FN <- midpoint(tree_FN)

# WPGMA (Weighted Pair Group Method with Arithmetic Mean)
tree_WPGMA <- wpgma(D_hamming)
tree_WPGMA <- midpoint(tree_WPGMA)

# UPGMC (Unweighted Pair Group Centroid Method)
tree_UPGMC <- upgma(D_hamming, method = "centroid")
tree_UPGMC <- midpoint(tree_UPGMC)

# WPGGMC (Weighted Pair Group Centroid Method)
tree_WPGGMC <- wpgma(D_hamming, method = "centroid")
tree_WPGGMC <- midpoint(tree_WPGGMC)

# UPGMA (Unweighted Pair Group Method with Arithmetic Mean)
tree_UPGMA <- upgma(D_hamming)
tree_UPGMA <- midpoint(tree_UPGMA)

# Ward's Minimum Variance Method
tree_WARD <- upgma(D_hamming, method = "ward.D2")
tree_WARD <- midpoint(tree_WARD)

# Define plotting function for ultrametric trees
plot_tree <- function(tree_plot, title_plot, max_x) {
  g <- ggtree(tree_plot, colour = "#00A499", size = 1) +
    geom_tiplab(size = 4, colour = "black", align = TRUE) +
    geom_nodepoint(size = 3, colour = "#c7254e") +
    labs(title = title_plot, size = 6) +
    xlim(0, max_x) +
    theme(
      axis.line = element_blank(),
      axis.text = element_blank(),
      panel.grid.major = element_blank(),
      panel.grid.minor = element_blank(),
      plot.margin = margin(t = 0, r = -5, b = 0, l = 0, unit = "pt"),
      legend.position = "top"
    ) +
    geom_text(aes(x = branch, label = round(branch, 2)), size = 3, colour = "black", vjust = -0.5, hjust = 0)
  return(g)
}

# Generate visualisations
g_NNC    <- plot_tree(tree_NNC,    "Nearest Neighbour Clustering (NNC)", 0.3)
g_FN     <- plot_tree(tree_FN,     "Furthest Neighbour (FN)", 0.3)
g_WPGMA  <- plot_tree(tree_WPGMA,  "WPGMA", 0.3)
g_UPGMC  <- plot_tree(tree_UPGMC,  "UPGMC", 0.3)
g_WPGGMC <- plot_tree(tree_WPGGMC, "WPGGMC", 0.3)
g_UPGMA  <- plot_tree(tree_UPGMA,  "UPGMA", 0.3)
g_WARD   <- plot_tree(tree_WARD,   "Ward's Minimum Variance", 0.3)

# Arrange all trees in a grid layout
grid.arrange(g_NNC, g_FN, g_WPGMA, g_UPGMC, g_WPGGMC, g_UPGMA, g_WARD, nrow = 4, ncol = 2)

suppressWarnings({
  suppressMessages({
    options(verbose = FALSE)
    
    # Distance Calculation
    # Compute the Hamming distance matrix for the given aligned sequences
    # This will serve as a measure of pairwise sequence dissimilarity
    D_hamming = dist.hamming(phyDat_msa_sample)

    # Base Tree Construction
    # Compute the Neighbor Joining (NJ) tree as a base tree
    # NJ is a distance-based method for constructing phylogenetic trees
    nj_tree = nj(D_hamming)

    # Correct Negative Branch Lengths
    # Negative branch lengths are not biologically meaningful,
    # hence set any negative branch lengths to 0
    nj_tree$edge.length[which(nj_tree$edge.length < 0)] = 0

    # Rooting the Tree
    # Root the NJ tree at its midpoint to ensure ultrametricity,
    # and convert any multifurcating nodes to bifurcating nodes
    nj_tree = midpoint(multi2di(nj_tree))

    # Parsimony optimisation
    # optimise the parsimony of the NJ tree
    # The optim.parsimony function refines the tree topology to minimize the parsimony score
    par_tree = optim.parsimony(nj_tree, phyDat_msa_sample)

    
    #Edge estimation
    par_tree=acctran(par_tree, phyDat_msa_sample)
    par_tree=midpoint(par_tree)
    
    # Calculating and Printing Parsimony
    # Calculate and print the parsimony score of the optimised tree
    # Lower parsimony scores indicate more parsimonious (simpler) trees
    print(paste("Parsimony Score for the optimised tree:", parsimony(par_tree, phyDat_msa_sample)))
  })
})

g_par_tree = plot_tree(par_tree, "Parsimony tree",2.8)
grid.arrange(g_par_tree, nrow = 1, ncol = 1)

# Suppress warnings and messages
suppressWarnings({
  suppressMessages({
    options(verbose = FALSE)
    
    # Distance Calculation
    # Compute the Hamming distance matrix for the given aligned sequences
    # This will serve as a measure of pairwise sequence dissimilarity
    D_hamming = dist.hamming(phyDat_msa_sample)

    # Base Tree Construction
    # Compute the Neighbor Joining (NJ) tree as a base tree
    # NJ is a distance-based method for constructing phylogenetic trees
    nj_tree = nj(D_hamming)

    # Correct Negative Branch Lengths
    # Negative branch lengths are not biologically meaningful,
    # hence set any negative branch lengths to 0
    nj_tree$edge.length[which(nj_tree$edge.length < 0)] = 0
    
    # Model Testing
    # Conduct a model test to determine the best-fitting evolutionary model based on AIC
    mt = modelTest(phyDat_msa_sample, nj_tree, model = c("JC", "K80", "HKY", "GTR"))
    var_mt = attr(mt, "env")
    evolutionary_model = eval(get(mt$Model[which.min(mt$AIC)], var_mt), var_mt)
    
    # Initialize PML Tree
    # Create an initial PML (Phylogenetic Maximum Likelihood) tree using the NJ tree and selected evolutionary model
    initial_pml_tree = pml(tree = nj_tree, data = phyDat_msa_sample, model = evolutionary_model$model,
                  bf = evolutionary_model$bf, Q = evolutionary_model$Q, inv = evolutionary_model$inv,
                  k = evolutionary_model$k, shape = evolutionary_model$shape)
    
    # optimise PML Tree
    # optimise the PML tree through various parameters and tree rearrangement methods
    pml_tree = optim.pml(initial_pml_tree, model = evolutionary_model$model, optInv = T, optGamma = F,
                         rearrangement = "stochastic", optBf = F, optQ = F,
                         optEdge = F, optNni = TRUE, control = pml.control(trace = 0))
  })
})

# Visualize the optimised PML Tree
# Create a visualization of the optimised PML tree with a specified title
pml_tree$tree = multi2di(pml_tree$tree)
g_pml_tree = plot_tree(pml_tree$tree, "Likelihood tree", 2.8)
grid.arrange(g_pml_tree, nrow = 1, ncol = 1)

# Print the selected evolutionary model with a descriptive message
cat("The parameters for this evolutionary model are:\n")

print(pml_tree)

library("kmer")

kmer_distance = round(kdistance(phyDat_msa_sample, k = 5),2)

# Print Kmers distance matrix as an HTML table with a 20% width
kable(as.matrix(kmer_distance), format = "html", table.attr = "style='width:20%;'")

# Compute the Neighbor Joining tree
nj_km_tree = nj(kmer_distance)
nj_km_tree$edge.length[which(nj_km_tree$edge.length<0)]=0
nj_km_tree = midpoint(multi2di(nj_km_tree))


# Visualize the optimised PML Tree
# Create a visualization of the optimised PML tree with a specified title
g_nj_km_tree = plot_tree(nj_km_tree, "NJ kmer tree", 2.8)
grid.arrange(g_nj_km_tree, nrow = 1, ncol = 1)

# ----------------------------- Bootstrapped NJ Tree with Node Support -----------------------------

# Calculate the Hamming distance matrix for the given phylogenetic alignment
distance <- dist.hamming(phyDat_msa_sample)

# Build a Neighbor Joining (NJ) tree based on the distance matrix
nj_tree <- nj(distance)

# Resolve multifurcations (convert multifurcating nodes into bifurcating)
nj_tree <- multi2di(nj_tree)

# Midpoint root the tree
nj_tree <- midpoint(nj_tree)

# Define a bootstrap function that builds an NJ tree from a resampled dataset
bootstrap_fun <- function(x) nj(dist.hamming(x))

# Perform bootstrap analysis using 100 trees
bootstrapped_trees <- boot.phylo(nj_tree, as.DNAbin(phyDat_msa_sample), bootstrap_fun, trees = TRUE)

# Add bootstrap support values to the internal nodes of the original NJ tree
bootstrapped_conf_tree <- addConfidences(nj_tree, bootstrapped_trees$trees)

# Replace NA support value at the root with 1 (if present)
na_index <- which(is.na(bootstrapped_conf_tree$node.label))
bootstrapped_conf_tree$node.label[na_index] <- 1

# ----------------------------- Tree Visualization -----------------------------

# Define a color palette for bootstrap values (gradient from red to green)
color_palette <- colorRampPalette(c("red3", "mediumseagreen"))(100)

# Map node support values (scaled 0–100) to the color palette
color_indices <- round(bootstrapped_conf_tree$node.label * 100)
node_colors <- color_palette[color_indices]

# Define tip label colors (highlight a specific tip, e.g., tip 111)
tip_colors <- rep("black", length(bootstrapped_conf_tree$tip.label))

# Plot the tree with bootstrap node support
g <- ggtree(bootstrapped_conf_tree, layout = "dendrogram", ladderize = TRUE)
g <- g + geom_tiplab(size = 4, color = tip_colors)  # Add tip labels with specified colors
g <- g + theme(plot.margin = margin(t = 0, r = 0, b = 150, l = 0, unit = "pt"), legend.position = "top")
g <- g + geom_nodepoint(size = 4, color = node_colors)  # Add colored nodes by support
g <- g + labs(title = "Representative consensus tree (medoid) with node support (100 trees)")
g <- g + theme(plot.title = element_text(size = 20))

# Create a data frame with values from 1 to 100 for color mapping
color_data = data.frame(x = 1:100, y = 1)
color_data$colors = color_palette

# Create a color bar plot using ggplot2
g1 = ggplot(color_data, aes(x = x, y = y, fill = colors)) +
  geom_bar(stat = "identity", width = 1, position = "identity") + scale_fill_manual(values = rev(color_data$colors)) +
  labs(title = NULL,
       x = "% of node support",
       y = NULL) +
  theme_minimal() + theme(axis.text.y = element_blank(),  # Remove y-axis labels
                          axis.title.y = element_blank(),
                          panel.grid = element_blank())  # Remove y-axis title
g1 = g1 + guides(fill = "none")

# Arrange the tree plot and the color bar plot
grid.arrange(g, g1, nrow = 2, ncol = 1, heights = c(0.8, 0.2))

# ----------------------------- Load Libraries -----------------------------
library(clValid)      
library(factoextra)
library(ggplot2)

# ----------------------------- KMeans Method Function -----------------------------
kmeans_method <- function(tree_distance_matrix) {
  possible_k <- seq(2, 3)
  set.seed(2)
  clustering_quality_results <- NULL
  for (i in possible_k) {
    clusters <- kmeans(tree_distance_matrix, i)
    silhouette_score <- mean(silhouette(clusters$cluster, tree_distance_matrix)[,3])
    tmp <- c("Kmeans", i, round(silhouette_score, 3))
    clustering_quality_results <- rbind(clustering_quality_results, tmp)
  }
  colnames(clustering_quality_results) <- c("Method", "k", "Silhouette")
  clustering_quality_results <- data.frame(clustering_quality_results)
  clustering_quality_results$Silhouette <- as.numeric(clustering_quality_results$Silhouette)
  return(clustering_quality_results)
}

# ----------------------------- RF Distance -----------------------------
trees <- c(nj_tree, bootstrapped_trees$trees)
names(trees) <- c("original", paste0("bootstrap_", seq_along(bootstrapped_trees$trees)))

tree_distance_bm <- RF.dist(trees, normalize = TRUE, check.labels = FALSE)
tree_distance_bm <- round(as.matrix(tree_distance_bm), 2)

# ----------------------------- Clustering -----------------------------
quality_kmeans <- kmeans_method(tree_distance_bm)
optimal_k <- as.numeric(quality_kmeans$k[which.max(quality_kmeans$Silhouette)])
set.seed(2)
kmeans_results <- kmeans(tree_distance_bm, optimal_k)

# ----------------------------- Visualización Mejorada -----------------------------
# Extraer coordenadas y asignar nombres
coordinates <- fviz_cluster(kmeans_results, data = tree_distance_bm)$data
coordinates$tree_name <- rownames(tree_distance_bm)
coordinates$label <- ifelse(coordinates$tree_name == "original", "Original NJ", "Bootstrap")

# Visualización con etiqueta al NJ original
g <- fviz_cluster(kmeans_results, data = tree_distance_bm, ellipse.type = "norm",
                  labelsize = 10, repel = TRUE, show.clust.cent = FALSE,
                  ggtheme = theme_classic(), pointsize = 3, shape = "circle",
                  palette = c("#00A499", "#695c59", "lightgreen")) +
     labs(title = "Phylogenetic Forest - 4 Sequences") +
     theme(axis.text = element_text(size = 12), legend.position = "bottom") +
     geom_text(data = coordinates[coordinates$label == "Original NJ", ],
               aes(x = x, y = y, label = label), color = "black", size = 4, fontface = "bold")

# Plot final
plot(g)