---
title: "Supplementary Code 1"
subtitle: "Similarity-based homolog classification algorithm"
format:
  html:
    toc: true
    number-sections: true
    code-fold: false
execute:
  eval: false
editor: source
  markdown: 
    wrap: 72
---

## Purpose

This script implements the similarity-based procedure used to assign coherent alpha and beta labels to homologous chromosome territories across a population of nuclei. It calculates pairwise matching energies, constructs one greedy minimum-energy path from each possible starting nucleus, propagates homolog labels along each path, and returns the candidate paths with their cumulative energies.

## Required input

The input table must contain the following columns:

| Column                  | Description                                                     |
|-------------------------|-----------------------------------------------------------------|
| `IdentificativoCellula` | Unique nucleus identifier                                       |
| `Object_Name`           | Chromosome territory identifier, for example `Chr06` or `Chr12` |
| `Pattern`               | Experimental-condition identifier                               |
| `Omologo`               | Provisional homolog label (`A` or `B`)                          |
| `Replica`               | Replicate identifier                                            |
| `X_circ`, `Y_circ`      | Normalized in-plane chromosome-territory coordinates            |

## Output

The function `run_homolog_energy_analysis()` returns a list containing:

-   `paths`: all chromosome-specific candidate paths, grouped by their starting nucleus;
-   `data`: a combined data frame containing the path index, frame number, propagated alpha/beta labels, pairwise energies, and cumulative energies.

## R code

```{r}
# ==============================================================================
# Supplementary Code 1
# Similarity-based homolog classification algorithm
#
# Purpose
# -------
# This script implements the similarity-based procedure used to assign coherent
# alpha and beta labels to homologous chromosome territories across a population
# of nuclei. The algorithm:
#   1. Calculates pairwise homolog-matching energies between nuclei.
#   2. Builds a greedy, minimum-energy path from each possible starting nucleus.
#   3. Propagates alpha/beta labels along every path.
#   4. Returns the candidate paths and their cumulative energies for downstream
#      identification of the globally lowest-energy path.
#
#
# Required input columns
# ----------------------
# IdentificativoCellula : unique nucleus identifier
# Object_Name           : chromosome territory identifier (e.g., "Chr06", "Chr12")
# Pattern               : experimental condition identifier
# Omologo               : provisional homolog label ("A" or "B")
# Replica               : replicate identifier
# X_circ, Y_circ        : normalized in-plane chromosome-territory coordinates
# ==============================================================================

suppressPackageStartupMessages({
  library(dplyr)
  library(tidyr)
  library(purrr)
})

# ------------------------------------------------------------------------------
# 1. Pairwise distances and energies
# ------------------------------------------------------------------------------

# Calculate the Euclidean distance between two homolog coordinates selected from
# two different nuclei.
calculate_distance <- function(data, cell_id_1, cell_id_2, homolog_1, homolog_2) {
  coordinates_1 <- data[
    data$IdentificativoCellula == cell_id_1,
    c(paste0("X_circ_", homolog_1), paste0("Y_circ_", homolog_1))
  ]
  coordinates_2 <- data[
    data$IdentificativoCellula == cell_id_2,
    c(paste0("X_circ_", homolog_2), paste0("Y_circ_", homolog_2))
  ]

  # Use shared coordinate names before calculating the Euclidean distance.
  names(coordinates_1) <- c("X", "Y")
  names(coordinates_2) <- c("X", "Y")

  distance <- as.numeric(dist(rbind(coordinates_1, coordinates_2)))
  return(distance)
}

# Convert a pairwise distance into the elastic energy used by the path-building
# procedure. The default stiffness is retained from the original implementation.
calculate_elastic_energy <- function(distance, stiffness = 1) {
  return(stiffness * distance^2)
}

# For a reference nucleus, calculate all four possible pairwise distances and
# energies to every other nucleus: AA, AB, BB, and BA.
calculate_system_energy <- function(data, reference_cell_id, stiffness = 1) {

  results <- purrr::map(1:nrow(data), function(row_index) {
    cell_id <- data$IdentificativoCellula[row_index]

    if (cell_id != reference_cell_id) {

      distance_AA <- calculate_distance(
        data, reference_cell_id, cell_id, "A", "A"
      )
      energy_AA <- calculate_elastic_energy(distance_AA, stiffness)

      distance_BB <- calculate_distance(
        data, reference_cell_id, cell_id, "B", "B"
      )
      energy_BB <- calculate_elastic_energy(distance_BB, stiffness)

      distance_AB <- calculate_distance(
        data, reference_cell_id, cell_id, "A", "B"
      )
      energy_AB <- calculate_elastic_energy(distance_AB, stiffness)

      distance_BA <- calculate_distance(
        data, reference_cell_id, cell_id, "B", "A"
      )
      energy_BA <- calculate_elastic_energy(distance_BA, stiffness)

      return(data.frame(
        ReferenceCellID = reference_cell_id,
        CellID = cell_id,

        Distance_AA = distance_AA,
        Energy_AA = energy_AA,

        Distance_AB = distance_AB,
        Energy_AB = energy_AB,

        Distance_BB = distance_BB,
        Energy_BB = energy_BB,

        Distance_BA = distance_BA,
        Energy_BA = energy_BA
      ))
    }
  })

  results <- do.call(rbind, results)
  return(results)
}

# Compare the two possible homolog pairings:
#   direct pairing:  A -> A and B -> B
#   swapped pairing: A -> B and B -> A
#
# The original strict comparison is retained. Therefore, exact ties are assigned
# to the swapped pairing (AB/BA), as in the source implementation.
select_homolog_pairing <- function(energy_data) {

  pairing_data <- energy_data %>%
    mutate(
      Pairing_A = if_else(
        (Energy_AA + Energy_BB) < (Energy_AB + Energy_BA),
        "Energy_AA",
        "Energy_AB"
      ),
      SelectedEnergy_A = if_else(
        Pairing_A == "Energy_AA",
        Energy_AA,
        Energy_AB
      ),
      Pairing_B = if_else(
        Pairing_A == "Energy_AA",
        "Energy_BB",
        "Energy_BA"
      ),
      SelectedEnergy_B = if_else(
        Pairing_A == "Energy_AA",
        Energy_BB,
        Energy_BA
      ),
      A = substr(Pairing_A, 8, 8),
      A_match = substr(Pairing_A, 9, 9),
      B = substr(Pairing_B, 8, 8),
      B_match = substr(Pairing_B, 9, 9),
      SelectedEnergy_A = as.numeric(SelectedEnergy_A),
      SelectedEnergy_B = as.numeric(SelectedEnergy_B),
      Energy_AA_AB_BB_BA = SelectedEnergy_A + SelectedEnergy_B
    ) %>%
    select(
      ReferenceCellID,
      CellID,
      SelectedEnergy_A,
      SelectedEnergy_B,
      Energy_AA_AB_BB_BA,
      A,
      A_match,
      B,
      B_match
    ) %>%
    rename(
      Energy_AA_AB = SelectedEnergy_A,
      Energy_BB_BA = SelectedEnergy_B
    )

  return(pairing_data)
}

# ------------------------------------------------------------------------------
# 2. Construction of one minimum-energy path
# ------------------------------------------------------------------------------

# Build a greedy minimum-energy path beginning at one selected nucleus.
# Each path is constructed independently for every candidate starting nucleus.
build_minimum_energy_path <- function(data, start_cell_id, chromosomes) {

  # Reformat chromosome-specific coordinates as one row per nucleus.
  wide_data <- data %>%
    pivot_wider(
      id_cols = c("IdentificativoCellula", "Pattern", "Replica"),
      names_from = Object_Name,
      values_from = c(X_circ_A, Y_circ_A, X_circ_B, Y_circ_B),
      names_sep = "_"
    )

  # Store chromosome-specific energy tables during each path iteration.
  energy_by_chromosome <- vector("list", length = length(chromosomes))
  names(energy_by_chromosome) <- chromosomes

  # Initialize the path with the user-selected starting nucleus.
  path_data <- data %>%
    filter(IdentificativoCellula == start_cell_id) %>%
    mutate(
      Energy_AA_AB_BB_BA = 0,
      Energy_AA_AB = 0,
      Energy_BB_BA = 0,
      A = "A",
      A_match = "A",
      B = "B",
      B_match = "B"
    ) %>%
    pivot_wider(
      id_cols = c("IdentificativoCellula", "Pattern", "Replica"),
      names_from = Object_Name,
      values_from = c(
        X_circ_A,
        Y_circ_A,
        X_circ_B,
        Y_circ_B,
        Energy_AA_AB_BB_BA,
        Energy_AA_AB,
        Energy_BB_BA,
        A,
        A_match,
        B,
        B_match
      ),
      names_sep = "_"
    ) %>%
    mutate(TotalEnergy = 0)

  # Track nuclei already included in the current path.
  used_cell_ids <- c(start_cell_id)

  # Select one new nucleus at each step until all nuclei are included.
  for (iteration in 1:(length(unique(data$IdentificativoCellula)) - 1)) {

    # Calculate candidate pairing energies for every chromosome.
    for (chromosome_index in 1:length(chromosomes)) {

      chromosome_data <- data %>%
        filter(
          Object_Name == chromosomes[chromosome_index] &
            !(
              IdentificativoCellula %in% used_cell_ids &
                !(IdentificativoCellula %in% start_cell_id)
            )
        )

      energy_data <- calculate_system_energy(
        chromosome_data,
        start_cell_id,
        stiffness = 0.5
      )

      pairing_data <- select_homolog_pairing(energy_data)

      energy_by_chromosome[[chromosomes[chromosome_index]]] <- pairing_data

      # Append the chromosome identifier to energy and homolog-pairing columns.
      colnames(energy_by_chromosome[[chromosomes[chromosome_index]]]) <- c(
        "ReferenceCellID",
        "CellID",
        paste("Energy_AA_AB", chromosomes[chromosome_index], sep = "_"),
        paste("Energy_BB_BA", chromosomes[chromosome_index], sep = "_"),
        paste("Energy_AA_AB_BB_BA", chromosomes[chromosome_index], sep = "_"),
        paste("A", chromosomes[chromosome_index], sep = "_"),
        paste("A_match", chromosomes[chromosome_index], sep = "_"),
        paste("B", chromosomes[chromosome_index], sep = "_"),
        paste("B_match", chromosomes[chromosome_index], sep = "_")
      )
    }

    # Join the chromosome-specific energy tables by reference and candidate IDs.
    merged_energy <- Reduce(
      function(x, y) merge(x, y, by = c("ReferenceCellID", "CellID")),
      energy_by_chromosome
    )

    # Sum energy across chromosomes for each candidate nucleus.
    energy_columns <- grep(
      "Energy_AA_AB_BB_BA",
      names(merged_energy),
      value = TRUE
    )

    merged_energy$TotalEnergy <- rowSums(
      merged_energy[energy_columns],
      na.rm = TRUE
    )

    # Exclude nuclei already included in the path.
    candidate_ids <- setdiff(merged_energy$CellID, used_cell_ids)

    # Select the candidate with the lowest total energy.
    minimum_energy_row <- merged_energy %>%
      filter(CellID %in% candidate_ids) %>%
      filter(TotalEnergy == min(TotalEnergy))

    # Append the selected nucleus and its pairing information to the path.
    new_frame <- wide_data %>%
      filter(IdentificativoCellula == minimum_energy_row$CellID)

    new_frame <- cbind(
      new_frame,
      minimum_energy_row[
        ,
        !colnames(minimum_energy_row) %in%
          c("ReferenceCellID", "CellID")
      ]
    )

    path_data <- rbind(path_data, new_frame)

    # Use the selected nucleus as the reference at the next iteration.
    start_cell_id <- as.character(minimum_energy_row$CellID)
    used_cell_ids <- c(used_cell_ids, start_cell_id)
  }

  return(path_data)
}

# ------------------------------------------------------------------------------
# 3. Propagation of alpha and beta labels along one path
# ------------------------------------------------------------------------------

# Resolve the provisional A/B labels into coherent alpha/beta labels along a
# single minimum-energy path.
classify_homologs_alpha_beta <- function(path_data, chromosomes) {

  # Separate the multi-chromosome path into chromosome-specific data frames.
  chromosome_paths <- lapply(chromosomes, function(chromosome) {
    path_data %>%
      select(-contains(chromosomes[chromosomes != chromosome])) %>%
      mutate(Object_Name = chromosome) %>%
      rename_with(
        ~ ifelse(
          grepl(chromosome, .),
          substr(., 1, nchar(.) - 6),
          .
        )
      )
  })

  names(chromosome_paths) <- chromosomes

  for (chromosome in chromosomes) {

    chromosome_path <- chromosome_paths[[chromosome]]

    # Initialize output columns.
    chromosome_path$Alpha <- NA
    chromosome_path$Beta <- NA

    chromosome_path$X_circ_Alpha <- NA
    chromosome_path$Y_circ_Alpha <- NA
    chromosome_path$X_circ_Beta <- NA
    chromosome_path$Y_circ_Beta <- NA

    chromosome_path$Distance_Alpha <- NA
    chromosome_path$Energy_Alpha <- NA
    chromosome_path$Distance_Beta <- NA
    chromosome_path$Energy_Beta <- NA

    # The first path frame provides the arbitrary alpha/beta seed labels.
    chromosome_path$Alpha[1] <- "A"
    chromosome_path$Beta[1] <- "B"

    chromosome_path$X_circ_Alpha[1] <- chromosome_path$X_circ_A
    chromosome_path$Y_circ_Alpha[1] <- chromosome_path$Y_circ_A
    chromosome_path$X_circ_Beta[1] <- chromosome_path$X_circ_B
    chromosome_path$Y_circ_Beta[1] <- chromosome_path$Y_circ_B

    chromosome_path$Distance_Alpha[1] <- 0
    chromosome_path$Energy_Alpha[1] <- 0
    chromosome_path$Distance_Beta[1] <- 0
    chromosome_path$Energy_Beta[1] <- 0

    # Retained for traceability with the original implementation. These flags
    # are not used by the subsequent calculations.
    alpha_switched <- FALSE
    beta_switched <- FALSE

    # Propagate homolog labels through consecutive path frames.
    for (row_index in 2:(nrow(chromosome_path))) {

      chromosome_path$Alpha[row_index] <- as.character(
        chromosome_path[
          row_index,
          paste(chromosome_path$Alpha[row_index - 1], "_match", sep = "")
        ]
      )

      chromosome_path$Beta[row_index] <- as.character(
        chromosome_path[
          row_index,
          paste(chromosome_path$Beta[row_index - 1], "_match", sep = "")
        ]
      )

      chromosome_path$X_circ_Alpha[row_index] <- as.numeric(
        chromosome_path[
          row_index,
          paste("X_circ_", chromosome_path$Alpha[row_index], sep = "")
        ]
      )

      chromosome_path$Y_circ_Alpha[row_index] <- as.numeric(
        chromosome_path[
          row_index,
          paste("Y_circ_", chromosome_path$Alpha[row_index], sep = "")
        ]
      )

      chromosome_path$X_circ_Beta[row_index] <- as.numeric(
        chromosome_path[
          row_index,
          paste("X_circ_", chromosome_path$Beta[row_index], sep = "")
        ]
      )

      chromosome_path$Y_circ_Beta[row_index] <- as.numeric(
        chromosome_path[
          row_index,
          paste("Y_circ_", chromosome_path$Beta[row_index], sep = "")
        ]
      )

      distance_alpha <- calculate_distance(
        chromosome_path,
        chromosome_path$IdentificativoCellula[row_index - 1],
        chromosome_path$IdentificativoCellula[row_index],
        "Alpha",
        "Alpha"
      )
      chromosome_path$Distance_Alpha[row_index] <- distance_alpha

      energy_alpha <- calculate_elastic_energy(
        distance_alpha,
        stiffness = 0.5
      )
      chromosome_path$Energy_Alpha[row_index] <- energy_alpha

      distance_beta <- calculate_distance(
        chromosome_path,
        chromosome_path$IdentificativoCellula[row_index - 1],
        chromosome_path$IdentificativoCellula[row_index],
        "Beta",
        "Beta"
      )
      chromosome_path$Distance_Beta[row_index] <- distance_beta

      energy_beta <- calculate_elastic_energy(
        distance_beta,
        stiffness = 0.5
      )
      chromosome_path$Energy_Beta[row_index] <- energy_beta
    }

    chromosome_path$Energy_AlphaBeta <-
      chromosome_path$Energy_Alpha + chromosome_path$Energy_Beta

    chromosome_paths[[chromosome]] <- chromosome_path
  }

  return(chromosome_paths)
}

# ------------------------------------------------------------------------------
# 4. Construction of all candidate paths
# ------------------------------------------------------------------------------

# Build and classify one path for each possible starting nucleus.
calculate_all_minimum_energy_paths <- function(
    data,
    cell_id_list,
    chromosomes
) {

  all_paths <- list()
  counter <- 1

  for (cell_id in cell_id_list) {

    start_cell_id <- as.character(cell_id)

    path_data <- build_minimum_energy_path(
      data,
      start_cell_id,
      chromosomes
    )

    classified_path <- classify_homologs_alpha_beta(
      path_data,
      chromosomes
    )

    # Add frame number and cumulative energy summaries.
    classified_path <- lapply(
      purrr::imap(chromosomes, ~ classified_path[[.x]]),
      function(path_table) {
        path_table %>%
          mutate(
            Frame = row_number(),
            CumulativeEnergy_Alpha = cumsum(Energy_Alpha),
            CumulativeEnergy_Beta = cumsum(Energy_Beta),
            CumulativeEnergy_AlphaBeta = cumsum(Energy_AlphaBeta),
            CumulativeTotalEnergy = cumsum(TotalEnergy)
          )
      }
    )

    names(classified_path) <- chromosomes
    all_paths[[cell_id]] <- classified_path

    print(paste(counter, "/", length(cell_id_list), sep = ""))
    counter <- counter + 1
  }

  return(all_paths)
}

# ------------------------------------------------------------------------------
# 5. Analysis wrapper
# ------------------------------------------------------------------------------

# Prepare input data, generate all candidate paths, and combine the resulting
# chromosome-specific tables into one output data frame.
run_homolog_energy_analysis <- function(
    dataset,
    chromosomes,
    pattern,
    cell_id_list
) {

  data <- dataset %>%
    select(
      IdentificativoCellula,
      Object_Name,
      Pattern,
      Omologo,
      Replica,
      X_circ,
      Y_circ
    ) %>%
    filter(
      Object_Name %in% chromosomes,
      Pattern == pattern,
      Replica == "A" | Replica == "B" | Replica == "C"
    ) %>%
    pivot_wider(
      id_cols = c("IdentificativoCellula", "Object_Name", "Pattern", "Replica"),
      names_from = Omologo,
      values_from = c(X_circ, Y_circ),
      names_sep = "_"
    ) %>%
    select(
      c("IdentificativoCellula", "Object_Name", "Pattern", "Replica"),
      starts_with("X_circ_A"),
      starts_with("Y_circ_A"),
      starts_with("X_circ_B"),
      starts_with("Y_circ_B")
    )

  all_paths <- calculate_all_minimum_energy_paths(
    data,
    cell_id_list,
    chromosomes
  )

  path_names <- names(all_paths)

  # Add the numerical path index to each chromosome-specific data frame.
  all_paths <- lapply(seq_along(all_paths), function(path_index) {

    all_paths[[path_index]] <- lapply(
      all_paths[[path_index]],
      function(path_table) {
        path_table %>%
          mutate(Path = path_index)
      }
    )

    return(all_paths[[path_index]])
  })

  names(all_paths) <- path_names

  # Combine chromosome-specific tables within, and then across, all paths.
  combined_output <- lapply(seq_along(all_paths), function(path_index) {
    as.data.frame(do.call(rbind, all_paths[[path_index]]))
  })

  combined_output <- as.data.frame(do.call(rbind, combined_output))

  results <- list(
    paths = all_paths,
    data = combined_output
  )

  beepr::beep(8)

  return(results)
}

# ==============================================================================
# 6. Interactive execution example
# ==============================================================================
# This section is retained as an interactive workflow. Before running it, replace
# the selected data file, output directory, pattern, chromosomes, and output file
# name as appropriate for the analysis.

# Select the Imaris-exported dataset.
data_file_path <- tcltk::tk_choose.files()
chromosome_dataset <- rio::import(data_file_path)
remove(data_file_path)

# Select the directory for saving the analysis output.
output_directory <- tcltk::tk_choose.dir()

# Define the candidate starting nuclei for the selected experimental condition.
cell_id_list <- chromosome_dataset %>%
  filter(Pattern == "F0.M0.0") %>%
  pull(IdentificativoCellula) %>%
  unique() %>%
  as.character()

# Calculate candidate minimum-energy paths for chromosomes 6 and 12.
analysis_results <- run_homolog_energy_analysis(
  dataset = chromosome_dataset,
  chromosomes = c("Chr06", "Chr12"),
  pattern = "F0.M0.0",
  cell_id_list = cell_id_list
)

# Save the output. Update the file name to reflect the analysed condition.
save(
  analysis_results,
  file = paste(
    output_directory,
    "/",
    "01_Demo_Top_EnergyClassification.RData",
    sep = "/"
  )
)

```
