
rm(list = ls())
library(tidyverse)
library(maftools)
library(tidyverse)
library(ggpubr)
library(egg)

## 从maf中提取对应的"Tumor_Sample_Barcode"
maf1 <- read.maf("TCGA.LGG.somatic.maf")
maf2 <- read.maf("TCGA.GBM.somatic.maf")
maf <- merge_mafs(list(maf1, maf2))

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()


## diff
Barcode <- maf@variants.per.sample[["Tumor_Sample_Barcode"]]
ID <- str_sub(Barcode, 1, 16)
data <- data.frame(ID, Barcode)

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

surv = dplyr::filter(dat, Cohort %in% c("TCGA_LGG", "TCGA_GBM"))
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"))

cli <- inner_join(surv, data, by="ID")


table(cli$GRSTME)
cli.low <- subset(cli, GRSTME=="GRSLow/TMEHigh")$Barcode
cli.high <- subset(cli, GRSTME=="GRSHigh/TMELow")$Barcode

# 使用subsetMaf构建不同分组的MAF对象
maf.low <- subsetMaf(maf=maf, tsb=cli.low, isTCGA=FALSE)
maf.high <- subsetMaf(maf=maf, tsb=cli.high, isTCGA=FALSE)

# 瀑布图
oncoplot(maf=maf.low, top=15, borderCol=NULL)
oncoplot(maf=maf.high, top=15, borderCol=NULL)

# 使用mafCompare比较差异突变基因
fvsm <- mafCompare(m1=maf.low, m2=maf.high, m1Name="GRSlow/TMEhigh", m2Name="GRShigh/TMElow", minMut=5)
forestPlot(mafCompareRes=fvsm, pVal=0.000518436, color=c("maroon", "royalblue"), geneFontSize=0.8)
write.table(fvsm$results, file="low_vs_high.tsv", quote=FALSE, row.names=FALSE, sep="\t")


## TMB
tmb <- tmb(maf = maf, captureSize = 50, logScale = TRUE)
colnames(tmb)[c(1,3)] = c("ID", "TMB")
tmb$ID = str_sub(tmb$ID, 1, 16)
write.table(tmb, file="glomia_TMB.csv", sep=",", row.names=F, quote=F)

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


library(egg)
library(ggpubr)
p2 <- ggplot(data, aes(x=GRSTME, y=TMB)) +
  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("#00BA38", "#619CFF", "#F8766D")) +
  scale_y_continuous(limits = c(0, 3))

p2

##
exp = read.table("Meta_TCGA.exp.txt", header = T, sep = "\t", check.names = F, row.names = 1)
table(surv$GRSTME)
ID = intersect(colnames(exp), surv$ID[surv$GRSTME %in% c("GRSLow/TMEHigh", "GRSHigh/TMELow")])

exp = exp[,ID]
surv = tibble::column_to_rownames(surv, "ID")
cli1 = surv[ID,]

p = identical(rownames(cli1),colnames(exp));p
if(!p) exp = exp[,match(rownames(cli1),colnames(exp))]

group = cli1$GRSTME
Group = factor(group,
               levels = c("GRSLow/TMEHigh","GRSHigh/TMELow"))

library(limma)
design=model.matrix(~Group)
fit=lmFit(exp, design)
fit=eBayes(fit)
deg=topTable(fit, coef=2, number = Inf)

library(dplyr)
deg <- mutate(deg, ID=rownames(deg))
write.csv(deg, "diff_gene.csv", row.names = F)


BiocManager::install("clusterProfiler")
library("clusterProfiler")
library("org.Hs.eg.db")
library("enrichplot")
library("ggplot2")

rt = read.csv("diff_gene.csv")
rt = rt[rt$adj.P.Val < 0.05,]

## ID转换
gene <- bitr(rt$ID, 
             fromType="SYMBOL", toType="ENTREZID", 
             OrgDb="org.Hs.eg.db")

## 去重合并
gene <- dplyr::distinct(gene, SYMBOL, .keep_all=TRUE)
gene_df <- data.frame(logFC=rt$logFC, SYMBOL = rt$ID)
gene_df <- merge(gene_df, gene, by="SYMBOL")

## 排序很重要
geneList <- gene_df$logFC
names(geneList) = gene_df$ENTREZID
geneList = sort(geneList, decreasing = TRUE)

## KEGG GSEA
R.utils::setOption( "clusterProfiler.download.method",'auto' )
kk2 <- gseKEGG(geneList     = geneList,
               organism     = 'hsa',
               minGSSize    = 10,
               pvalueCutoff = 0.05,
               verbose      = FALSE)
kk = as.data.frame(kk2)
write.csv(kk2, file="gseaKEGG.csv", row.names = F)

##
hallmarks = read.gmt("h.all.v2023.1.Hs.entrez.gmt")
kegg = read.gmt("c2.cp.kegg.v2023.1.Hs.entrez.gmt")
go = read.gmt("c5.go.bp.v2023.1.Hs.entrez.gmt")
ms = rbind(hallmarks, kegg)

go$term = str_remove(go$term,"GOBP_")
go$term = tolower(go$term)

gsea <- GSEA(geneList, 
             minGSSize = 10,
             maxGSSize = 500,
             eps = 1e-10,
             pvalueCutoff = 0.05,
             pAdjustMethod = "BH",
             TERM2GENE = hallmarks,
             by = "fgsea")
res <- data.frame(gsea)
write.csv(res, "hallmarks.csv", row.names = F)

##
gseap1 <- gseaplot2(gsea,
                    "KEGG_CYTOKINE_CYTOKINE_RECEPTOR_INTERACTION",#富集的ID编号
                    #title = up_gsea$Description[i],#标题
                    color = "red", #GSEA线条颜色
                    base_size = 20,#基础字体大小
                    rel_heights = c(1.5, 0.5, 1),#副图的相对高度
                    subplots = 1:3,   #要显示哪些副图 如subplots=c(1,3) #只要第一和第三个图
                    ES_geom = "line", #enrichment score用线还是用点"dot"
                    pvalue_table = T) #显示pvalue等信息




