#if (!require("BiocManager"))
#    install.packages("BiocManager")
BiocManager::install("maftools")

library(maftools)
maf = read.maf(maf = 'TCGA.GBM.somatic.maf')

tmb <- tmb(maf = maf, captureSize = 50, logScale = TRUE)
tmb
write.table(tmb, file="glomia_TMB.csv", sep=",", row.names=F, quote=F)

pdf(file="summary.pdf",width=7,height=6)
plotmafSummary(maf = maf, rmOutlier = TRUE, addStat = 'median', dashboard = TRUE, titvRaw = FALSE)
dev.off()

pdf(file="waterfall.pdf",width=7,height=6)
oncoplot(maf = maf, top = 30, fontSize = 0.5 ,showTumorSampleBarcodes = F )
dev.off()

pdf(file="interaction.pdf",width=7,height=6)
somaticInteractions(maf = maf, top = 25, pvalue = c(0.05, 0.001))
dev.off()

pdf(file="Genecloud.pdf",width=7,height=6)
geneCloud(input = maf, minMut = 5)
dev.off()

maf # 直接输入MAF对象可以查看MAF文件的基本信息
getSampleSummary(maf) # 显示样品的统计
getGeneSummary(maf) # 显示基因的统计
getClinicalData(maf) # 显示样品关联的临床数据
getFields(maf) # 显示MAF文件中的所有字段
write.mafSummary(maf=maf, basename="lihc") # 将详细统计结果输出到文件


library(tidyverse)
library(ggpubr)
library(egg)
tmb = read.csv("glomia_TMB.csv")

hrs = read.csv("0105 Glycolysis.lasson.csv")
tme = read.csv("00202 TME.lasson.csv")
tme = tme[,-c(2,3,4)]
cli = dplyr::inner_join(hrs, tme, by = "ID")

data = dplyr::inner_join(cli, tmb, by="ID")

surv = data
surv$GRS <- ifelse(surv$GRSscore > median(surv$GRSscore),'GRSHigh', 'GRSLow')
surv$TME <- ifelse(surv$TMEscore > median(surv$TMEscore),'TMELow', 'TMEHigh') 

surv$group = paste(surv$GRS, surv$TME, sep = "/")
table(surv$group)

surv$GRSTME = forcats::fct_collapse(surv$group, Mixed = c("GRSHigh/TMEHigh", "GRSLow/TMELow"))
surv$GRSTME = factor(surv$GRSTME, 
                     levels = c("GRSLow/TMEHigh", "Mixed", "GRSHigh/TMELow"),
                     labels = c("GRSLow/TMEHigh", "Mixed", "GRSHigh/TMELow"))



colnames(dat)
my_com = list(c("GRSHigh/TMELow", "Mixed"), 
              c("Mixed", "GRSLow/TMEHigh"), 
              c("GRSHigh/TMELow", "GRSLow/TMEHigh"))



library(egg)
library(ggpubr)

p2 <- ggplot(surv, aes(x=GRSTME, y=total_perMB)) +
  geom_violin(aes(fill=GRSTME)) +
  geom_boxplot(width = 0.2, position = position_dodge(width = 1.0)) +
  theme_article() +
  theme(legend.position = "none",
        axis.text = element_text(colour = "black")) +
  labs(x="", y="Tumor mutation burden (TMB)") +
  stat_compare_means() +
  scale_fill_manual(values = c("#00BFC4", "#7CAE00", "#F8766D")) +
  scale_y_continuous(limits = c(0, 3))

p2




