query_general_id <- function() {
  tcga_clinical_fine <- load_data("tcga_clinical_fine")
  pcawg_info_fine <- load_data("pcawg_info_fine")
  ccle_info_fine <- load_data("ccle_info_fine")
  #### TCGA general value: index/immune/pathway/phenotype
  ## tumor_index
  tcga_index_value <- list(
    `Tumor Purity` = load_data("tcga_purity") %>% dplyr::rename("Sample" = "sample"),
    `Tumor Stemness` = load_data("tcga_stemness") %>% dplyr::rename("Sample" = "sample"),
    `Tumor Mutation Burden` = load_data("tcga_tmb") %>% dplyr::rename("Sample" = "Tumor_Sample_ID"),
    `Microsatellite Instability` = load_data("tcga_MSI"),
    `Genome Instability` = load_data("tcga_genome_instability") %>% dplyr::rename("Sample" = "sample")
  )
  # colnames(tcga_index_value$`Tumor Mutation Burden`)[3] = "sample"
  tcga_index_value$`Microsatellite Instability` <- load_data("tcga_gtex") %>%
    dplyr::rename("Sample" = "sample") %>%
    dplyr::mutate(Barcode = stringr::str_sub(.data$Sample, 1, 12)) %>%
    dplyr::select("Barcode", "Sample") %>%
    dplyr::inner_join(tcga_index_value$`Microsatellite Instability`, by = "Barcode")

  ## tumor_immune
  tcga_TIL <- load_data("tcga_TIL")
  colnames(tcga_TIL)[1] <- "Sample"
  dim(tcga_TIL)
  # [1] 11070   120
  tcga_immune_value <- list(
    `CIBERSORT` = tcga_TIL %>% dplyr::select("Sample", ends_with("CIBERSORT")) %>% dplyr::rename_with(~ gsub("_CIBERSORT", "", .x)),
    `CIBERSORT-ABS` = tcga_TIL %>% dplyr::select("Sample", ends_with("CIBERSORT-ABS")) %>% dplyr::rename_with(~ gsub("_CIBERSORT-ABS", "", .x)),
    `EPIC` = tcga_TIL %>% dplyr::select("Sample", ends_with("EPIC")) %>% dplyr::rename_with(~ gsub("_EPIC", "", .x)),
    `MCPCOUNTER` = tcga_TIL %>% dplyr::select("Sample", ends_with("MCPCOUNTER")) %>% dplyr::rename_with(~ gsub("_MCPCOUNTER", "", .x)),
    `QUANTISEQ` = tcga_TIL %>% dplyr::select("Sample", ends_with("QUANTISEQ")) %>% dplyr::rename_with(~ gsub("_QUANTISEQ", "", .x)),
    `TIMER` = tcga_TIL %>% dplyr::select("Sample", ends_with("TIMER")) %>% dplyr::rename_with(~ gsub("_TIMER", "", .x)),
    `XCELL` = tcga_TIL %>% dplyr::select("Sample", ends_with("XCELL")) %>% dplyr::rename_with(~ gsub("_XCELL", "", .x))
  )

  ## tumor_pathway
  tcga_PW <- load_data("tcga_PW") %>%
    as.data.frame() %>%
    tibble::rownames_to_column("Sample")
  dim(tcga_PW)
  # [1] 10534   501
  tcga_pathway_value <- list(
    `HALLMARK` = tcga_PW %>% dplyr::select("Sample", contains("HALLMARK")) %>% dplyr::rename_with(~ gsub("HALLMARK_", "", .x)),
    `KEGG` = tcga_PW %>% dplyr::select("Sample", contains("KEGG")) %>% dplyr::rename_with(~ gsub("KEGG_", "", .x)),
    `IOBR` = tcga_PW %>% dplyr::select("Sample", contains("IOBR")) %>% dplyr::rename_with(~ gsub("IOBR_", "", .x))
  )

  ## tumor_phenotype (user upload custom metadata)
  tcga_phenotype_value <- list(
    `Clinical Phenotype` = tcga_clinical_fine,
    `Custom metadata` = NULL
  )


  #### PCAWG general value: index/immune/pathway/phenotype
  ## tumor_index
  pcawg_index_value <- list(
    `Tumor Purity` = load_data("pcawg_purity") %>%
      dplyr::filter(.data$icgc_specimen_id %in% pcawg_info_fine$Sample) %>%
      dplyr::rename("Sample" = "icgc_specimen_id")
  )

  ## tumor_immune
  pcawg_TIL <- load_data("pcawg_TIL")
  colnames(pcawg_TIL)[1] <- "Sample"
  dim(pcawg_TIL)
  # [1] 1466   120
  pcawg_immune_value <- list(
    `CIBERSORT` = pcawg_TIL %>% dplyr::select("Sample", ends_with("CIBERSORT")) %>% dplyr::rename_with(~ gsub("_CIBERSORT", "", .x)),
    `CIBERSORT-ABS` = pcawg_TIL %>% dplyr::select("Sample", ends_with("CIBERSORT-ABS")) %>% dplyr::rename_with(~ gsub("_CIBERSORT-ABS", "", .x)),
    `EPIC` = pcawg_TIL %>% dplyr::select("Sample", ends_with("EPIC")) %>% dplyr::rename_with(~ gsub("_EPIC", "", .x)),
    `MCPCOUNTER` = pcawg_TIL %>% dplyr::select("Sample", ends_with("MCPCOUNTER")) %>% dplyr::rename_with(~ gsub("_MCPCOUNTER", "", .x)),
    `QUANTISEQ` = pcawg_TIL %>% dplyr::select("Sample", ends_with("QUANTISEQ")) %>% dplyr::rename_with(~ gsub("_QUANTISEQ", "", .x)),
    `TIMER` = pcawg_TIL %>% dplyr::select("Sample", ends_with("TIMER")) %>% dplyr::rename_with(~ gsub("_TIMER", "", .x)),
    `XCELL` = pcawg_TIL %>% dplyr::select("Sample", ends_with("XCELL")) %>% dplyr::rename_with(~ gsub("_XCELL", "", .x))
  )

  ## tumor_pathway
  pcawg_PW <- load_data("pcawg_PW") %>%
    as.data.frame() %>%
    tibble::rownames_to_column("Sample")
  dim(pcawg_PW)
  # [1] 1466   501
  pcawg_pathway_value <- list(
    `HALLMARK` = pcawg_PW %>% dplyr::select("Sample", contains("HALLMARK")) %>% dplyr::rename_with(~ gsub("HALLMARK_", "", .x)),
    `KEGG` = pcawg_PW %>% dplyr::select("Sample", contains("KEGG")) %>% dplyr::rename_with(~ gsub("KEGG_", "", .x)),
    `IOBR` = pcawg_PW %>% dplyr::select("Sample", contains("IOBR")) %>% dplyr::rename_with(~ gsub("IOBR_", "", .x))
  )

  pcawg_phenotype_value <- list(
    `Clinical Phenotype` = pcawg_info_fine,
    `Custom metadata` = NULL
  )

  #### CCLE general value: index/immune/pathway/phenotype
  ## tumor_index
  ccle_index_value <- list(
    `Tumor Purity` = load_data("ccle_absolute") %>%
      dplyr::rename("Sample" = "Cell Line") %>%
      dplyr::filter(.data$Sample %in% ccle_info_fine$Sample)
  )
  ccle_phenotype_value <- list(
    `Clinical Phenotype` = ccle_info_fine,
    `Custom metadata` = NULL
  )

  tcga_value_option <- list(
    "Tumor index" = tcga_index_value, "Immune Infiltration" = tcga_immune_value,
    "Pathway activity" = tcga_pathway_value, "Phenotype data" = tcga_phenotype_value
  )
  pcawg_value_option <- list(
    "Tumor index" = pcawg_index_value, "Immune Infiltration" = pcawg_immune_value,
    "Pathway activity" = pcawg_pathway_value, "Phenotype data" = pcawg_phenotype_value
  )
  ccle_value_option <- list("Tumor index" = ccle_index_value, "Phenotype data" = ccle_phenotype_value)

  #### TCGA general id: molecule/index/immune/pathway/phenotype
  ## molecule_profile
  # pancan_identifiers <- readRDS(
  #   system.file(
  #     "extdata", "pancan_identifier_list.rds",
  #     package = "UCSCXenaShiny"
  #   )
  # )
  tcga_id <- load_data("v2_tcga_id")
  tcga_id.list <- split(tcga_id$L3, tcga_id$L2)

  tcga_molecule_id <- list(
    `mRNA Expression` = list(all = tcga_id.list[["Gene"]], default = "TP53"),
    # `Transcript Expression` = list(all = load_data("transcript_identifier"), default = "ENST00000000233"),
    `Transcript Expression` = list(all = tcga_id.list[["Transcript"]], default = "ENST00000000233"),
    `DNA Methylation` = list(all = tcga_id.list[["Gene"]], default = "TP53"),
    `Protein Expression` = list(all = tcga_id.list[["Protein"]], default = "P53"),
    `miRNA Expression` = list(all = tcga_id.list[["miRNA"]], default = "hsa-miR-769-3p"),
    `Mutation status` = list(all = tcga_id.list[["Gene"]], default = "TP53"),
    `Copy Number Variation` = list(all = tcga_id.list[["Gene"]], default = "TP53")
  )

  ## tumor_index
  tcga_index_id <- list(
    `Tumor Purity` = list(all = sort(colnames(tcga_index_value$`Tumor Purity`)[3:7]), default = "ESTIMATE"),
    `Tumor Stemness` = list(all = sort(colnames(tcga_index_value$`Tumor Stemness`)[2:6]), default = "RNAss"),
    `Tumor Mutation Burden` = list(all = sort(colnames(tcga_index_value$`Tumor Mutation Burden`)[4:5]), default = "Non_silent_per_Mb"),
    `Microsatellite Instability` = list(
      all = sort(setdiff(
        colnames(tcga_index_value$`Microsatellite Instability`)[3:21],
        c("MSI_intronic", "MSI_intronic_profiled", "MSI_noncoding", "MSI_noncoding_profiled")
      )),
      default = "Total_nb_MSI_events"
    ),
    `Genome Instability` = list(all = sort(colnames(tcga_index_value$`Genome Instability`)[2:6]), default = "ploidy")
  )

  ## tumor_immune
  tcga_immune_id <- lapply(tcga_immune_value, function(x) {
    list(all = sort(colnames(x)[-1]))
  })
  tcga_immune_id$`CIBERSORT`$default <- "Monocyte"
  tcga_immune_id$`CIBERSORT-ABS`$default <- "Monocyte"
  tcga_immune_id$`EPIC`$default <- "Macrophage"
  tcga_immune_id$`MCPCOUNTER`$default <- "Monocyte"
  tcga_immune_id$`QUANTISEQ`$default <- "Monocyte"
  tcga_immune_id$`TIMER`$default <- "Monocyte"
  tcga_immune_id$`XCELL`$default <- "Monocyte"

  ## tumor_pathway
  tcga_pathway_id <- lapply(tcga_pathway_value, function(x) {
    list(all = sort(colnames(x)[-1]))
  })
  tcga_pathway_id$HALLMARK$default <- "APOPTOSIS"
  tcga_pathway_id$KEGG$default <- "CELL_CYCLE"
  tcga_pathway_id$IOBR$default <- "Biotin_Metabolism"

  ## tumor_phenotype
  tcga_phenotype_id <- list(
    `Clinical Phenotype` = list(all = sort(colnames(tcga_phenotype_value$`Clinical Phenotype`[3:8])), default = "Code"),
    `Custom metadata` = list(all = NULL, default = NULL)
  )


  ## merge
  tcga_id_option <- list(
    "Molecular profile" = tcga_molecule_id,
    "Tumor index" = tcga_index_id,
    "Immune Infiltration" = tcga_immune_id,
    "Pathway activity" = tcga_pathway_id,
    "Phenotype data" = tcga_phenotype_id
  )

  #### PCAWG ID
  pcawg_id_option <- tcga_id_option
  # pcawg_id_referrence = load_data("pcawg_identifier")
  pcawg_id_referrence <- load_data("v2_pcawg_id")
  pcawg_id_option$`Molecular profile` <- list(
    `mRNA Expression` = list(all = sort(pcawg_id_referrence$L3[pcawg_id_referrence$L2 == "Gene"]), default = "TP53"),
    `Promoter Activity` = list(all = sort(pcawg_id_referrence$L3[pcawg_id_referrence$L2 == "Promoter"]), default = "prmtr.1"),
    `Gene Fusion` = list(all = sort(pcawg_id_referrence$L3[pcawg_id_referrence$L2 == "Fusion"]), default = "SAMD11"),
    `miRNA Expression` = list(all = sort(pcawg_id_referrence$L3[pcawg_id_referrence$L2 == "miRNA"]), default = "hsa-let-7a-2-3p"),
    `APOBEC Mutagenesis` = list(all = sort(pcawg_id_referrence$L3[pcawg_id_referrence$L2 == "Muta"]), default = "A3A_or_A3B")
  )
  pcawg_id_option$`Tumor index` <- list(
    `Tumor Purity` = list(
      all = c("purity", "ploidy", "purity_conf_mad", "wgd_status", "wgd_uncertain"),
      default = "purity"
    )
  )
  pcawg_id_option$`Phenotype data` <- list(
    `Clinical Phenotype` = list(
      all = c("Age", "Gender", "Type"),
      default = "Age"
    ),
    `Custom metadata` = list(all = NULL, default = NULL)
  )

  #### CCLE ID
  ccle_id_option <- list()
  # ccle_id_referrence = load_data("ccle_identifier")
  ccle_id_referrence <- load_data("v2_ccle_id")
  ccle_id_option$`Molecular profile` <- list(
    `mRNA Expression` = list(all = sort(ccle_id_referrence$L3[ccle_id_referrence$L2 == "Gene"]), default = "TP53"),
    `Protein Expression` = list(all = sort(ccle_id_referrence$L3[ccle_id_referrence$L2 == "Protein"]), default = "14-3-3_beta"),
    `Copy Number Variation` = list(all = sort(ccle_id_referrence$L3[ccle_id_referrence$L2 == "Gene"]), default = "TP53"),
    `Mutation status` = list(all = sort(ccle_id_referrence$L3[ccle_id_referrence$L2 == "Gene"]), default = "TP53")
  )
  ccle_id_option$`Tumor index` <- list(
    `Tumor Purity` = list(
      all = c("Purity", "Ploidy", "Genome Doublings", "Lineage"),
      default = "Purity"
    )
  )
  ccle_id_option$`Immune Infiltration` <- list(NULL)
  ccle_id_option$`Pathway activity` <- list(NULL)
  ccle_id_option$`Phenotype data` <- list(
    `Clinical Phenotype` = list(
      all = c("Gender", "Histology", "Type"),
      default = "Gender"
    ),
    `Custom metadata` = list(all = NULL, default = NULL)
  )

  res_merge <- list(
    value = list(tcga_value_option, pcawg_value_option, ccle_value_option),
    id = list(tcga_id_option, pcawg_id_option, ccle_id_option)
  )

  return(res_merge)
}







#' Download data for shiny general analysis
#'
#' @param L1 level 1  main datatype
#' @param L2 level 2  sub datatype
#' @param L3 level 3  identifier
#' @param opt_pancan      molecular datasets parameters
#' @param custom_metadata user customized metadata
#' @param database        one of c("toil","pcawg","ccle")
#' @param tpc_value_nonomics     non-omics matrix data of one database
#'
#' @examples
#' \dontrun{
#' general_value_id <- UCSCXenaShiny:::query_general_id()
#' tcga_value_option <- general_value_id[["value"]][[1]]
#' tcga_index_value <- tcga_value_option[["Tumor index"]]
#' tcga_immune_value <- tcga_value_option[["Immune Infiltration"]]
#' tcga_pathway_value <- tcga_value_option[["Pathway activity"]]
#' tcga_phenotype_value <- tcga_value_option[["Phenotype data"]]
#' 
#' clinical_phe = tcga_phenotype_value[["Clinical Phenotype"]]
#' x_data = UCSCXenaShiny:::query_general_value(
#'            "Molecular profile", "mRNA Expression", "TP53", "toil",
#'            tcga_index_value, tcga_immune_value, tcga_pathway_value,
#'            clinical_phe)
#'
#' y_data = UCSCXenaShiny:::query_general_value(
#'            "Immune Infiltration", "CIBERSORT", "Monocyte", "toil",
#'            tcga_index_value, tcga_immune_value, tcga_pathway_value,
#'            clinical_phe)
#' }


query_general_value <- function(L1, L2, L3,
                                database = c("toil", "pcawg", "ccle"),
                                tpc_value_nonomics = NULL,
                                opt_pancan = NULL, custom_metadata = NULL) {
  database <- match.arg(database)

  if (L1 == "Molecular profile") {
    # L2 = "mRNA Expression"
    # L3 = "TP53"
    L2_label <- switch(L2,
      `mRNA Expression` = "mRNA",
      `Transcript Expression` = "transcript",
      `DNA Methylation` = "methylation",
      `Protein Expression` = "protein",
      `miRNA Expression` = "miRNA",
      `Mutation status` = "mutation",
      `Copy Number Variation` = "cnv",
      `Promoter Activity` = "promoter",
      `Gene Fusion` = "fusion",
      `APOBEC Mutagenesis` = "APOBEC"
    )
    if (is.null(opt_pancan)) {
      opt_pancan <- .opt_pancan
    }
    x_data <- query_pancan_value(
      molecule = L3,
      data_type = L2_label,
      database = database,
      opt_pancan = opt_pancan
    )
    if (database == "ccle" & L2 == "Mutation status") {
      x_data <- x_data[, c("sampleID", "genes")]
      x_data <- x_data %>%
        dplyr::select("sampleID", "genes") %>%
        dplyr::mutate(genes = ifelse(is.na(.data$genes), 0, 1)) %>%
        dplyr::arrange(.data$sampleID, desc(.data$genes)) %>%
        dplyr::distinct(.data$sampleID, .keep_all = T) %>%
        dplyr::mutate(id = L3, .before = 1) %>%
        dplyr::rename("Sample" = "sampleID", "value" = "genes") %>%
        dplyr::mutate(level2 = "Mutation status")
    } else {
      if (is.list(x_data)) x_data <- x_data[[1]]
      x_data <- data.frame(
        "id" = L3,
        "Sample" = names(x_data), "value" = as.numeric(x_data),
        "level2" = L2
      )
    }
  } else {
    if (is.null(tpc_value_nonomics)) {
      if(database=="toil"){
        tpc_value_nonomics <- load_data(paste0("v2_tcga_value_nonomics"))
      } else {
        tpc_value_nonomics <- load_data(paste0("v2_", database, "_value_nonomics"))
      }
    }
    if (L2 != "Custom metadata") {
      L1_label <- switch(L1,
        `Tumor index` = "Index",
        `Immune Infiltration` = "Immune",
        `Pathway activity` = "Pathway",
        `Phenotype data` = "Phenotype"
      )
      L2_label <- switch(L2,
        # "Tumor index"
        `Tumor Purity` = "Purity",
        `Tumor Stemness` = "Stem",
        `Tumor Mutation Burden` = "TMB",
        `Microsatellite Instability` = "MSI",
        `Genome Instability` = "GI",
        # "Immune Infiltration"
        `CIBERSORT` = "CIB",
        `CIBERSORT-ABS` = "CIB.ABS",
        `EPIC` = "EPIC",
        `MCPCOUNTER` = "MCP",
        `QUANTISEQ` = "Quant",
        `TIMER` = "TIMER",
        `XCELL` = "XCELL",
        # "Pathway activity"
        `HALLMARK` = "HM",
        `KEGG` = "KEGG",
        `IOBR` = "IOBR",
        # "Phenotype data"
        `Clinical Phenotype` = "Clinical"
      )
      nonomics_id <- paste(L1_label, L2_label, L3, sep = "+")
      stopifnot(nonomics_id %in% colnames(tpc_value_nonomics))
      x_data <- tpc_value_nonomics[, c("Sample", nonomics_id)]
    } else {
      stopifnot(!is.null(custom_metadata))
      x_data <- custom_metadata[, c("Sample", L3)] %>%
        dplyr::filter(.data$Sample %in% tpc_value_nonomics$Sample)
    }

    colnames(x_data) <- c("Sample", "value")
    x_data <- x_data %>%
      dplyr::mutate(id = L3, .before = 1) %>%
      dplyr::mutate(level2 = L2)
  }
  x_data <- as.data.frame(x_data)
  # 检查重复样本
  if (is.numeric(x_data$value)) {
    x_data <- x_data %>%
      # dplyr::group_by(.data$Sample, .data$id, .data$level2) %>%
      # dplyr::summarise(value = mean(.data$value)) %>%
      dplyr::filter(!is.na(.data$value)) %>%
      as.data.frame()
  }
  x_data
  # id                  Sample value          level2
  # 1 TP53 GTEX-S4Q7-0003-SM-3NM8M 4.785 mRNA Expression
  # 2 TP53         TCGA-19-1787-01 5.887 mRNA Expression
  # 3 TP53         TCGA-S9-A7J2-01 5.517 mRNA Expression
  # 4 TP53 GTEX-QV31-1626-SM-2S1QC 4.431 mRNA Expression
  # 5 TP53         TCGA-G3-A3CH-11 2.382 mRNA Expression
  # 6 TP53         TCGA-B5-A5OE-01 5.765 mRNA Expression
}









#' Quick molecule analysis and report generation
#'
#' @inheritParams query_pancan_value
#' @param out_dir path to save analysis result and report, default is '.'
#' @param out_report logical value wheather to generate html report
#'
#' @return a list.
#' @export
mol_quick_analysis <- function(molecule, data_type, out_dir = ".", out_report = FALSE) {
  if (dir.exists(out_dir)) {
    cat("Warning: The file path already exists. Existing file may be overwritten.\n")
  } else {
    dir.create(out_dir)
  }

  print("##### Step1: Query the moleluce value... #####")
  mol_data <- query_pancan_value(molecule, data_type = data_type, database = "toil")
  if (is.list(mol_data)) mol_data <- mol_data[[1]]
  
  if (all(is.na(mol_data))) {
    message("All NAs returned, return NULL instead.")
    return(NULL)
  }

  print(paste0("=== ", "Clinical phenotype"))
  mol_data_df <- suppressMessages(data.frame(
    id = molecule,
    value = mol_data
  ) %>%
    tibble::rownames_to_column("Sample") %>%
    dplyr::inner_join(load_data("tcga_gtex"), by = c("Sample" = "sample")))

  print(paste0("---> ", "Tumor&Normal"))
  comp_TN <- suppressMessages(mol_data_df %>%
    dplyr::filter(!.data$tissue %in% c("MESO", "UVM")) %>%
    dplyr::group_by(.data$tissue) %>%
    dplyr::summarise(wilcox_p = stats::wilcox.test(value ~ type2)$p.value))
  comp_TN <- suppressMessages(mol_data_df %>%
    dplyr::filter(!.data$tissue %in% c("MESO", "UVM")) %>%
    dplyr::group_by(.data$tissue, .data$type2) %>%
    dplyr::summarise(value = mean(.data$value)) %>%
    tidyr::pivot_wider(names_from = "type2", values_from = "value") %>%
    dplyr::inner_join(comp_TN))

  print(paste0("---> ", "Age"))
  mol_data_Age <- suppressMessages(mol_data_df[, 1:4] %>%
    dplyr::inner_join(load_data("tcga_clinical_fine")[, c("Sample", "Age")]) %>%
    dplyr::filter(!is.na(.data$Age)) %>%
    dplyr::mutate(Age = ifelse(.data$Age >= 60, "Old", "Young")))
  valid_types <- suppressMessages(mol_data_Age %>%
    dplyr::group_by(.data$tissue) %>%
    dplyr::count(.data$Age) %>% dplyr::count(.data$tissue) %>%
    dplyr::filter(n == 2) %>% dplyr::pull(.data$tissue))
  mol_data_Age <- mol_data_Age[mol_data_Age$tissue %in% valid_types, ]
  comp_Age <- suppressMessages(mol_data_Age %>%
    dplyr::group_by(.data$tissue) %>%
    dplyr::summarise(wilcox_p = stats::wilcox.test(value ~ Age)$p.value))
  comp_Age <- suppressMessages(mol_data_Age %>%
    dplyr::group_by(.data$tissue, .data$Age) %>%
    dplyr::summarise(value = mean(.data$value)) %>%
    tidyr::pivot_wider(names_from = "Age", values_from = "value") %>%
    dplyr::inner_join(comp_Age))

  print(paste0("---> ", "Gender"))
  mol_data_Gender <- suppressMessages(mol_data_df[, 1:4] %>%
    dplyr::inner_join(load_data("tcga_clinical_fine")[, c("Sample", "Gender")]) %>%
    dplyr::filter(!is.na(.data$Gender)))
  valid_types <- suppressMessages(mol_data_Gender %>%
    dplyr::group_by(.data$tissue) %>%
    dplyr::count(.data$Gender) %>% dplyr::count(.data$tissue) %>%
    dplyr::filter(n == 2) %>% dplyr::pull(.data$tissue))
  mol_data_Gender <- mol_data_Gender[mol_data_Gender$tissue %in% valid_types, ]
  comp_Gender <- suppressMessages(mol_data_Gender %>%
    dplyr::group_by(.data$tissue) %>%
    dplyr::summarise(wilcox_p = stats::wilcox.test(value ~ Gender)$p.value))
  comp_Gender <- suppressMessages(mol_data_Gender %>%
    dplyr::group_by(.data$tissue, .data$Gender) %>%
    dplyr::summarise(value = mean(.data$value)) %>%
    tidyr::pivot_wider(names_from = "Gender", values_from = "value") %>%
    dplyr::inner_join(comp_Gender))

  print(paste0("---> ", "Stage"))
  mol_data_Stage <- suppressMessages(mol_data_df[, 1:4] %>%
    dplyr::inner_join(load_data("tcga_clinical_fine")[, c("Sample", "Stage_ajcc")]) %>%
    dplyr::filter(!is.na(.data$Stage_ajcc)))
  valid_types <- suppressMessages(mol_data_Stage %>%
    dplyr::group_by(.data$tissue) %>%
    dplyr::count(.data$Stage_ajcc) %>% dplyr::count(.data$tissue) %>%
    dplyr::filter(n >= 2) %>% dplyr::pull(.data$tissue))
  mol_data_Stage <- mol_data_Stage[mol_data_Stage$tissue %in% valid_types, ]
  comp_Stage <- suppressMessages(mol_data_Stage %>%
    dplyr::group_by(.data$tissue) %>%
    dplyr::summarise(aov_p = summary(stats::aov(value ~ Stage_ajcc))[[1]][1, 5]))
  comp_Stage <- suppressMessages(mol_data_Stage %>%
    dplyr::group_by(.data$tissue, .data$Stage_ajcc) %>%
    dplyr::summarise(value = mean(.data$value)) %>%
    tidyr::pivot_wider(names_from = "Stage_ajcc", values_from = "value") %>%
    dplyr::inner_join(comp_Stage))

  phe_res <- list(
    comp_TN = comp_TN %>%
      tidyr::pivot_longer(c("normal", "tumor")) %>% dplyr::rename("P.value" = "wilcox_p"),
    comp_Age = comp_Age %>%
      tidyr::pivot_longer(c("Old", "Young")) %>% dplyr::rename("P.value" = "wilcox_p"),
    comp_Gender = comp_Gender %>%
      tidyr::pivot_longer(c("FEMALE", "MALE")) %>% dplyr::rename("P.value" = "wilcox_p"),
    comp_Stage = comp_Stage %>%
      tidyr::pivot_longer(starts_with("Stage")) %>% dplyr::rename("P.value" = "aov_p")
  ) %>%
    do.call(rbind, .) %>%
    as.data.frame()


  tcga_clinical <- load_data("tcga_clinical")
  mol_data_df <- data.frame(
    id = molecule,
    value = mol_data
  ) %>%
    tibble::rownames_to_column("Sample") %>%
    dplyr::filter(.data$Sample %in% tcga_clinical$sample) %>%
    dplyr::mutate(type = tcga_clinical$type[match(.data$Sample, tcga_clinical$sample)]) %>%
    dplyr::filter(grepl("01$", .data$Sample)) %>% # -01 primary tumor
    dplyr::arrange(.data$type, .data$Sample)

  # 依据每种癌症的中位数进行分组
  if (data_type == "mutation") {
    mol_data_df <- mol_data_df %>%
      dplyr::mutate(group = ifelse(.data$value == 1, "Mutate", "Wild")) %>%
      as.data.frame()
  } else {
    mol_data_df <- mol_data_df %>%
      dplyr::group_by(.data$type) %>%
      dplyr::mutate(group = ifelse(.data$value > median(.data$value), "Higher", "Lower")) %>%
      as.data.frame()
  }

  general_value_id <- query_general_id()
  tcga_value_option <- general_value_id[["value"]][[1]]
  tcga_id_option <- general_value_id[["id"]][[1]]

  ## 相关性分析结果
  print("##### Step2: Execute correlation analysis...  #####")
  cor_res <- lapply(c("Tumor index", "Immune Infiltration", "Pathway activity"), function(L1_type) {
    # L1_type = "Pathway activity"
    print(paste0("=== ", L1_type))
    tcga_id <- tcga_id_option[[L1_type]]
    tcga_value <- tcga_value_option[[L1_type]]
    L2_types <- names(tcga_id)

    cor_L1 <- lapply(L2_types, function(L2_type) {
      # L2_type = L2_types[3]
      print(paste0("---> ", L2_type))
      L3_ids <- tcga_id[[L2_type]]$all
      value_df <- tcga_value[[L2_type]] %>% as.data.frame()
      L3_ids <- L3_ids[sapply(L3_ids, function(id) {
        class(value_df[, id])
      }) != "character"]

      value_df_long <- value_df[, c("Sample", L3_ids)] %>%
        tidyr::pivot_longer(-.data$Sample, names_to = "L3", values_to = "target")

      data_merge <- suppressMessages(dplyr::inner_join(mol_data_df, value_df_long) %>%
        dplyr::filter(!is.na(.data$target)))

      valid_L3 <- data_merge %>%
        dplyr::distinct(.data$group, .data$L3) %>%
        dplyr::count(.data$L3) %>%
        dplyr::filter(n == 2) %>%
        dplyr::pull(.data$L3)

      data_merge <- data_merge %>%
        dplyr::filter(.data$L3 %in% valid_L3)

      cor_L2 <- suppressMessages(data_merge %>%
        dplyr::group_by(.data$L3, .data$type) %>%
        dplyr::summarise(
          spearman_r = stats::cor.test(.data$value, .data$target, method = "spearman")$estimate,
          spearman_p = stats::cor.test(.data$value, .data$target, method = "spearman")$p.value
        ) %>%
        dplyr::mutate(L1 = L1_type, .before = 1) %>%
        dplyr::mutate(L2 = L2_type, .before = 2))
      return(cor_L2)
    }) %>% do.call(rbind, .)
    cor_L1
  }) %>% do.call(rbind, .)


  ## 生存分析结果
  print("##### Step4: Execute survival analysis...  #####")
  sur_res <- lapply(c("OS", "DSS", "DFI", "PFI"), function(event) {
    # event = "OS"
    # print(event)
    print(paste0("=== ", event))
    tcga_surv_sub <- load_data("tcga_surv") %>%
      dplyr::select(sample, contains(event)) %>%
      na.omit()
    colnames(tcga_surv_sub) <- c("Sample", "status", "time")

    data_merge <- suppressMessages(dplyr::inner_join(mol_data_df, tcga_surv_sub))


    sur_merge <- lapply(unique(data_merge$type), function(type) {
      # type = "BRCA"
      # print(type)
      sur_dat <- data_merge[data_merge$type == type, ]
      if (length(unique(sur_dat$group)) == 1) {
        return(NA)
      }
      surv_diff <- survival::survdiff(survival::Surv(time, status) ~ group, data = sur_dat)
      p.val <- 1 - stats::pchisq(surv_diff$chisq, length(surv_diff$n) - 1)
      hr <- exp(survival::coxph(survival::Surv(time, status) ~ group, data = sur_dat)$coefficients)
      c(p.val, hr)
    }) %>%
      do.call(rbind, .) %>%
      as.data.frame() %>%
      dplyr::mutate(type = unique(data_merge$type), .before = 1) %>%
      # tibble::rownames_to_column("type") %>%
      dplyr::rename("logrank_p" = "V1") %>%
      dplyr::mutate(P.direction = ifelse(.data$groupLower < 1, "Higher Risk", "Lower Risk")) %>%
      dplyr::mutate("event" = event, .before = 1) %>%
      dplyr::select(!.data$groupLower)
    sur_merge
  }) %>% do.call(rbind, .)


  utils::write.csv(phe_res,
    row.names = FALSE,
    file = file.path(out_dir, "molecule_clinical_result.csv")
  )
  utils::write.csv(cor_res,
    row.names = FALSE,
    file = file.path(out_dir, "molecule_correlation_result.csv")
  )
  utils::write.csv(sur_res,
    row.names = FALSE,
    file = file.path(out_dir, "molecule_survival_result.csv")
  )

  # phe_res = read.csv(file.path(tempdir(),"molecule_clinical_result.csv"))
  # sur_res = read.csv(file.path(tempdir(),"molecule_survival_result.csv"))
  # cor_res = read.csv(file.path(tempdir(),"molecule_correlation_result.csv"))




  if (out_report) {
    params <- list(
      "sur_res" = sur_res,
      "cor_res" = cor_res,
      "phe_res" = phe_res,
      "id_name" = molecule,
      "id_type" = data_type
    )
    print("##### Step5: Render html report...  #####")
    file.copy(
      from = system.file("rmd", "report_template.Rmd", package = "UCSCXenaShiny"),
      to = file.path(out_dir, "report_template.Rmd"), overwrite = TRUE
    )
    rmarkdown::render(file.path(out_dir, "report_template.Rmd"), "html_document",
      file.path(out_dir, "report_result.html"),
      params = params,
      envir = new.env(parent = globalenv())
    )
    file.remove(file.path(out_dir, "report_template.Rmd"))
  }
  print("All analysis is done!")

  return(list(sur_res = sur_res, cor_res = cor_res, phe_res = phe_res))
}
