
#if (!requireNamespace("BiocManager", quietly = TRUE))
#    install.packages("BiocManager")
#BiocManager::install("limma", version = "3.8")

library("limma")

geneName="FTH1"                                                           #需要提取的基因

library("ggplot2")

#MyPal <- c(brewer.pal(12,"Set3")[c(3,6,8,10,11)],brewer.pal(8,"Set2")[-6])
MyPal <- c("#E41A1C","#377EB8","#4DAF4A","#FF7F00","#7CC767","#223D6C","#D20A13","#FFD121","#088247","#11AA4D")#（"#0571b0"-蓝,"#ca0020"-红）
bioCol=MyPal


#Fig 4###################################
###--- Risk clinical analysis  ---### ---#--------------------------------------------------------------------------------
##########################################
options(stringsAsFactors=F)

library(limma)
library(ggpubr)


##RiskMedian----
inputFilecli="The data for Figure 4-clinical information.txt"             #输入文件  #filename.OS.riskCutoff
rtcli=read.table(inputFilecli,sep="\t",header=T,check.names=F)

filename.OS.riskMedian<- paste("The data for Figure 4-risk and OS.txt",sep = "")
rtrisk=read.table(filename.OS.riskMedian,sep="\t",header=T,check.names=F)
#log2
rtrisk$FTH1=log2(rtrisk$FTH1+1)

rtcli=as.matrix(rtcli)
rtrisk=as.matrix(rtrisk)
rownames(rtcli)=rtcli[,1]
rownames(rtrisk)=rtrisk[,1]
rtcli[duplicated(rownames(rtcli)),]
rtrisk[duplicated(rownames(rtrisk)),]
rtmerge=merge(rtcli,rtrisk,all=FALSE)
duplicated(rtmerge$id)
rtmerge=rtmerge[,-((ncol(rtmerge)-3):(ncol(rtmerge)-2))]
rtmerge=rtmerge[,-(as.numeric(ncol(rtmerge)))]

#risk order
FTH1orderltoh = order(rtmerge[,"FTH1"],decreasing=F)
rtmerge.order=rtmerge[FTH1orderltoh,]
show(rtmerge.order$FTH1)

filename.Tcga.Risk.Acli.merge <- paste("The data for Figure 4-clinical information-FTH1order.txt",sep = "")
write.table(rtmerge.order,file=filename.Tcga.Risk.Acli.merge,sep="\t",quote=F,row.names=F)
rt=read.table(filename.Tcga.Risk.Acli.merge,sep="\t",header=T,check.names=F)

#



###ggpubr
library(plyr)
library(ggplot2)
library(ggpubr)

#设置比较组
for(cli in colnames(rtmerge.order[,2:(ncol(rtmerge.order)-1)])){
  
  #cli = colnames(rt)[11]
  clinical=cli
  gene="FTH1" 
  
  data=rt[c(gene,clinical)]
  bioCol=c("#377EB8","#E41A1C","#FF9900","#6E568C","#7CC767","#223D6C","#D20A13","#FFD121","#088247","#11AA4D")
  bioCol=bioCol[1:length(unique(data[,clinical]))]
  
  #data=data[(data[,gene] > -4),]
  
  colnames(data)=c("gene","clinical")
  data=data[(data[,"clinical"]!="unknow"),]
  data=na.omit(data)
  
  #设置比较组
  group=levels(factor(data$clinical))
  data$clinical=factor(data$clinical, levels=group)
  comp=combn(group,2)
  my_comparisons=list()
  
  for(i in 1:ncol(comp)){my_comparisons[[i]]<-comp[,i]}
  
  
  #绘制箱线图
  boxplot=ggboxplot(data, x="clinical", y="gene", color="clinical",
                    xlab=clinical,
                    ylab="Risk score",
                    legend.title=clinical,
                    palette=bioCol
  )+ 
    stat_compare_means(comparisons=my_comparisons)
  
  
  pdf(file=paste0("FTH1-",clinical,".pdf"),width=4,height=6)
  print(boxplot)
  dev.off()
}

