# Clear the workspace
rm(list = ls())

# Load packages
library(tidyverse)
library(maftools)

# Read data for plotting
drawdata <- read.table("score.group.txt", header = TRUE, sep = "\t")
drawdata$id <- substr(drawdata$id, 1, 12)

# Read MAF files
laml.maf1 <- data.table::fread(paste0("TCGA.", "COAD", ".varscan.gz"), data.table = FALSE)
laml.maf2 <- data.table::fread(paste0("TCGA.", "READ", ".varscan.gz"), data.table = FALSE)
laml.maf <- rbind(laml.maf1, laml.maf2)

# Modify Tumor_Sample_Barcode
laml.maf$Tumor_Sample_Barcode <- substr(laml.maf$Tumor_Sample_Barcode, 1, 12)
laml.maf$Tumor_Sample_Barcode <- str_replace_all(laml.maf$Tumor_Sample_Barcode, "-", ".")

newdata <- drawdata
colnames(newdata)[colnames(newdata) == "id"] <- "ID"
value <- newdata$ID[which(newdata$ID %in% laml.maf$Tumor_Sample_Barcode)]
newdata <- newdata %>%
  dplyr::filter(ID %in% value)

Highexpr <- newdata$ID[which(newdata[,"group"] == "High")]
Lowexpr <- newdata$ID[which(newdata[,"group"] == "Low")]

High.laml.maf <- laml.maf %>%
  dplyr::filter(Tumor_Sample_Barcode %in% Highexpr)
Low.laml.maf <- laml.maf %>%
  dplyr::filter(Tumor_Sample_Barcode %in% Lowexpr)

High.laml <- read.maf(maf = High.laml.maf)
Low.laml <- read.maf(maf = Low.laml.maf)

vc_cols <- RColorBrewer::brewer.pal(n = 8, name = 'Paired')
names(vc_cols) <- c(
  'Frame_Shift_Del',
  'Missense_Mutation',
  'Nonsense_Mutation',
  'Multi_Hit',
  'Frame_Shift_Ins',
  'In_Frame_Ins',
  'Splice_Site',
  'In_Frame_Del'
)

# Create waterfall plots for High and Low groups
pdf(file = paste0("High-score-group-mutation-waterfall-plot.pdf"), width = 8, height = 6)
oncoplot(maf = High.laml, colors = vc_cols, top = 10)
dev.off()

pdf(file = paste0("Low-score-group-mutation-waterfall-plot.pdf"), width = 8, height = 6)
oncoplot(maf = Low.laml, colors = vc_cols, top = 10)
dev.off()

# Perform mutation comparison between High and Low groups
pt.vs.rt <- mafCompare(m1 = High.laml, m2 = Low.laml, m1Name = "High", m2Name = "Low", minMut = 0)
write.csv(pt.vs.rt$results, file = "High-vs-Low-score-group-differential-mutation-gene-comparison.csv")

# Select the top 12 genes for display
pt.vs.rt$results <- pt.vs.rt$results[1:12,]

# Create a forest plot for differential mutation gene comparison
pdf(file = paste0("High-vs-Low-score-group-differential-mutation-gene-comparison.pdf"), width = 8, height = 8)
maftools::forestPlot(mafCompareRes = pt.vs.rt, pVal = 1, geneFontSize = 0.8)
dev.off()