####TMB
tmb=read.table('TMB.txt',header = T,sep = '\t',check.names = F)

load('GPR_TME_combined_group.Rdata')
subtype=data
subtype$Sample=rownames(subtype)
subtype$Sample=stringr::str_sub(subtype$Sample,1,12)
# 查找重复的ID


# 删除重复的行
subtype <- subtype[!duplicated(subtype$Sample), ]

# 将ID作为行名


ss=intersect(subtype$Sample,tmb$id)

rownames(tmb)=tmb$id
rownames(subtype)=NULL
rownames(subtype)=subtype$Sample

tmb=tmb[ss,]
subtype=subtype[ss,]

rt=cbind(tmb,subtype)
library(ggplot2)
library(ggsci)
library(ggpubr)
rt$TMB=log2(rt$TMB+1)
rt$subtype=factor(rt$group,levels = c('GPR_high+TME_low','Mixed','GPR_low+TME_high'))
ggviolin(data=rt,x = 'subtype',y='TMB',fill='subtype',palette=c("#EFC000","#0073C2","#6E568C"))+theme_classic() +
  theme(text = element_text(size=15, colour = "black")) + 
  stat_summary(fun.data = "mean_sdl",  fun.args = list(mult = 1),  geom = "pointrange", color = "black")+stat_compare_means()
library(ggstatsplot)
#install.packages('ggstatsplot')
ggstatsplot::ggbetweenstats(data=rt,x = 'subtype',y='TMB',xlab = 'subtype_TCGA',ylab = 'log2(TMB+1)',ggtheme = theme_bw())+scale_fill_jco()


rt$tmb_group=ifelse(rt$TMB>median(rt$TMB),'TMB_high','TMB_low')
rt=rt[rt$subtype != 'Mixed',]
rt$GPR_TME_TMB=paste0(rt$subtype ,'+',rt$tmb_group)


library(survival)
library(survminer)
# 这是为了统一代码，直接复制
data=rt
data$group=data$GPR_TME_TMB
diff=survdiff(Surv(futime, fustat) ~ group, data = data)
length=length(levels(factor(data[,"group"])))
pValue=1-pchisq(diff$chisq, df=length-1)
if(pValue<0.001){
  pValue="p<0.001"
}else{
  pValue=paste0("p=",sprintf("%.03f",pValue))
}
fit <- survfit(Surv(futime, fustat) ~ group, data = data)
bioCol=c("#EFC000","#0073C2","#6E568C","#7CC767","#223D6C","#D20A13","#FFD121","#088247","#11AA4D")
bioCol=bioCol[1:length]

p=ggsurvplot(fit, 
             data=data,
             conf.int=F,
             pval=pValue,
             pval.size=6,
             legend.title='GPR-TME score',
             legend.labs=levels(factor(data[,"group"])),
             legend = c(0.7, 0.85),
             font.legend=12,
             xlab="Time(Months)",
             palette = bioCol,
             surv.median.line = "hv",
             risk.table=F,
             cumevents=F,
             risk.table.height=.25)
p    


###maf整理#################################################
gc()
load('GPR_TME_combined_group.Rdata')
subtype=data
subtype$Sample=rownames(subtype)

subtype$Subtype=subtype$group
subtype=subtype[subtype$Subtype != 'Mixed',]
subtype=subtype[,c('Sample','Subtype')]

subtype1=subtype[subtype$Subtype=='GPR_high+TME_low',]
subtype2=subtype[subtype$Subtype=='GPR_low+TME_high',]


write.table(subtype,file ='group_for_maf.txt',col.names = NA,quote = F,sep = '\t')



options(stringsAsFactors = F) 
install.packages("maftools")
require(maftools) 
require(dplyr)


### 突变数据的可视化
##加分型
library(data.table)
data=fread('input.maf')
data$Tumor_Sample_Barcode=stringr::str_sub(data$Tumor_Sample_Barcode,1,16)
data1=data[data$Tumor_Sample_Barcode %in% rownames(subtype1),]
data2=data[data$Tumor_Sample_Barcode %in% rownames(subtype2),]

write.table(data1,file='maf_TCGA_hl.maf',row.names = F,sep='\t',quote = F)
write.table(data2,file='maf_TCGA_lh.maf',row.names = F,sep='\t',quote = F)

subtype=read.table('group_for_maf.txt',header = T)
colnames(subtype)[1]='Tumor_Sample_Barcode'

write.table(subtype,'maf_ann.txt',sep = '\t',row.names = F,quote = F)


#颜色
ann_colors=list()
col=c("#efc000","#6E568C")

names(col)=c("GPR_high+TME_low", "GPR_low+TME_high")
ann_colors[["Subtype"]]=col

#低评分组瀑布图
maf=read.maf(maf="maf_TCGA_hl.maf", clinicalData="maf_ann.txt")
oncoplot(maf=maf, clinicalFeatures="Subtype",top = 20, annotationColor=ann_colors, keepGeneOrder=T)

#高评分组瀑布图
maf=read.maf(maf="maf_TCGA_lh.maf", clinicalData="maf_ann.txt")
oncoplot(maf=maf, clinicalFeatures="Subtype", top=20, annotationColor=ann_colors, keepGeneOrder=T)


### 检查点###############################################
##########################################################
gc()
load('KIRC_tpm.Rdata')
checkpoint=read.table('checkpoint.txt',header = F)
rt=exprSet_tcga_mRNA[rownames(exprSet_tcga_mRNA) %in% checkpoint$V1,]

rt=as.data.frame(t(rt))
#读取分组文件
load('GPR_TME_combined_group.Rdata')
subtype=data
subtype$Subtype=factor(subtype$group,levels = c('GPR_high+TME_low','Mixed','GPR_low+TME_high'))

rt=rt[rownames(subtype),]
rt$group=subtype$Subtype

rt=tidyr::pivot_longer(rt,cols=c(-'group'),names_to='gene',values_to='expression')
ggboxplot(
  rt,
  x = "gene",
  y = "expression",
  color = "black",notch = T,
  fill = "group",
  xlab = "Subtype_TCGA",
  ylab = "Expression", palette=c("#EFC000","#0073C2","#6E568C")
) +
  stat_compare_means(
    aes(group = group),
    label = "p.signif", 
    method = "wilcox.test",
    hide.ns = T,
    size = 4.5
  ) +
  theme(axis.text.x = element_text(
    angle =60,
    hjust = 1,
    vjust = 1
  ))

# hla
hla=c('HLA-B','HLA-C','HLA-G','HLA-F','HLA-DOB')
rt=exprSet_tcga_mRNA[hla,]

rt=as.data.frame(t(rt))
#读取分组文件
load('GPR_TME_combined_group.Rdata')
subtype=data
subtype$Subtype=factor(subtype$group,levels = c('GPR_high+TME_low','Mixed','GPR_low+TME_high'))

rt=rt[rownames(subtype),]
rt$group=subtype$Subtype

rt=tidyr::pivot_longer(rt,cols=c(-'group'),names_to='gene',values_to='expression')
ggboxplot(
  rt,
  x = "gene",
  y = "expression",
  color = "black",notch = T,
  fill = "group",
  xlab = "Subtype_TCGA",
  ylab = "Expression", palette=c("#EFC000","#0073C2","#6E568C")
) +
  stat_compare_means(
    aes(group = group),
    label = "p.signif", 
    method = "wilcox.test",
    hide.ns = T,
    size = 4.5
  ) +
  theme(axis.text.x = element_text(
    angle =60,
    hjust = 1,
    vjust = 1
  ))

setwd("G:\\肾细胞癌")
#############################TIDE###############################
load('KIRC_tpm.Rdata')
TIDE <- exprSet_tcga_mRNA#行是基因，列是样本
load('GPR_TME_combined_group.Rdata')
TIDE=TIDE[,rownames(data)]
# 这里为了得到比较好的结果，采用two direction median centered
TIDE <- sweep(TIDE,2, apply(TIDE,2,median,na.rm=T))
TIDE <- sweep(TIDE,1, apply(TIDE,1,median,na.rm=T))
write.table(TIDE,"TIDE_input.txt",sep = "\t",row.names = T,col.names = NA,quote = F)

## 读入结果
tide_rt=read.csv('ABS.CSV.csv',row.names = 1)
# 确认一致
tide_rt=tide_rt[rownames(data),]
identical(rownames(tide_rt),rownames(data))

data$Response=tide_rt$Responder
dev.off()
data$group=factor(data$group,levels = c('GPR_high+TME_low','Mixed','GPR_low+TME_high'))
ggstatsplot::ggbarstats(data = data,x='Response',y='group')

library(ggpubr)
ggboxplot(
  data,
  x = "Response",
  y = 'gpr_score',
  notch = T,
  fill = "Response",add = 'jitter',
  xlab = "ICB response",ylab = 'GPR_score'
  , palette=c("#EFC000","#0073C2","#6E568C")
) +stat_compare_means()

load('TMEscore_and_group.Rdata')

## 读入结果
tide_rt=read.csv('ABS.CSV.csv',row.names = 1)
# 确认一致
tide_rt=tide_rt[rownames(data_immune),]
identical(rownames(tide_rt),rownames(data_immune))

data_immune$Response=tide_rt$Responder
library(ggpubr)
ggboxplot(
  data_immune,
  x = "Response",
  y = 'TME_score',
  notch = T,
  fill = "Response",add = 'jitter',
  xlab = "ICB response",ylab = 'TME_score'
  , palette=c("#EFC000","#0073C2","#6E568C")
) +stat_compare_means()


#### 差异基因的比较(Hl和LH)

data1=data[data$group != 'Mixed',]
rt_hl_lh=exprSet_tcga_mRNA
rt_hl_lh=rt_hl_lh[,rownames(data1)]
group_list=as.factor(as.character(data1$group))
group_list=factor(group_list,levels = c('GPR_high+TME_low','GPR_low+TME_high'))

library(limma)
design=model.matrix(~ group_list)
fit=lmFit(rt_hl_lh,design)
fit=eBayes(fit) 
allDiff=topTable(fit,adjust='fdr',coef=2,number=Inf,p.value=1) 

allDiff_up=allDiff[allDiff$logFC>0.5 & allDiff$adj.P.Val<0.05,]
allDiff_up_gpr_tme=allDiff_up
library(clusterProfiler)
library(org.Hs.eg.db)
library(clusterProfiler)
gene = bitr(rownames(allDiff_up), fromType="SYMBOL", toType="UNIPROT", OrgDb="org.Hs.eg.db")
## 去重
gene <- dplyr::distinct(gene,SYMBOL,.keep_all=TRUE)

gene_df <- data.frame(SYMBOL = rownames(allDiff),logFC=allDiff$logFC)

gene_df <- merge(gene,gene_df,by="SYMBOL")
gene_df <- dplyr::distinct(gene_df,UNIPROT,.keep_all=TRUE)

prote=gene_df[,-1]
write.table(prote,file ='prote_up1.tsv',sep = '\t',row.names = F,col.names = F,quote = F)



allDiff_down=allDiff[allDiff$logFC< -0.5 & allDiff$adj.P.Val<0.05,]

library(clusterProfiler)
library(org.Hs.eg.db)
library(clusterProfiler)
gene = bitr(rownames(allDiff_down), fromType="SYMBOL", toType="UNIPROT", OrgDb="org.Hs.eg.db")
## 去重
gene <- dplyr::distinct(gene,SYMBOL,.keep_all=TRUE)

gene_df <- data.frame(SYMBOL = rownames(allDiff),logFC=allDiff$logFC)

gene_df <- merge(gene,gene_df,by="SYMBOL")
gene_df <- dplyr::distinct(gene_df,UNIPROT,.keep_all=TRUE)

prote=gene_df[,-1]
## 注意负号变绝对值
prote$logFC=abs(prote$logFC)
write.table(prote,file ='prote_down1.tsv',sep = '\t',row.names = F,col.names = F,quote = F)



#### 差异基因的比较(是否免疫治疗有效)
## 只需要接着做
group_list=as.factor(as.character(data1$Response))
group_list=factor(group_list,levels = c('False','True'))

library(limma)
design=model.matrix(~ group_list)
fit=lmFit(rt_hl_lh,design)
fit=eBayes(fit) 
allDiff=topTable(fit,adjust='fdr',coef=2,number=Inf,p.value=1) 

allDiff_up=allDiff[allDiff$logFC>0.5 & allDiff$adj.P.Val<0.05,]
ss=intersect(rownames(allDiff_up_gpr_tme),rownames(allDiff_up))
library(clusterProfiler)
library(org.Hs.eg.db)
library(clusterProfiler)
gene = bitr(rownames(allDiff_up), fromType="SYMBOL", toType="UNIPROT", OrgDb="org.Hs.eg.db")
## 去重
gene <- dplyr::distinct(gene,SYMBOL,.keep_all=TRUE)

gene_df <- data.frame(SYMBOL = rownames(allDiff),logFC=allDiff$logFC)

gene_df <- merge(gene,gene_df,by="SYMBOL")
gene_df <- dplyr::distinct(gene_df,UNIPROT,.keep_all=TRUE)

prote=gene_df[,-1]
write.table(prote,file ='prote_up2.tsv',sep = '\t',row.names = F,col.names = F,quote = F)



allDiff_down=allDiff[allDiff$logFC< -0.5 & allDiff$adj.P.Val<0.05,]

library(clusterProfiler)
library(org.Hs.eg.db)
library(clusterProfiler)
gene = bitr(rownames(allDiff_down), fromType="SYMBOL", toType="UNIPROT", OrgDb="org.Hs.eg.db")
## 去重
gene <- dplyr::distinct(gene,SYMBOL,.keep_all=TRUE)

gene_df <- data.frame(SYMBOL = rownames(allDiff),logFC=allDiff$logFC)

gene_df <- merge(gene,gene_df,by="SYMBOL")
gene_df <- dplyr::distinct(gene_df,UNIPROT,.keep_all=TRUE)

prote=gene_df[,-1]
## 注意负号变绝对值
prote$logFC=abs(prote$logFC)
write.table(prote,file ='prote_down2.tsv',sep = '\t',row.names = F,col.names = F,quote = F)


