# Description: This file contains the wrapper functions for the methods that
# are used to compute the different layers of the multilayer network. The
# functions are called from the compute_*_network functions in layers.R
# For now, only the compute_atac_peak_network function has wrapper functions
# for the different methods. The other methods are still directly implemented
# in the compute_*_network functions in layers.R

#' @title Cicero wrapper function for the compute_atac_peak_network function
#'
#' @description This function is a wrapper for the compute_atac_peak_network
#' function in layers.R. It computes the peak network from scATAC-seq data
#' using Cicero. It returns a data frame with the peak network. The data frame
#' also contains the coaccess score for each edge. The coaccess score is the
#' probability that two peaks are accessible in the same cell. The coaccess 
#' score is computed by Cicero. Edges are filtered based on the coaccess score.
#' Only edges with a coaccess score > threshold are kept.
#'
#' @param hummus A hummus object
#' @param atac_assay The name of the assay containing the scATAC-seq data
#' @param genome The genome object
#' @param window The window size used by Cicero to compute the coaccess score
#' @param number_cells_per_clusters The number of cells per cluster used by
#' Cicero to compute the coaccess score
#' @param sample_num The number of samples used by Cicero to compute the
#' coaccess score
#' @param seed The seed used by Cicero to compute the coaccess score
#' @param verbose The verbosity level
#' @param threshold The threshold used to filter edges based on the coaccess
#' score
#' @param reduction_method The method used by monocle3 to reduce the dimension
#' of the scATAC-seq data before defining the pseudocells. The default is UMAP.
#'
#' @return A data frame containing the peak network
#' @export
#'
run_cicero_wrapper <- function(
    hummus,
    atac_assay,
    genome,
    window,
    number_cells_per_clusters,
    sample_num,
    seed,
    verbose,
    threshold,
    reduction_method = "UMAP"
    ) {

    # functions that need to be renamed
    int_elementMetadata <- SingleCellExperiment::int_elementMetadata
    counts <- SingleCellExperiment::counts

    # obtain chromosome sizes
    chromosome_sizes <- data.frame(V1 = genome@seqinfo@seqnames,
                                   V2 = genome@seqinfo@seqlengths)

    # Get scATAC-seq data
    scATAC <- as.matrix(hummus@assays[[atac_assay]]@counts)
    # Matrix to edgelist
    acc <- reshape2::melt(scATAC)
    colnames(acc) <- c("V1", "V2", "V3")

    # Prepare cicero input
    input_cds <- cicero::make_atac_cds(acc, binarize = TRUE) # Create CDS object
    set.seed(seed)
    # It is required that there is no empty cell
    if (length(which(colSums(as.matrix(monocle3::exprs(input_cds))) == 0)) == 0
    ) {
    # Calculating size factors using default method = mean-geometric-mean-total
      input_cds <- monocle3::estimate_size_factors(input_cds)
      # Preprocessing using LSI
      input_cds <- monocle3::preprocess_cds(input_cds, method = "LSI")
      # Dimensionality reduction using UMAP
      input_cds <- monocle3::reduce_dimension(
                                    input_cds,
                                    reduction_method = reduction_method,
                                    preprocess_method = "LSI")
    } else {
      print("Error: there is at least one cell with no signal.")
    }
    # Get reduced (UMAP) coordinates
    umap_coords <- SingleCellExperiment::reducedDims(input_cds)$UMAP
    # Compute pseudocells
    cicero_cds <- cicero::make_cicero_cds(
      input_cds,  # Create a Cicero CDS object
      reduced_coordinates = umap_coords,
      k = number_cells_per_clusters,  #number neighbors/ Default = 50
      summary_stats = NULL,         # Default
      size_factor_normalize = TRUE, # Default
      silent = FALSE)               # Default

    cicero <- cicero::run_cicero(
      cds = cicero_cds, # Infer peak-links
      genomic_coords = chromosome_sizes,
      window = window,             # Default = 5e+05
      silent = FALSE,             # Default = FALSE
      sample_num = sample_num) # Default = 100

    # Remove NAs, double edges, and edges with coaccess score <=0
    # Check for coaccess = NA
    if (length(which(is.na(cicero$coaccess))) > threshold) {
      cicero <- cicero[which(!is.na(cicero$coaccess)), ]  # Remove NAs
    }
    cicero$temp <- NA  # Helper column to check and remove double edges
    my_cols <- which(as.character(cicero$Peak1) <= as.character(cicero$Peak2))
    cicero$temp[my_cols] <- paste(cicero$Peak1[my_cols],
                                  cicero$Peak2[my_cols],
                                  sep = ";")

    my_cols <- which(as.character(cicero$Peak1) > as.character(cicero$Peak2))
    cicero$temp[my_cols] <- paste(cicero$Peak2[my_cols],
                                  cicero$Peak1[my_cols],
                                  sep = ";")

    # Sort table according to temp-column (each entry appears double)
    cicero <- cicero[with(cicero, order(temp, decreasing = TRUE)), ]
    rownames(cicero) <- c(1:dim(cicero)[1])
    A <- as.character(cicero$Peak1[seq(1, dim(cicero)[1], 2)])
    Anum <- round(cicero$coaccess[seq(1, dim(cicero)[1], 2)], 10)
    B <- as.character(cicero$Peak2[seq(2, dim(cicero)[1], 2)])
    Bnum <- round(cicero$coaccess[seq(2, dim(cicero)[1], 2)], 10)
    # length(which(A==B & Anum==Bnum)) 
    # Each edge appears twice with same coaccess score (rounded to 10 digits)
    cicero <- cicero[seq(1, dim(cicero)[1], 2), ] # Remove double edges
    cicero$temp <- NULL # Remove helper column
    cicero <- cicero[with(cicero, order(cicero$coaccess,
                                        decreasing = TRUE)), ]  # Sort
    rownames(cicero) <- c(1:dim(cicero)[1])
    cicero$Peak1 <- gsub("_", "-", cicero$Peak1)
    # Peak names 2x"-" to match bipartites
    cicero$Peak2 <- gsub("_", "-", cicero$Peak2)
    # Peak names 2x"-" to match bipartites ? 2x"-" or 2x"_"

    peak_network <- cicero[which(cicero$coaccess > threshold), ]
    # Remove edges with coaccess score <= threshold

    if (verbose > 0) {
      cat("\n", dim(peak_network)[1], "peak edges with a coaccess score >",
          threshold, "were found.\n")
    }

    # Return peak network including edges with positive coaccess score
    return(peak_network)
}


#' @title Omnipath wrapper function for the compute_tf_network function
#' @description This function is a wrapper for the compute_tf_network function
#' in layers.R. It computes the TF network from using Omnipath database.
#' It returns a data frame with the TF network. The data frame is not weighted
#' and does not contain scores for the edges.
#' @param hummus A hummus object
  # Get TF-TF interactions from Omnipath
run_omnipath_wrapper <- function(
  hummus = hummus,
  organism = organism,
  tfs = tfs,
  gene_assay = gene_assay,
  source_target = source_target,
  verbose = 1) {

  TF_PPI <- OmnipathR::import_post_translational_interactions(
    organism = organism, partners = tfs, source_target = source_target
  )

  if (verbose > 0) {
    cat("\tNumber of edges from Omnipath:", nrow(TF_PPI),
    "\nWill now be filtered to only those corresponding to specified tfs")
  }

  if (is.na(tfs)) {
    # Get tfs list
    tfs <- get_tfs(hummus = hummus,
              assay = gene_assay,
              store_tfs = FALSE,
              output_file = NULL,
              verbose = verbose)
  } else if (typeof(tfs) != "character") {
      stop("'tfs' argument needs to be a vector of characters
      (e.g.: c('MYC', 'JAK1')).")
  }

  # add filtering if element is not a TF expressed in the dataset
  if (source_target == "AND") {
    TF_PPI <- TF_PPI[which(TF_PPI$source_genesymbol %in% tfs &
                           TF_PPI$target_genesymbol %in% tfs), ]
  } else if (source_target == "OR") {
    TF_PPI <- TF_PPI[which(TF_PPI$source_genesymbol %in% tfs |
                           TF_PPI$target_genesymbol %in% tfs), ]
  }
  # Get only source and target columns
  tf_network <- TF_PPI[, c(3, 4)]

  # Convert to data.frame from tibble
  tf_network <- as.data.frame(tf_network)

  # Check if there is any TF-TF edges otherwise add a fake node
  # and connect all TFs to it (to allow HuMMuS to run without impacting result)
  if (nrow(tf_network) == 0) {
    cat("No TF-TF edges from Omnipath for the given parameters.
        You can try to change the source_target parameter to 'OR' to get
        TF-other protein interactions. Or try to import a network  
        computed externally. Right now, a network with all TFs connected
        to a fake node is created, for HuMMuS analysis.\n It has no biological
        meaning but will allow to run the pipeline as if no edges were present.
        \n")
    tf_network <- run_tf_null_wrapper(
      hummus = hummus,
      organism = organism,
      tfs = tfs,
      gene_assay = gene_assay,
      verbose = )

  }
  return(tf_network)
}

#' @title tf_null wrapper function for the tf_network function
#' @description This function is a wrapper for the tf_network function
#' 
#' @param hummus A hummus object
#' 
#' 
run_tf_null_wrapper <- function(
  hummus = hummus,
  organism = organism,
  tfs = tfs,
  gene_assay = gene_assay,
  verbose = 1) {
    
  if (verbose > 0) {
    cat("Creating a fake TF network with all TFs connected to a fake node.\n")
  }

  if (is.na(tfs)) {
    # Get tfs list
    tfs <- get_tfs(hummus = hummus,
              assay = gene_assay,
              store_tfs = FALSE,
              output_file = NULL,
              verbose = verbose)
  } else if (typeof(tfs) != "character") {
      stop("'tfs' argument needs to be a vector of characters
      (e.g.: c('MYC', 'JAK1')).")
  }


  FAKE_NODE <- "fake_node"
  tf_network <- data.frame(source = c(), target = c())
  for (tf in tfs) {
    tf_network <- rbind(tf_network, data.frame(source = tf, target = FAKE_NODE))
  }
  return(tf_network)
}