Lesson 01 - Introduction to Tree Inference from Molecular Sequences Using R

# Distance Calculation
# Calculate the Hamming distance matrix for the given aligned sequences
# This serves as a measure of pairwise sequence dissimilarity for tree construction
D_hamming = dist.hamming(phyDat_msa_sample)

# Nearest Neighbor Clustering (NNC)
# Construct a phylogenetic tree using the Nearest Neighbor Clustering method
# This method groups sequences based on the nearest (smallest) pairwise distance
tree_NNC = upgma(D_hamming, "single")
tree_NNC = midpoint(tree_NNC)

# Furthest Neighbor (FN)
# Construct a phylogenetic tree using the Furthest Neighbor method
# This method groups sequences based on the furthest (largest) pairwise distance
tree_FN = upgma(D_hamming, "complete")
tree_FN = midpoint(tree_FN)

# Weighted Pair Group Method with Arithmetic Mean (WPGMA)
# Construct a phylogenetic tree using the WPGMA method
# This method considers all pairwise distances for clustering and calculates average distances
tree_WPGMA = wpgma(D_hamming)
tree_WPGMA = midpoint(tree_WPGMA)

# Unweighted Pair-Group Centroid Method (UPGMC)
# Construct a phylogenetic tree using the UPGMC method
# This method clusters sequences based on the centroid distance without considering the number of sequences in each cluster
tree_UPGMC = upgma(D_hamming, "centroid")
tree_UPGMC = midpoint(tree_UPGMC)

# Weighted Pair-Group Centroid Method (WPGGMC)
# Construct a phylogenetic tree using the WPGGMC method
# This method clusters sequences based on the centroid distance and considers the number of sequences in each cluster
tree_WPGGMC = wpgma(D_hamming, "centroid")
tree_WPGGMC = midpoint(tree_WPGGMC)

# Unweighted Pair Group Method with Arithmetic Mean (UPGMA)
# Construct a phylogenetic tree using the UPGMA method
# This method is similar to WPGMA but does not consider the number of sequences in each cluster for centroid calculation
tree_UPGMA = upgma(D_hamming)
tree_UPGMA = midpoint(tree_UPGMA)

# Ward's Minimum Variance Method
# Construct a phylogenetic tree using Ward's method
# This method minimizes the total within-cluster variance
tree_WARD = upgma(D_hamming, "ward")
tree_WARD = midpoint(tree_WARD)

# Function to plot trees
plot_tree = function(tree_plot, title_plot, max_x) {
  g = ggtree(tree_plot, color = "#00A499", size = 1)
  
  # Customize the appearance of the plot and tip labels
  g = g + geom_tiplab(size = 4, color = "black", align = TRUE) +
    geom_nodepoint(size = 3, color = "#c7254e") +
    labs(title = title_plot, size = 6) +
    xlim(0, max_x) +
    theme(
      # Remove axis lines and text
      axis.line = element_blank(),
      axis.text = element_blank(),
      # Remove all grids
      panel.grid.major = element_blank(),
      panel.grid.minor = element_blank(),
      # Adjust margins and legend position
      plot.margin = margin(t = 0, r = -5, b = 0, l = 0, unit = "pt"),
      legend.position = 'top'
    )
  
  # Add labels for branch lengths with a smaller font size
  g = g + geom_text(aes(x = branch, label = round(branch, 2)), size = 3, color = "black", vjust = -0.5, hjust = 0)
  
  return(g)
}

# Generate plots for each tree
g_NNC    = plot_tree(tree_NNC, "Nearest Neighbor Clustering (NNC)", 0.3)
g_FN     = plot_tree(tree_FN, "Furthest Neighbor (FN)", 0.3)
g_WPGMA  = plot_tree(tree_WPGMA, "Weighted Pair Group Method with Arithmetic Mean (WPGMA)", 0.3)
g_UPGMC  = plot_tree(tree_UPGMC, "Unweighted Pair-Group Centroid Method (UPGMC)", 0.3)
g_WPGGMC = plot_tree(tree_WPGGMC, "Weighted Pair-Group Centroid Method (WPGGMC)", 0.3)
g_UPGMA  = plot_tree(tree_UPGMA, "Unweighted Pair Group Method with Arithmetic Mean (UPGMA)", 0.3)
g_ward   = plot_tree(tree_WARD, "Ward's Minimum Variance Method", 0.3)

# Arrange the plots in a grid
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 Optimization
    # Optimize 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 optimized tree
    # Lower parsimony scores indicate more parsimonious (simpler) trees
    print(paste("Parsimony Score for the optimized 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)
    
    # Optimize PML Tree
    # Optimize 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 Optimized PML Tree
# Create a visualization of the optimized 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 Optimized PML Tree
# Create a visualization of the optimized 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)

# Calculate the Hamming distance for the given phylogenetic data
distance = dist.hamming(phyDat_msa_sample)

# Construct a Neighbor Joining tree based on the calculated distance
nj_tree = NJ(distance)

# Resolve multifurcations in the tree
nj_tree = multi2di(nj_tree)

# Place the root of the tree at its midpoint
nj_tree = midpoint(nj_tree)

# Define a function that computes the NJ tree based on Hamming distance
f = function(x) nj(dist.hamming(x))

# Perform bootstrap analysis on the tree
bootstrapped_trees = boot.phylo(nj_tree, as.DNAbin(phyDat_msa_sample), f, trees = T)

# Add bootstrap confidence values to the nodes of the original tree
bootstrapped_conf_tree = addConfidences(nj_tree, bootstrapped_trees$trees)

# Define a color palette for visualizing bootstrap support
color_palette = colorRampPalette(c("red3", "mediumseagreen"))(100)
color_positions = round(bootstrapped_conf_tree$node.label * 97, 0)
colors = color_palette[color_positions]

# Create a plot of the bootstrapped tree using ggtree
g = ggtree(bootstrapped_conf_tree, layout = "dendrogram", ladderize = T)
g = g + geom_tiplab(size = 4, color = "black")  # Add labels to terminal nodes
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 = colors) + theme(legend.position='top')
g = g + labs(title = "NJ tree Bootstrapping (100 trees)")

# 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 necessary libraries
library("clValid")          # For validation of clustering results
library("factoextra")       # For visualizing results

# Define a function for applying the KMEANS method
kmeans_method <- function(tree_distance_matrix) {
  
  # Analysis will be performed for a range of 2 to 3 clusters
  possible_k <- seq(2, 3)
  
  # Set seed for reproducibility
  set.seed(2)
  
  # Initialize a variable to store clustering quality indices
  clustering_quality_results <- NULL
  
  # Loop over possible k values to calculate Dunn, Connectivity and Silhouette indices
  for (i in possible_k) {
    # Perform clustering using KMEANS for different numbers of clusters
    clusters <- kmeans(tree_distance_matrix, i)
    
    # Calculate Silhouette score for the current number of clusters
    silhouette_score <- mean(silhouette(clusters$cluster, tree_distance_matrix)[,3])
    
    # Store results in a temporary variable
    tmp <- c("Kmeans", i, round(silhouette_score, 3))
    
    # Append results to the global variable clustering_quality_results
    clustering_quality_results <- rbind(clustering_quality_results, tmp)
  }
  
  # Set column names for the results data frame
  colnames(clustering_quality_results) <- c("Method", "k", "Silhouette")
  rownames(clustering_quality_results) <- NULL
  clustering_quality_results <- data.frame(clustering_quality_results)
  clustering_quality_results$Silhouette <- as.numeric(clustering_quality_results$Silhouette)
  
  # Return the clustering_quality_results data frame
  return(clustering_quality_results)
}

# Calculate Robinson-Foulds distance
trees <- c(nj_tree, bootstrapped_trees$trees)
tree_distance_bm <- RF.dist(trees, normalize = T, check.labels = F)
tree_distance_bm <- round(tree_distance_bm, 2)
tree_distance_bm <- as.matrix(tree_distance_bm)

# Apply KMEANS method and determine optimal number of clusters
quality_kmeans <- kmeans_method(tree_distance_bm)
optimal_k <- which.max(quality_kmeans$Silhouette)

# Perform kmeans clustering with optimal number of clusters
set.seed(2)
kmeans_results <- kmeans(tree_distance_bm, optimal_k)

# Visualization
coordinates <- fviz_cluster(kmeans_results, data = tree_distance_bm)$data
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"),)
g <- g + labs(title = "Phylogenetic Forest - 4 Sequences")
g <- g + theme(axis.text = element_text(size = 12), legend.position = "bottom")

# Plot the result
plot(g)