#' TCGA Expression Data Processing

#'

#' This function processes expression data and phenotype information, separates tumor and normal samples,

#' and saves the results into different files. It's specifically designed for data obtained from TCGA.

#'

#' @param counts_file_path File path to the counts data (usually in the form of a large matrix with gene expression data).

#' @param gene_probes_file_path File path containing the gene probes data.

#' @param phenotype_file_path File path to the phenotype data, which includes various sample attributes.

#' @param output_file_path Path where the output files, distinguished between tumor and normal, will be saved.

#'

#' @return A list containing matrices for tumor and normal expression data.

#'

#' @note IMPORTANT: This function assumes that the input files follow a specific format and structure, typically found in TCGA data releases.

#' Users should verify their data's compatibility. Additionally, the function does not perform error checking on the data's content,

#' which users should handle through proper preprocessing.

#'

#' @note CRITICAL: The 'output_file_path' parameter must end with '.rds' to be properly recognized by the function. It is also highly recommended

#' that the path includes specific identifiers related to the target samples, as the function will create further subdivisions in the specified

#' path for tumor or normal tissues. Please structure the 'output_file_path' following this pattern: './your_directory/your_sample_type.exp.rds'.

#'

#' @importFrom dplyr distinct filter

#' @importFrom utils read.table

#' @importFrom rlang .data

#' @export

#' @author Dongyue Yu

#'

#' @examples

#' counts_file <- system.file("extdata", "TCGA-SKCM.htseq_counts_test.tsv", package = "TransProR")

#' gene_probes_file <- system.file("extdata",

#'                                 "TCGA_gencode.v22.annotation.gene.probeMap_test",

#'                                 package = "TransProR")

#' phenotype_file <- system.file("extdata", "TCGA-SKCM.GDC_phenotype_test.tsv", package = "TransProR")

#' ouput_file <- file.path(tempdir(), "SKCM_Skin_TCGA_exp_test.rds")

#'

#' SKCM_exp <- get_tcga_exp(

#'   counts_file_path = counts_file,

#'   gene_probes_file_path = gene_probes_file,

#'   phenotype_file_path = phenotype_file,

#'   output_file_path = ouput_file

#' )

#' head(SKCM_exp[["tumor_tcga_data"]])[1:5, 1:5]

#' head(SKCM_exp[["normal_tcga_data"]], n = 10) # Because there is only one column.

get_tcga_exp <- function(counts_file_path,

                          gene_probes_file_path,

                          phenotype_file_path,

                          output_file_path) {

  # Load expression matrix

  # count_data <- data.table::fread(counts_file_path, header = TRUE, sep = '\t', data.table = FALSE)

  count_data <- utils::read.table(counts_file_path,

                          header = TRUE,

                          sep = '\t',

                          stringsAsFactors = FALSE,

                          check.names = FALSE)

  # Load gene ID conversion information

  # gene_probes <- data.table::fread(gene_probes_file_path, header = TRUE, sep = '\t', data.table = FALSE)

  gene_probes <- utils::read.table(gene_probes_file_path,

                            header = TRUE,

                            sep = '\t',

                            stringsAsFactors = FALSE,

                            check.names = FALSE)

  # Keep only necessary columns

  gene_probes <- gene_probes[, c(1, 2)]

  # Merge gene ID information with expression matrix

  count_probes_merged <- merge(gene_probes, count_data, by.x = "id", by.y = "Ensembl_ID")

  # Remove duplicates

  count_data_unique <- dplyr::distinct(count_probes_merged, .data$gene, .keep_all = TRUE)

  # Set gene names as row names

  rownames(count_data_unique) <- count_data_unique$gene

  count_data_final <- count_data_unique[, -c(1,2)]  # Remove extra columns

  # Load clinical information

  # phenotype_data <- data.table::fread(phenotype_file_path, header = TRUE, sep = '\t', data.table = FALSE)

  # Load clinical information with proper column types

  phenotype_data <- utils::read.table(phenotype_file_path,

                              header = TRUE,

                              sep = '\t',

                              stringsAsFactors = FALSE,

                              check.names = FALSE,

                              colClasses = list(

                                withdrawn = "logical",

                                releasable.project = "logical",

                                is_ffpe.samples = "logical",

                                oct_embedded.samples = "logical"

                              ))

  rownames(phenotype_data) <- phenotype_data$submitter_id.samples  # Set sample names as row names

  # Check first if there are any "Metastatic", "Primary Tumor", or "normal" samples in your data

  table(phenotype_data$sample_type.samples)

  # Create a data frame for tumor samples

  tumor_samples <- phenotype_data %>%

    dplyr::filter(grepl("Metastatic|Primary Tumor", .data$sample_type.samples, ignore.case = TRUE))

  # Check for 'Primary Tumor' or 'Metastatic' samples

  if (nrow(tumor_samples) == 0) {

    message("No 'Primary Tumor' or 'Metastatic' samples found.\n")

  } else {

    message("Number of 'Primary Tumor' or 'Metastatic' samples: ", nrow(tumor_samples), "\n")

  }

  # Create a data frame for normal samples

  normal_samples <- phenotype_data %>%

    dplyr::filter(grepl("normal", .data$sample_type.samples, ignore.case = TRUE))

  if (nrow(normal_samples) == 0) {

    message("No 'normal' samples found.\n")

  } else {

    message("Number of 'normal' samples:", nrow(normal_samples), "\n")

  }

  # Get the intersection of clinical information and expression matrix

  tumor_common_samples <- intersect(rownames(tumor_samples), colnames(count_data_final))

  # Extract corresponding expression matrix

  tumor_expression_data <- count_data_final[, tumor_common_samples]

  # Get the intersection of clinical information and expression matrix

  normal_common_samples <- intersect(rownames(normal_samples), colnames(count_data_final))

  # Extract corresponding expression matrix

  normal_expression_data <- count_data_final[, normal_common_samples, drop = FALSE]

  # Save results

  tumor_output_path <- gsub("\\.rds$", "_tumor.rds", output_file_path)

  normal_output_path <- gsub("\\.rds$", "_normal.rds", output_file_path)

  saveRDS(tumor_expression_data, file = tumor_output_path)

  saveRDS(normal_expression_data, file = normal_output_path)

  result <- list(

    tumor_tcga_data = tumor_expression_data,

    normal_tcga_data = normal_expression_data

  )

  return(result)

}