############################################################
# GSE121212 xCell analysis
#
# Input files:
#   GSE121212_input1_normalized_log2_filled_gene_counts.txt
#   GSE121212_input2_normalized_log2_low_signal_cutoff_filled_gene_counts.txt
#
# File format:
#   Row 1: status
#   Row 2: patient
#   Row 3: sample ID
#   Row 4 onward:
#       Column 1 = Entrez Gene ID
#       Column 2 onward = expression values
#
# Analyses:
#   1. Convert Entrez Gene IDs to HGNC Gene Symbols
#   2. Run xCell on all 40 samples together
#   3. Compare non-lesional and lesional samples without pairing
#   4. Compare non-lesional and lesional samples with pairing
#   5. Compare xCell results between Input 1 and Input 2
############################################################


############################
# 1. Package installation
############################

cran_packages <- c(
  "remotes",
  "dplyr",
  "tidyr",
  "readr",
  "tibble",
  "ggplot2"
)

for (pkg in cran_packages) {
  if (!requireNamespace(pkg, quietly = TRUE)) {
    install.packages(pkg)
  }
}

if (!requireNamespace("BiocManager", quietly = TRUE)) {
  install.packages("BiocManager")
}

bioconductor_packages <- c(
  "AnnotationDbi",
  "org.Hs.eg.db"
)

for (pkg in bioconductor_packages) {
  if (!requireNamespace(pkg, quietly = TRUE)) {
    BiocManager::install(pkg, ask = FALSE, update = FALSE)
  }
}

if (!requireNamespace("xCell", quietly = TRUE)) {
  remotes::install_github("dviraran/xCell")
}


############################
# 2. Load packages
############################

library(xCell)
library(dplyr)
library(tidyr)
library(readr)
library(tibble)
library(ggplot2)
library(AnnotationDbi)
library(org.Hs.eg.db)


############################
# 3. File settings
############################

input1_file <-
  "GSE121212_input1_normalized_log2_filled_gene_counts.txt"

input2_file <-
  paste0(
    "GSE121212_input2_normalized_log2_",
    "low_signal_cutoff_filled_gene_counts.txt"
  )

output_directory <- "GSE121212_xCell_results"

if (!dir.exists(output_directory)) {
  dir.create(output_directory, recursive = TRUE)
}


############################
# 4. Read the Subio-style
#    expression text file
############################

read_subio_expression_file <- function(file_path) {

  if (!file.exists(file_path)) {
    stop("Input file was not found: ", file_path)
  }

  message("Reading: ", file_path)

  # Read every row without treating the first row as a header.
  raw_data <- read.delim(
    file_path,
    header = FALSE,
    sep = "\t",
    check.names = FALSE,
    stringsAsFactors = FALSE,
    quote = "",
    comment.char = "",
    fill = TRUE
  )

  if (nrow(raw_data) < 4) {
    stop(
      "The input file must contain three metadata rows ",
      "followed by expression data."
    )
  }

  if (ncol(raw_data) < 3) {
    stop("The input file does not contain enough sample columns.")
  }

  # Extract the first three metadata rows.
  status_values  <- trimws(as.character(unlist(raw_data[1, -1])))
  patient_values <- trimws(as.character(unlist(raw_data[2, -1])))
  sample_ids     <- trimws(as.character(unlist(raw_data[3, -1])))

  metadata <- data.frame(
    sample_id = sample_ids,
    patient   = patient_values,
    status    = status_values,
    stringsAsFactors = FALSE
  )

  # Standardize status labels.
  metadata$status <- tolower(metadata$status)
  metadata$status <- gsub("_", "-", metadata$status)
  metadata$status <- gsub(" ", "-", metadata$status)

  expected_statuses <- c("non-lesional", "lesional")

  unexpected_statuses <- setdiff(
    unique(metadata$status),
    expected_statuses
  )

  if (length(unexpected_statuses) > 0) {
    stop(
      "Unexpected status label(s): ",
      paste(unexpected_statuses, collapse = ", ")
    )
  }

  if (anyDuplicated(metadata$sample_id)) {
    duplicated_ids <- unique(
      metadata$sample_id[duplicated(metadata$sample_id)]
    )

    stop(
      "Duplicated sample ID(s): ",
      paste(duplicated_ids, collapse = ", ")
    )
  }

  # Check that every patient has one sample from each status.
  pair_check <- metadata |>
    count(patient, status, name = "n") |>
    tidyr::pivot_wider(
      names_from = status,
      values_from = n,
      values_fill = 0
    )

  if (
    !all(c("non-lesional", "lesional") %in% colnames(pair_check))
  ) {
    stop("Both non-lesional and lesional samples are required.")
  }

  invalid_pairs <- pair_check |>
    filter(
      .data[["non-lesional"]] != 1 |
        .data[["lesional"]] != 1
    )

  if (nrow(invalid_pairs) > 0) {
    print(invalid_pairs)

    stop(
      "One or more patients do not have exactly one ",
      "non-lesional and one lesional sample."
    )
  }

  # Expression data start from the fourth row.
  expression_data <- raw_data[-c(1, 2, 3), , drop = FALSE]

  entrez_ids <- trimws(as.character(expression_data[[1]]))

  # Convert all expression columns to numeric.
  expression_values <- expression_data[, -1, drop = FALSE]

  expression_values[] <- lapply(
    expression_values,
    function(x) {
      suppressWarnings(as.numeric(as.character(x)))
    }
  )

  expression_matrix <- as.matrix(expression_values)
  storage.mode(expression_matrix) <- "numeric"

  colnames(expression_matrix) <- metadata$sample_id
  rownames(expression_matrix) <- entrez_ids

  # Remove blank Gene IDs.
  valid_gene_id <- !is.na(entrez_ids) & entrez_ids != ""

  expression_matrix <- expression_matrix[
    valid_gene_id,
    ,
    drop = FALSE
  ]

  entrez_ids <- entrez_ids[valid_gene_id]

  # Remove rows containing non-finite expression values.
  finite_rows <- apply(
    expression_matrix,
    1,
    function(x) all(is.finite(x))
  )

  if (any(!finite_rows)) {
    warning(
      sum(!finite_rows),
      " gene row(s) containing NA, NaN or infinite values were removed."
    )

    expression_matrix <- expression_matrix[
      finite_rows,
      ,
      drop = FALSE
    ]

    entrez_ids <- entrez_ids[finite_rows]
  }

  # Remove duplicated Entrez IDs by keeping the maximum value
  # for each sample.
  expression_df <- as.data.frame(
    expression_matrix,
    check.names = FALSE
  )

  expression_df$ENTREZID <- entrez_ids

  expression_df <- expression_df |>
    group_by(ENTREZID) |>
    summarise(
      across(
        all_of(metadata$sample_id),
        ~ max(.x, na.rm = TRUE)
      ),
      .groups = "drop"
    )

  list(
    expression_entrez = expression_df,
    metadata = metadata
  )
}


############################
# 5. Convert Entrez IDs to
#    HGNC Gene Symbols
############################

convert_entrez_to_symbol <- function(expression_entrez) {

  entrez_ids <- as.character(expression_entrez$ENTREZID)

  gene_symbols <- AnnotationDbi::mapIds(
    org.Hs.eg.db,
    keys = entrez_ids,
    column = "SYMBOL",
    keytype = "ENTREZID",
    multiVals = "first"
  )

  converted_data <- expression_entrez
  converted_data$SYMBOL <- unname(gene_symbols[entrez_ids])

  mapped <- (
    !is.na(converted_data$SYMBOL) &
      converted_data$SYMBOL != ""
  )

  message(
    "Entrez IDs supplied: ",
    nrow(converted_data)
  )

  message(
    "Entrez IDs mapped to Gene Symbols: ",
    sum(mapped)
  )

  message(
    "Unmapped Entrez IDs: ",
    sum(!mapped)
  )

  converted_data <- converted_data[mapped, , drop = FALSE]

  sample_columns <- setdiff(
    colnames(converted_data),
    c("ENTREZID", "SYMBOL")
  )

  # Multiple Entrez IDs mapped to the same Gene Symbol are combined.
  # The maximum expression value is retained for each sample.
  symbol_data <- converted_data |>
    dplyr::select(
      SYMBOL,
      dplyr::all_of(sample_columns)
    ) |>
    dplyr::group_by(SYMBOL) |>
    dplyr::summarise(
      dplyr::across(
        dplyr::all_of(sample_columns),
        ~ max(.x, na.rm = TRUE)
      ),
      .groups = "drop"
    )

  expression_matrix <- as.matrix(
    symbol_data[, sample_columns, drop = FALSE]
  )

  storage.mode(expression_matrix) <- "numeric"
  rownames(expression_matrix) <- symbol_data$SYMBOL

  if (anyDuplicated(rownames(expression_matrix))) {
    stop("Duplicated Gene Symbols remain after aggregation.")
  }

  message(
    "Unique Gene Symbols: ",
    nrow(expression_matrix)
  )

  expression_matrix
}


############################
# 6. Run xCell
############################

run_xcell <- function(
  file_path,
  input_name,
  output_directory
) {

  imported <- read_subio_expression_file(file_path)

  expression_matrix <- convert_entrez_to_symbol(
    imported$expression_entrez
  )

  metadata <- imported$metadata

  # Ensure that expression columns follow the metadata order.
  expression_matrix <- expression_matrix[
    ,
    metadata$sample_id,
    drop = FALSE
  ]

  message(
    "Running xCell for ",
    input_name,
    " with ",
    nrow(expression_matrix),
    " genes and ",
    ncol(expression_matrix),
    " samples."
  )

  # Run all 40 samples together.
  xcell_scores <- xCell::xCellAnalysis(
    expression_matrix,
    rnaseq = TRUE
  )

  xcell_scores <- as.matrix(xcell_scores)

  # Save the original cell type × sample matrix.
  write.table(
    data.frame(
      cell_type = rownames(xcell_scores),
      xcell_scores,
      check.names = FALSE
    ),
    file = file.path(
      output_directory,
      paste0(input_name, "_xCell_scores.txt")
    ),
    sep = "\t",
    quote = FALSE,
    row.names = FALSE
  )

  # Convert to sample × cell type format.
  score_by_sample <- as.data.frame(
    t(xcell_scores),
    check.names = FALSE
  )

  score_by_sample$sample_id <- rownames(score_by_sample)

  score_by_sample <- metadata |>
    left_join(
      score_by_sample,
      by = "sample_id"
    )

  write.table(
    score_by_sample,
    file = file.path(
      output_directory,
      paste0(input_name, "_xCell_scores_with_metadata.txt")
    ),
    sep = "\t",
    quote = FALSE,
    row.names = FALSE
  )

  list(
    expression_matrix = expression_matrix,
    metadata = metadata,
    xcell_scores = xcell_scores,
    score_by_sample = score_by_sample
  )
}


############################
# 7. Statistical comparison
############################

analyse_xcell_scores <- function(
  xcell_result,
  input_name,
  output_directory
) {

  score_data <- xcell_result$score_by_sample

  cell_types <- setdiff(
    colnames(score_data),
    c("sample_id", "patient", "status")
  )

  results <- lapply(
    cell_types,
    function(cell_type) {

      temporary <- score_data |>
        dplyr::select(
          sample_id,
          patient,
          status,
          score = all_of(cell_type)
        )

      non_lesional_scores <- temporary |>
        filter(status == "non-lesional") |>
        pull(score)

      lesional_scores <- temporary |>
        filter(status == "lesional") |>
        pull(score)

      # Unpaired Welch t-test.
      unpaired_test <- tryCatch(
        t.test(
          lesional_scores,
          non_lesional_scores,
          paired = FALSE,
          var.equal = FALSE
        ),
        error = function(e) NULL
      )

      # Arrange paired values by patient.
      paired_data <- temporary |>
        dplyr::select(patient, status, score) |>
        pivot_wider(
          names_from = status,
          values_from = score
        ) |>
        arrange(patient) |>
        mutate(
          paired_difference =
            .data[["lesional"]] -
            .data[["non-lesional"]]
        )

      paired_test <- tryCatch(
        t.test(
          paired_data$lesional,
          paired_data$`non-lesional`,
          paired = TRUE
        ),
        error = function(e) NULL
      )

      mean_non_lesional <- mean(
        paired_data$`non-lesional`,
        na.rm = TRUE
      )

      mean_lesional <- mean(
        paired_data$lesional,
        na.rm = TRUE
      )

      mean_paired_difference <- mean(
        paired_data$paired_difference,
        na.rm = TRUE
      )

      median_paired_difference <- median(
        paired_data$paired_difference,
        na.rm = TRUE
      )

      n_increased <- sum(
        paired_data$paired_difference > 0,
        na.rm = TRUE
      )

      n_decreased <- sum(
        paired_data$paired_difference < 0,
        na.rm = TRUE
      )

      n_unchanged <- sum(
        paired_data$paired_difference == 0,
        na.rm = TRUE
      )

      direction_consistency <- max(
        n_increased,
        n_decreased
      ) / nrow(paired_data)

      data.frame(
        cell_type = cell_type,
        n_patients = nrow(paired_data),
        mean_non_lesional = mean_non_lesional,
        mean_lesional = mean_lesional,
        unpaired_mean_difference =
          mean_lesional - mean_non_lesional,
        unpaired_p_value = if (is.null(unpaired_test)) {
          NA_real_
        } else {
          unpaired_test$p.value
        },
        mean_paired_difference = mean_paired_difference,
        median_paired_difference = median_paired_difference,
        paired_p_value = if (is.null(paired_test)) {
          NA_real_
        } else {
          paired_test$p.value
        },
        patients_increased = n_increased,
        patients_decreased = n_decreased,
        patients_unchanged = n_unchanged,
        direction_consistency = direction_consistency,
        stringsAsFactors = FALSE
      )
    }
  )

  results <- bind_rows(results)

  results$unpaired_adjusted_p_value <- p.adjust(
    results$unpaired_p_value,
    method = "BH"
  )

  results$paired_adjusted_p_value <- p.adjust(
    results$paired_p_value,
    method = "BH"
  )

  results <- results |>
    arrange(paired_p_value)

  write.table(
    results,
    file = file.path(
      output_directory,
      paste0(input_name, "_statistical_comparison.txt")
    ),
    sep = "\t",
    quote = FALSE,
    row.names = FALSE
  )

  # Save all patient-level paired differences.
  paired_differences <- score_data |>
    pivot_longer(
      cols = all_of(cell_types),
      names_to = "cell_type",
      values_to = "score"
    ) |>
    dplyr::select(patient, status, cell_type, score) |>
    pivot_wider(
      names_from = status,
      values_from = score
    ) |>
    mutate(
      paired_difference =
        .data[["lesional"]] -
        .data[["non-lesional"]]
    ) |>
    arrange(cell_type, patient)

  write.table(
    paired_differences,
    file = file.path(
      output_directory,
      paste0(input_name, "_patient_level_paired_differences.txt")
    ),
    sep = "\t",
    quote = FALSE,
    row.names = FALSE
  )

  list(
    statistics = results,
    paired_differences = paired_differences
  )
}


############################
# 8. Run both input files
############################

input1_result <- run_xcell(
  file_path = input1_file,
  input_name = "input1",
  output_directory = output_directory
)

input2_result <- run_xcell(
  file_path = input2_file,
  input_name = "input2",
  output_directory = output_directory
)

input1_analysis <- analyse_xcell_scores(
  xcell_result = input1_result,
  input_name = "input1",
  output_directory = output_directory
)

input2_analysis <- analyse_xcell_scores(
  xcell_result = input2_result,
  input_name = "input2",
  output_directory = output_directory
)


############################
# 9. Compare Input 1 and
#    Input 2 xCell scores
############################

compare_input_scores <- function(
  input1_result,
  input2_result,
  output_directory
) {

  scores1 <- input1_result$xcell_scores
  scores2 <- input2_result$xcell_scores

  common_cell_types <- intersect(
    rownames(scores1),
    rownames(scores2)
  )

  common_samples <- intersect(
    colnames(scores1),
    colnames(scores2)
  )

  scores1 <- scores1[
    common_cell_types,
    common_samples,
    drop = FALSE
  ]

  scores2 <- scores2[
    common_cell_types,
    common_samples,
    drop = FALSE
  ]

  comparison <- lapply(
    common_cell_types,
    function(cell_type) {

      values1 <- as.numeric(scores1[cell_type, ])
      values2 <- as.numeric(scores2[cell_type, ])

      pearson_correlation <- tryCatch(
        cor(
          values1,
          values2,
          method = "pearson",
          use = "complete.obs"
        ),
        error = function(e) NA_real_
      )

      spearman_correlation <- tryCatch(
        cor(
          values1,
          values2,
          method = "spearman",
          use = "complete.obs"
        ),
        error = function(e) NA_real_
      )

      data.frame(
        cell_type = cell_type,
        pearson_correlation = pearson_correlation,
        spearman_correlation = spearman_correlation,
        mean_input1 = mean(values1, na.rm = TRUE),
        mean_input2 = mean(values2, na.rm = TRUE),
        mean_difference_input2_minus_input1 =
          mean(values2 - values1, na.rm = TRUE),
        maximum_absolute_difference =
          max(abs(values2 - values1), na.rm = TRUE),
        stringsAsFactors = FALSE
      )
    }
  )

  comparison <- bind_rows(comparison) |>
    arrange(spearman_correlation)

  write.table(
    comparison,
    file = file.path(
      output_directory,
      "input1_vs_input2_xCell_score_comparison.txt"
    ),
    sep = "\t",
    quote = FALSE,
    row.names = FALSE
  )

  comparison
}

input_score_comparison <- compare_input_scores(
  input1_result = input1_result,
  input2_result = input2_result,
  output_directory = output_directory
)


############################
# 10. Compare statistical
#     conclusions
############################

statistical_comparison <- input1_analysis$statistics |>
  dplyr::select(
    cell_type,
    input1_unpaired_difference =
      unpaired_mean_difference,
    input1_unpaired_p =
      unpaired_p_value,
    input1_unpaired_adjusted_p =
      unpaired_adjusted_p_value,
    input1_paired_difference =
      mean_paired_difference,
    input1_paired_p =
      paired_p_value,
    input1_paired_adjusted_p =
      paired_adjusted_p_value,
    input1_patients_increased =
      patients_increased,
    input1_patients_decreased =
      patients_decreased
  ) |>
  full_join(
    input2_analysis$statistics |>
      dplyr::select(
        cell_type,
        input2_unpaired_difference =
          unpaired_mean_difference,
        input2_unpaired_p =
          unpaired_p_value,
        input2_unpaired_adjusted_p =
          unpaired_adjusted_p_value,
        input2_paired_difference =
          mean_paired_difference,
        input2_paired_p =
          paired_p_value,
        input2_paired_adjusted_p =
          paired_adjusted_p_value,
        input2_patients_increased =
          patients_increased,
        input2_patients_decreased =
          patients_decreased
      ),
    by = "cell_type"
  ) |>
  mutate(
    same_paired_direction =
      sign(input1_paired_difference) ==
      sign(input2_paired_difference),

    paired_significant_input1 =
      input1_paired_p < 0.05,

    paired_significant_input2 =
      input2_paired_p < 0.05,

    paired_significant_both =
      paired_significant_input1 &
      paired_significant_input2
  ) |>
  arrange(
    input2_paired_p,
    input1_paired_p
  )

write.table(
  statistical_comparison,
  file = file.path(
    output_directory,
    "input1_vs_input2_statistical_result_comparison.txt"
  ),
  sep = "\t",
  quote = FALSE,
  row.names = FALSE
)


############################
# 11. Create a paired plot
#     for selected cell types
############################

create_paired_plot <- function(
  xcell_result,
  cell_type,
  input_name,
  output_directory
) {

  score_data <- xcell_result$score_by_sample

  if (!cell_type %in% colnames(score_data)) {
    warning(
      "Cell type was not found: ",
      cell_type
    )

    return(invisible(NULL))
  }

  plot_data <- score_data |>
    dplyr::select(
      sample_id,
      patient,
      status,
      score = all_of(cell_type)
    ) |>
    mutate(
      status = factor(
        status,
        levels = c("non-lesional", "lesional")
      )
    )

  plot_object <- ggplot(
    plot_data,
    aes(
      x = status,
      y = score,
      group = patient
    )
  ) +
    geom_line(alpha = 0.55) +
    geom_point(size = 2) +
    labs(
      title = paste0(
        cell_type,
        " — ",
        input_name
      ),
      x = NULL,
      y = "xCell enrichment score"
    ) +
    theme_classic(base_size = 12)

  ggsave(
    filename = file.path(
      output_directory,
      paste0(
        input_name,
        "_",
        gsub("[^A-Za-z0-9]+", "_", cell_type),
        "_paired_plot.png"
      )
    ),
    plot = plot_object,
    width = 5,
    height = 5,
    dpi = 300
  )

  plot_object
}


############################
# 12. Example plots
############################

# Replace these names with cell types present in the xCell output.
example_cell_types <- c(
  "CD8+ T-cells",
  "Macrophages",
  "Fibroblasts"
)

for (cell_type in example_cell_types) {

  create_paired_plot(
    xcell_result = input1_result,
    cell_type = cell_type,
    input_name = "input1",
    output_directory = output_directory
  )

  create_paired_plot(
    xcell_result = input2_result,
    cell_type = cell_type,
    input_name = "input2",
    output_directory = output_directory
  )
}


############################
# 13. Save session information
############################

capture.output(
  sessionInfo(),
  file = file.path(
    output_directory,
    "R_sessionInfo.txt"
  )
)

message(
  "Analysis completed. Results were saved in: ",
  normalizePath(output_directory)
)