#PBS -l walltime=1000:00:00
#PBS -l nodes=1:ppn=36
#PBS -S /bin/bash
#PBS -l mem=20gb
#PBS -j oe
#PBS -M w_shixiang@163.com
#PBS -q pub_fast

# Please set PBS arguments above
cd /public/home/wangshx/wangshx/PCAWG-TCGA
module load apps/R/3.6.1

# Following are commands/scripts you want to run
# If you do not know how to set this, please check README.md file 
Rscript call_pcawg_sp.R

library(sigminer)

tally_X <- readRDS("pcawg_cn_tally_X.rds")

sigprofiler_extract(tally_X$nmf_matrix,
  output = "PCAWG_CN176X",
  range = 2:30,
  nrun = 100,
  init_method = "random",
  is_exome = FALSE,
  use_conda = TRUE
)

library(sigminer)
SP_PCAWG <- sigprofiler_import("../SP/PCAWG_CN176X/", order_by_expo = TRUE, type = "all")
saveRDS(SP_PCAWG, file = "../data/pcawg_cn_solutions_sp.rds")

SP_PCAWG <- readRDS("../data/pcawg_cn_solutions_sp.rds")

show_sig_number_survey(
  SP_PCAWG$all_stats %>%
    dplyr::rename(
      s = `Stability (Avg Silhouette)`,
      e = `Mean Cosine Distance`
    ) %>%
    dplyr::mutate(
      SignatureNumber = as.integer(gsub("[^0-9]", "", SignatureNumber))
    ),
  x = "SignatureNumber",
  left_y = "s", right_y = "e",
  left_name = "Stability (Avg Silhouette)",
  right_name = "Mean Cosine Distance",
  highlight = 14
)

pcawg_sigs <- SP_PCAWG$solution_list$S14
apply(pcawg_sigs$Exposure, 1, mean)

sig_names(pcawg_sigs)

colnames(pcawg_sigs$Signature) <- colnames(pcawg_sigs$Signature.norm) <- rownames(pcawg_sigs$Exposure) <- rownames(pcawg_sigs$Exposure.norm) <- pcawg_sigs$Stats$signatures$Signatures <- paste0("CNS", 1:14)

sig_names(pcawg_sigs)

saveRDS(pcawg_sigs, file = "../data/pcawg_cn_sigs_CN176_signature.rds")

sim <- get_sig_similarity(pcawg_sigs, pcawg_sigs)

pheatmap::pheatmap(sim$similarity, cluster_cols = F, cluster_rows = F, display_numbers = TRUE)

df_abs <- get_sig_exposure(pcawg_sigs) %>%
  tidyr::pivot_longer(cols = starts_with("CNS"), names_to = "sig", values_to = "activity") %>%
  dplyr::mutate(
    activity = log10(activity + 1),
    sig = factor(sig, levels = paste0("CNS", 1:14))
  )

df_rel <- get_sig_exposure(pcawg_sigs, type = "relative", rel_threshold = 0) %>%
  tidyr::pivot_longer(cols = starts_with("CNS"), names_to = "sig", values_to = "activity") %>%
  dplyr::mutate(
    sig = factor(sig, levels = paste0("CNS", 1:14))
  )

show_group_distribution(
  df_abs,
  gvar = "sig",
  dvar = "activity",
  order_by_fun = FALSE,
  g_angle = 90,
  ylab = "log10(activity+1)"
)

show_group_distribution(
  df_rel,
  gvar = "sig",
  dvar = "activity",
  order_by_fun = FALSE,
  g_angle = 90,
  ylab = "relative activity"
)

get_sample_from_sig <- function(dt, sig) {
  res <- head(dt[order(dt[[sig]], dt[[paste0("ABS_", sig)]], decreasing = TRUE)], 2L)
  res
}

pcawg_cn_obj <- readRDS("../data/pcawg_cn_obj.rds")
samp_summary <- pcawg_cn_obj@summary.per.sample
rel_activity <- get_sig_exposure(pcawg_sigs, type = "relative", rel_threshold = 0)
abs_activity <- get_sig_exposure(pcawg_sigs, type = "absolute")

rel_activity <- rel_activity[, lapply(.SD, function(x) {
  if (is.numeric(x)) round(x, 2) else x
})]
# all(rownames(abs_activity) == rownames(rel_activity))
colnames(abs_activity) <- c("sample", paste0("ABS_CNS", 1:14))
act <- merge(
  rel_activity, abs_activity,
  by = "sample"
)

dir.create("../output/enrich_samples/", showWarnings = FALSE)
for (i in paste0("CNS", 1:14)) {
  cat(paste0("Most enriched in ", i, "\n"))
  s <- get_sample_from_sig(act, i)
  print(s)
  plist <- show_cn_profile(pcawg_cn_obj,
    samples = s$sample,
    show_title = TRUE,
    return_plotlist = TRUE
  )
  plist <- purrr::map2(plist, s$sample, function(p, x) {
    s <- samp_summary[sample == x]
    text <- paste0(
      "n_of_seg:", s$n_of_seg, "\n",
      "n_of_amp:", s$n_of_amp, "\n",
      "n_of_del:", s$n_of_del, "\n",
      "rel:", act[sample == x][[i]], "\n",
      "abs:", act[sample == x][[paste0("ABS_", i)]], "\n"
    )
    p <- p + annotate("text",
      x = Inf, y = Inf, hjust = 1, vjust = 1,
      label = text, color = "gray50"
    )
    p
  })
  p <- cowplot::plot_grid(plotlist = plist, nrow = 2)
  print(p)
  ggplot2::ggsave(file.path("../output/enrich_samples/", paste0(i, ".pdf")),
    plot = p, width = 12, height = 6
  )
}

sel_samps <- c(
  "SP117655", "SP123894", "SP105708", "SP116525", "SP117309",
  "SP77443", "SP114903", "SP116539", "SP118048", "SP112887",
  "SP135354", "SP37636", "SP43664", "SP96163"
)
# table(pcawg_cn_obj@data[sample == "SP116539"]$chromosome)
# table(pcawg_cn_obj@data[sample == "SP112887"]$chromosome)
plist <- list()
for (s in sel_samps) {
  if (s == sel_samps[4]) {
    p <- show_cn_profile(pcawg_cn_obj,
      samples = s, chrs = c("chr12"),
      show_title = TRUE,
      return_plotlist = TRUE
    )
  } else if (s == sel_samps[8]) {
    p <- show_cn_profile(pcawg_cn_obj,
      samples = s, chrs = c("chr1", "chr2", "chr4", "chr5", "chr6"),
      show_title = TRUE,
      return_plotlist = TRUE
    )
  } else if (s == sel_samps[10]) {
    p <- show_cn_profile(pcawg_cn_obj,
      samples = s, chrs = paste0("chr", c(1:10, 14, 16, 18, 19, 21)),
      show_title = TRUE,
      return_plotlist = TRUE
    )
  } else if (s == sel_samps[13]) {
    p <- show_cn_profile(pcawg_cn_obj,
      samples = s, chrs = c("chr7"),
      show_title = TRUE,
      return_plotlist = TRUE
    )
  } else {
    p <- show_cn_profile(pcawg_cn_obj,
      samples = s,
      show_title = TRUE,
      return_plotlist = TRUE
    )
  }
  plist[[s]] <- p[[1]] + ggplot2::labs(x = NULL)
}

plist <- purrr::pmap(list(plist, sel_samps, paste0("CNS", 1:14)), function(p, x, y) {
  s <- samp_summary[sample == x]
  text <- paste0(
    "n_of_seg:", s$n_of_seg, "\n",
    "n_of_amp:", s$n_of_amp, "\n",
    "n_of_del:", s$n_of_del, "\n",
    "rel:", act[sample == x][[y]], "\n",
    "abs:", act[sample == x][[paste0("ABS_", y)]], "\n"
  )
  p <- p + annotate("text",
    x = Inf, y = Inf, hjust = 1, vjust = 1,
    label = text, color = "gray50"
  )
  p
})
p <- cowplot::plot_grid(plotlist = plist, ncol = 1)
ggplot2::ggsave(file.path("../output/enrich_samples/", "selected_samples.pdf"),
  plot = p, width = 12, height = 42
)

ss <- c("SP102561", "SP123964")

p <- show_cn_profile(pcawg_cn_obj,
  samples = ss,
  show_title = TRUE, ncol = 1
)
ggplot2::ggsave("../output/example_cn_profile.pdf",
  plot = p, width = 12, height = 6
)

tally_X <- readRDS("../data/pcawg_cn_tally_X.rds")

p <- show_catalogue(tally_X, mode = "copynumber", method = "X", style = "cosmic", samples = ss, by_context = TRUE, font_scale = 0.7)
ggplot2::ggsave("../output/example_sig_profile.pdf",
  plot = p, width = 15, height = 4
)

p <- show_sig_exposure(sig_exposure(pcawg_sigs)[, ss], hide_samps = FALSE) + ggplot2::ylab("Activity") + ggpubr::rotate_x_text(0, hjust = 0.5)
ggplot2::ggsave("../output/example_sig_activity.pdf",
  plot = p, width = 4, height = 5
)

show_group_distribution(
  df_abs,
  gvar = "cancer_type",
  dvar = "CNS1",
  order_by_fun = FALSE,
  g_angle = 90,
  point_size = 0.3
)

dir.create("../output/cancer-type-dist", showWarnings = F)
signames <- paste0("CNS", 1:14)
for (i in signames) {
  pxx <- show_group_distribution(df_abs,
    gvar = "cancer_type",
    dvar = i, order_by_fun = FALSE,
    ylab = i,
    g_angle = 90, point_size = 0.3
  )
  ggplot2::ggsave(file.path("../output/cancer-type-dist/", paste0("Absolute_activity_", i, ".pdf")),
    plot = pxx, width = 12, height = 6
  )
  pxx <- show_group_distribution(df_rel,
    gvar = "cancer_type",
    dvar = i, order_by_fun = FALSE,
    ylab = i,
    g_angle = 90, point_size = 0.3
  )
  ggplot2::ggsave(file.path("../output/cancer-type-dist/", paste0("Relative_activity_", i, ".pdf")),
    plot = pxx, width = 12, height = 6
  )
}
rm(pxx)

df <- df_rel %>%
  dplyr::mutate_at(dplyr::vars(dplyr::starts_with("CNS")), ~ ifelse(. > 0.05, 1L, 0L)) %>%
  tidyr::pivot_longer(
    cols = dplyr::starts_with("CNS"),
    names_to = "sig", values_to = "detectable"
  )

df2 <- df_rel %>%
  tidyr::pivot_longer(
    cols = dplyr::starts_with("CNS"),
    names_to = "sig", values_to = "expo"
  )

df3 <- df_abs %>%
  dplyr::mutate_at(dplyr::vars(dplyr::starts_with("CNS")), ~ ifelse(. > 15, 1L, 0L)) %>%
  tidyr::pivot_longer(
    cols = dplyr::starts_with("CNS"),
    names_to = "sig", values_to = "segment_detect"
  )

df <- dplyr::left_join(df, df2,
  by = c("sample", "cancer_type", "sig")
) %>% dplyr::left_join(., df3, by = c("sample", "cancer_type", "sig"))


df_type <- df %>%
  dplyr::group_by(cancer_type, sig) %>%
  dplyr::summarise(
    freq = sum(segment_detect), # directly use count
    expo = median(expo[detectable == 1]),
    n = n(),
    label = paste0(unique(cancer_type), " (n=", n, ")"),
    .groups = "drop"
  ) %>%
  dplyr::group_by(cancer_type) %>%
  dplyr::mutate(pro = freq / sum(freq))

df_type$expo <- ifelse(df_type$freq == 0, 0, df_type$expo)

mps <- unique(df_type[, c("cancer_type", "label")])
mpss <- mps$label
names(mpss) <- mps$cancer_type

summary(df_type$freq)

library(cowplot)

p <- ggplot(
  df_type,
  aes(x = cancer_type, y = factor(sig, levels = paste0("CNS", 1:14)))
) +
  geom_point(aes(size = pro, color = expo)) +
  theme_cowplot() +
  ggpubr::rotate_x_text(60) +
  scale_x_discrete(breaks = mps$cancer_type, labels = mps$label) +
  scale_size_continuous(
    limits = c(0.1, 1),
    breaks = c(0, 0.25, 0.5, 0.75, 1)
  ) +
  scale_color_stepsn(
    colors = viridis::viridis(5, direction = -1),
    breaks = c(0, 0.25, 0.5, 0.75, 1)
  ) +
  labs(
    x = NULL, y = "Copy number signatures",
    color = "Median activity\ndue to signature",
    size = "Proportion of tumors \nwith the signatures"
  )
p

pcawg_cns <- readRDS("../data/pcawg_cn_sigs_CN176_signature.rds") %>%
  .[["Exposure.norm"]] %>%
  t() %>%
  as.data.frame() %>%
  tibble::rownames_to_column(., var = "sample")
pcawg_sbs <- read_csv("../data/PCAWG/PCAWG_sigProfiler_SBS_signatures_in_samples.csv")
pcawg_dbs <- read_csv("../data/PCAWG/PCAWG_sigProfiler_DBS_signatures_in_samples.csv")
pcawg_id <- read_csv("../data/PCAWG/PCAWG_SigProfiler_ID_signatures_in_samples.csv")

# Use signature relative contribution for analysis
pcawg_sbs <- pcawg_sbs[, -c(1, 3)] %>%
  dplyr::rename(sample = `Sample Names`) %>%
  dplyr::select(-SBS43, -c(SBS45:SBS60))

pcawg_dbs <- pcawg_dbs[, -c(1, 3)] %>%
  dplyr::rename(sample = `Sample Names`)

pcawg_id <- pcawg_id[, -c(1, 3)] %>%
  dplyr::rename(sample = `Sample Names`)

df_cns <- df_rel %>%
  dplyr::mutate_at(dplyr::vars(dplyr::starts_with("CNS")), ~ ifelse(. > 0.05, 1L, 0L)) %>%
  tidyr::pivot_longer(
    cols = dplyr::starts_with("CNS"),
    names_to = "sig", values_to = "detectable"
  )

df2_cns <- df_rel %>%
  tidyr::pivot_longer(
    cols = dplyr::starts_with("CNS"),
    names_to = "sig", values_to = "expo"
  )

df3_cns <- df_abs %>%
  dplyr::mutate_at(dplyr::vars(dplyr::starts_with("CNS")), ~ ifelse(. > 15, 1L, 0L)) %>%
  tidyr::pivot_longer(
    cols = dplyr::starts_with("CNS"),
    names_to = "sig", values_to = "segment_detect"
  )

df_cns <- dplyr::left_join(df_cns, df2_cns,
  by = c("sample", "cancer_type", "sig")
) %>% dplyr::left_join(., df3_cns, by = c("sample", "cancer_type", "sig"))


df_type_cns <- df_cns %>%
  dplyr::group_by(cancer_type, sig) %>%
  dplyr::summarise(
    freq = sum(segment_detect), # directly use count
    expo = median(expo[detectable == 1]),
    n = n(),
    label = paste0(unique(cancer_type), " (n=", n, ")"),
    .groups = "drop"
  )

mps <- unique(df_type_cns[, c("cancer_type", "label")])
mpss <- mps$label
names(mpss) <- mps$cancer_type

df_detc_cns <- df_cns %>%
  dplyr::group_by(cancer_type, sample) %>%
  dplyr::summarise(
    signumber = sum(segment_detect),
    .groups = "drop"
  )

num_CNS <- df_detc_cns$signumber

# SBS
pcawg_sbs2 <- dplyr::inner_join(pcawg_sbs,
  df_rel[, c("sample", "cancer_type")],
  by = "sample"
)

df_sbs <- pcawg_sbs2 %>%
  dplyr::mutate_at(dplyr::vars(dplyr::starts_with("SBS")), ~ ifelse(. > 0, 1L, 0L)) %>%
  tidyr::pivot_longer(
    cols = dplyr::starts_with("SBS"),
    names_to = "sig", values_to = "detectable"
  )

df2_sbs <- pcawg_sbs2 %>%
  tidyr::pivot_longer(
    cols = dplyr::starts_with("SBS"),
    names_to = "sig", values_to = "expo"
  )

df_sbs <- dplyr::left_join(df_sbs, df2_sbs,
  by = c("sample", "cancer_type", "sig")
)

df_type_sbs <- df_sbs %>%
  dplyr::group_by(cancer_type, sig) %>%
  dplyr::summarise(
    freq = sum(detectable), # directly use count
    expo = median(expo[detectable == 1]),
    n = n(),
    label = paste0(unique(cancer_type), " (n=", n, ")"),
    .groups = "drop"
  )


df_detc_sbs <- df_sbs %>%
  dplyr::group_by(cancer_type, sample) %>%
  dplyr::summarise(
    signumber = sum(detectable),
    .groups = "drop"
  )
num_SBS <- df_detc_sbs$signumber

# DBS
pcawg_dbs2 <- dplyr::inner_join(pcawg_dbs,
  df_rel[, c("sample", "cancer_type")],
  by = "sample"
)

df_dbs <- pcawg_dbs2 %>%
  dplyr::mutate_at(dplyr::vars(dplyr::starts_with("DBS")), ~ ifelse(. > 0, 1L, 0L)) %>%
  tidyr::pivot_longer(
    cols = dplyr::starts_with("DBS"),
    names_to = "sig", values_to = "detectable"
  )

df2_dbs <- pcawg_dbs2 %>%
  tidyr::pivot_longer(
    cols = dplyr::starts_with("DBS"),
    names_to = "sig", values_to = "expo"
  )

df_dbs <- dplyr::left_join(df_dbs, df2_dbs,
  by = c("sample", "cancer_type", "sig")
)

df_type_dbs <- df_dbs %>%
  dplyr::group_by(cancer_type, sig) %>%
  dplyr::summarise(
    freq = sum(detectable), # directly use count
    expo = median(expo[detectable == 1]),
    n = n(),
    label = paste0(unique(cancer_type), " (n=", n, ")"),
    .groups = "drop"
  )


df_detc_dbs <- df_dbs %>%
  dplyr::group_by(cancer_type, sample) %>%
  dplyr::summarise(
    signumber = sum(detectable),
    .groups = "drop"
  )

num_DBS <- df_detc_dbs$signumber

# ID
# DBS
pcawg_id2 <- dplyr::inner_join(pcawg_id,
  df_rel[, c("sample", "cancer_type")],
  by = "sample"
)

df_id <- pcawg_id2 %>%
  dplyr::mutate_at(dplyr::vars(dplyr::starts_with("ID")), ~ ifelse(. > 0, 1L, 0L)) %>%
  tidyr::pivot_longer(
    cols = dplyr::starts_with("ID"),
    names_to = "sig", values_to = "detectable"
  )

df2_id <- pcawg_id2 %>%
  tidyr::pivot_longer(
    cols = dplyr::starts_with("ID"),
    names_to = "sig", values_to = "expo"
  )

df_id <- dplyr::left_join(df_id, df2_id,
  by = c("sample", "cancer_type", "sig")
)

df_type_id <- df_id %>%
  dplyr::group_by(cancer_type, sig) %>%
  dplyr::summarise(
    freq = sum(detectable), # directly use count
    expo = median(expo[detectable == 1]),
    n = n(),
    label = paste0(unique(cancer_type), " (n=", n, ")"),
    .groups = "drop"
  )


df_detc_id <- df_id %>%
  dplyr::group_by(cancer_type, sample) %>%
  dplyr::summarise(
    signumber = sum(detectable),
    .groups = "drop"
  )

num_ID <- df_detc_id$signumber

# sum
df_num_all <- dplyr::tibble(
  sig_type = rep(
    c("CNS", "SBS", "DBS", "ID"),
    c(length(num_CNS), length(num_SBS), length(num_DBS), length(num_ID))
  ),
  num = c(num_CNS, num_SBS, num_DBS, num_ID)
)
# saveRDS(df_num_all, file = "/home/tzy/projects/CNX-method/data/df_num_all.rds")

library(ggpubr)
num <- ggboxplot(df_num_all,
  x = "sig_type", y = "num",
  ylab = "Signature number",
  color = "sig_type",
  xlab = FALSE,
  legend = "none",
  palette = "jco",
  width = 0.3,
  outlier.size = 0.05
) +
  scale_color_manual(values = c("#0073C2", "#EFC000", "#CD534C", "#7AA6DC"))
num

cancer_num <- ggplot(data = df_detc_cns, aes(x = cancer_type, y = signumber, fill = cancer_type)) +
  geom_boxplot() +
  coord_flip() +
  labs(x = NULL, y = "Signature number") +
  theme_cowplot() +
  theme(legend.position = "none")
cancer_num

library(cola)
library(tidyverse)

act <- readRDS("../data/pcawg_cn_sigs_CN176_activity.rds")

df <- purrr::reduce(
  list(
    act$absolute,
    act$similarity
  ),
  dplyr::left_join,
  by = "sample"
)

mat <- df %>%
  filter(similarity > 0.75) %>%
  select(sample, starts_with("CNS")) %>%
  column_to_rownames("sample") %>%
  t()
mat_adj <- adjust_matrix(mat)
# 1 rows have been removed with too low variance (sd < 0.05 quantile)

rownames(mat_adj)

# Select suitable parameters ----------------------------------------------

ds <- colSums(mat_adj)
boxplot(ds)

set.seed(123)
select_samps <- sample(names(sort(ds)), 500)

boxplot(ds[select_samps])

rl_samp <- run_all_consensus_partition_methods(mat_adj[, select_samps], top_n = 13, mc.cores = 8, max_k = 10)
cola_report(rl_samp, output_dir = "../output/cola_report/pcawg_sigs_500_sampls", mc.cores = 8)

rm(rl_samp)

# Cluster number 2
knitr::include_graphics("cola-test-figs/tab-collect-stats-from-consensus-partition-list-1-1.png")

# Cluster number 3
knitr::include_graphics("cola-test-figs/tab-collect-stats-from-consensus-partition-list-2-1.png")

# Cluster number 4
knitr::include_graphics("cola-test-figs/tab-collect-stats-from-consensus-partition-list-3-1.png")

# Cluster number 5
knitr::include_graphics("cola-test-figs/tab-collect-stats-from-consensus-partition-list-4-1.png")

# Cluster number 6
knitr::include_graphics("cola-test-figs/tab-collect-stats-from-consensus-partition-list-5-1.png")

# Cluster number 7
knitr::include_graphics("cola-test-figs/tab-collect-stats-from-consensus-partition-list-6-1.png")

# Cluster number 8
knitr::include_graphics("cola-test-figs/tab-collect-stats-from-consensus-partition-list-7-1.png")

# Cluster number 9
knitr::include_graphics("cola-test-figs/tab-collect-stats-from-consensus-partition-list-8-1.png")

# Cluster number 10
knitr::include_graphics("cola-test-figs/tab-collect-stats-from-consensus-partition-list-9-1.png")

final <- run_all_consensus_partition_methods(
  mat_adj,
  top_value_method = "ATC",
  partition_method = "skmeans",
  top_n = 13, mc.cores = 8, max_k = 10
)
cola_report(final, output_dir = "../output/cola_report/pcawg_sigs_all_sampls", mc.cores = 8)

saveRDS(final, file = "../data/pcawg_cola_result.rds")

library(cola)
library(ComplexHeatmap)
library(tidyverse)

cluster_res <- readRDS("../data/pcawg_cola_result.rds")
cluster_res <- cluster_res["ATC", "skmeans"]
cluster_res

select_partition_number(cluster_res)

collect_classes(cluster_res)

pdf(file = "../output/cola_all_plots.pdf", width = 16, height = 10)
collect_plots(cluster_res)
dev.off()

pcawg_clusters <- get_classes(cluster_res, k = 5)
saveRDS(pcawg_clusters, file = "../data/pcawg_clusters.rds")

library(ezcox)
library(tidyverse)

cluster_df <- readRDS("../data/pcawg_clusters.rds") %>%
  as.data.frame() %>%
  tibble::rownames_to_column("sample")

tidy_info <- readRDS("../data/pcawg_sample_tidy_info.rds")

df_os <- tidy_info %>%
  dplyr::filter(keep) %>%
  dplyr::select(sample, donor_vital_status, donor_survival_time) %>%
  purrr::set_names(c("sample", "os", "time")) %>%
  na.omit() %>%
  mutate(
    os = ifelse(os == "deceased", 1, 0),
    time = time / 365
  ) %>%
  left_join(cluster_df, by = "sample")

colnames(df_os)[4] <- "cluster"
# Use the group with minimal CNV level as reference group
df_os$cluster <- paste0("subgroup", df_os$cluster)

nrow(df_os)

show_forest(df_os, covariates = "cluster", status = "os", add_caption = FALSE, merge_models = TRUE)

library(survival)
library(survminer)

sfit <- survfit(Surv(time, os) ~ cluster, data = df_os)
ggsurvplot(sfit,
  pval = TRUE,
  fun = "pct",
  xlab = "Time (in years)",
  # palette = "jco",
  legend.title = ggplot2::element_blank(),
  legend.labs = paste0("subgroup", 1:5),
  break.time.by = 5,
  risk.table = TRUE,
  tables.height = 0.4
)

pcawg_sbs <- read_csv("../data/PCAWG/PCAWG_sigProfiler_SBS_signatures_in_samples.csv")
pcawg_dbs <- read_csv("../data/PCAWG/PCAWG_sigProfiler_DBS_signatures_in_samples.csv")
pcawg_id <- read_csv("../data/PCAWG/PCAWG_SigProfiler_ID_signatures_in_samples.csv")

df_merged <- purrr::reduce(
  list(
    tidy_info,
    pcawg_sbs[, -c(1, 3)] %>% dplyr::rename(sample = `Sample Names`) %>%
      # Drop Possible sequencing artefact associated signatures
      dplyr::select(-SBS43, -c(SBS45:SBS60)),
    pcawg_dbs[, -c(1, 3)] %>% dplyr::rename(sample = `Sample Names`),
    pcawg_id[, -c(1, 3)] %>% dplyr::rename(sample = `Sample Names`),
    cluster_df[, 1:2] %>% purrr::set_names(c("sample", "cluster"))
  ),
  dplyr::left_join,
  by = "sample"
) %>%
  dplyr::filter(!is.na(cluster)) %>%
  mutate(cluster = paste0("subgroup", cluster))

colnames(df_merged) <- gsub("Abs_", "", colnames(df_merged))

df_others <- readRDS("../data/pcawg_tidy_anno_df.rds") %>%
  tibble::rownames_to_column("sample") %>%
  dplyr::left_join(
    cluster_df[, 1:2] %>% purrr::set_names(c("sample", "cluster")),
    by = "sample"
  ) %>%
  dplyr::mutate(
    Tumor_Stage = dplyr::case_when(
      Tumor_Stage == "I" ~ 1,
      Tumor_Stage == "II" ~ 2,
      Tumor_Stage == "III" ~ 3,
      Tumor_Stage == "IV" ~ 4
    ),
    isFemale = Sex == "female",
    isWGD = WGD_Status == "wgd",
    isHRD = HRD_Status == "hrd",
    cluster = paste0("subgroup", cluster)
  ) %>%
  dplyr::select(-c(Sex, WGD_Status, HRD_Status))

nm_others <- setdiff(colnames(df_others), c("sample", "cluster"))

enrich_result_others <- group_enrichment(
  df_others,
  grp_vars = "cluster",
  enrich_vars = nm_others,
  co_method = "wilcox.test"
)

enrich_result_others$enrich_var <- factor(enrich_result_others$enrich_var, nm_others)

p <- show_group_enrichment(
  enrich_result_others,
  fill_by_p_value = TRUE,
  cut_p_value = TRUE,
  return_list = T
)
p <- p$cluster + labs(x = NULL, y = NULL)
p + coord_flip()

repair.select <- c(
  "ATM", "BRCA1", "BRCA2", "CDK12", "CHEK2",
  "FANCA", "FANCC", "FANCD2", "FANCF",
  # "MSH2",
  "MSH3", "MSH4",
  # "MSH5","MSH6",
  "PALB2", "POLE", "RAD51B", "TP53"
)

pcawg_pathogenic <- dplyr::inner_join(
  read.csv2("../data/TableS3_panorama_driver_mutations_ICGC_samples.public.tsv.gz", sep = "\t"),
  read.csv2("../data/pcawg_sample_sheet.tsv", sep = "\t")
  %>%
    dplyr::select(., aliquot_id, icgc_specimen_id, dcc_project_code) %>%
    dplyr::rename(., sample_id = "aliquot_id"),
  by = "sample_id"
) %>%
  dplyr::select(., gene, icgc_specimen_id) %>%
  dplyr::rename(., Sample = "icgc_specimen_id") %>%
  subset(., gene %in% repair.select) %>%
  .[!duplicated(.), ] %>%
  dplyr::mutate(num = 1)

load("../data/pcawg_driver_info.RData")
df_mut <- dplyr::left_join(df_others[, c("sample", "cluster")],
  pcawg_driver_number,
  by = "sample"
)
colnames(df_mut)[3] <- "Driver_variants"
df_mut$Driver_variants <- ifelse(is.na(df_mut$Driver_variants), 0, df_mut$Driver_variants)

driver_gene_num <- pcawg_driver_gene %>%
  dplyr::group_by(sample, ROLE) %>%
  dplyr::summarise(driver_genes = n(), .groups = "drop") %>%
  tidyr::pivot_wider(names_from = "ROLE", values_from = "driver_genes", values_fill = 0)

df_mut <- dplyr::left_join(df_mut, driver_gene_num, by = "sample")
colnames(df_mut)[4:5] <- c("Act_genes", "LoF_genes")
df_mut$Act_genes <- ifelse(is.na(df_mut$Act_genes), 0, df_mut$Act_genes)
df_mut$LoF_genes <- ifelse(is.na(df_mut$LoF_genes), 0, df_mut$LoF_genes)

df_mut2 <- df_mut %>%
  dplyr::left_join(
    pcawg_driver_gene %>%
      tidyr::pivot_wider(names_from = "gene", values_from = "ROLE") %>%
      dplyr::mutate_at(vars(-sample), ~ !is.na(.)),
    by = "sample"
  )

## Filter out driver genes with < 10 mutation
gene_sum <- sort(colSums(df_mut2[, -c(1:5)], na.rm = TRUE), decreasing = TRUE)


# 基因按突变数量排序

df_mut2 <- df_mut2 %>% dplyr::select_at(c(colnames(df_mut2)[1:5], names(gene_sum[gene_sum >= 10])))

saveRDS(df_mut2, file = "../data/pcawg_mut_df.rds")

nm_mut <- setdiff(colnames(df_mut2), c("sample", "cluster"))

enrich_result_mut <- group_enrichment(
  df_mut2,
  grp_vars = "cluster",
  enrich_vars = nm_mut,
  co_method = "wilcox.test"
)

keep_df <- enrich_result_mut[, .(count = sum(p_value < 0.05)), by = enrich_var]
drop_vars <- keep_df[count == 0]$enrich_var

enrich_result_mut2 <- enrich_result_mut[!enrich_var %in% drop_vars]

enrich_result_mut2$enrich_var <- factor(enrich_result_mut2$enrich_var, setdiff(nm_mut, drop_vars))

keep_vars <- setdiff(nm_mut, drop_vars)[-c(1:3)]
labels <- c(nm_mut[1:3], paste0(keep_vars, " (", gene_sum[keep_vars], ")"))

p <- show_group_enrichment(
  enrich_result_mut2,
  fill_by_p_value = TRUE,
  cut_p_value = TRUE,
  return_list = T
)
p <- p$cluster + labs(x = NULL, y = NULL)
p + coord_flip() + scale_x_discrete(labels = labels)

rm(list = ls())