#GO analysis
ego <- enrichGO(gene = DEG$ENTREZID,
                OrgDb = org.Hs.eg.db, 
                ont = "all",
                pAdjustMethod = "BH",
                minGSSize = 1,
                pvalueCutoff =0.05, 
                qvalueCutoff =0.05,
                readable = TRUE)

ego_res <- ego@result
save(ego,ego_res,file = "GO_SPP1_DEG.Rdata")

##Air bubble diagram
dotplot(ego, showCategory = 20)
##Category Display
barplot(ego, drop = TRUE, showCategory =10,split="ONTOLOGY") + 
  facet_grid(ONTOLOGY~., scale='free')
dotplot(ego,showCategory = 10,split="ONTOLOGY") + 
  facet_grid(ONTOLOGY~., scale='free')

####KEGG analysis####
setwd("1_KEGG")
#install.packages("tidyverse")
#install.packages("BiocManager")
BiocManager::install('clusterProfiler')
#BiocManager::install('org.Hs.eg.db')
library(tidyverse)
library("BiocManager")
library(org.Hs.eg.db)
library(clusterProfiler)
DEG <- as.data.frame(res)%>% 
  arrange(padj) %>% 
  dplyr::filter(abs(log2FoldChange) > 1, padj < 0.05)

DEG <- DEG %>% rownames_to_column("Gene")

genelist <- bitr(DEG$Gene, fromType="EYA1",
                 toType="ENTREZID", OrgDb='org.Hs.eg.db')
DEG <- inner_join(DEG,genelist,by=c("Gene"="EYA1"))

#KEGG analysis
kk <- enrichKEGG(gene         = DEG$ENTREZID,
                 organism     = 'hsa',
                 pvalueCutoff = 0.1,
                 qvalueCutoff =0.1)
kk_res <- kk@result
save(kk,kk_res,file = "KEGG_1_DEG.Rdata")

load("KEGG_1_DEG.Rdata")

#Histogram
barplot(kk, showCategory = 20,color = "pvalue")
#Air bubble diagram
dotplot(kk, showCategory = 20)

dev.off()

#Plotting the ROC curve
plot(ROC2,
     col="red",   
     xlab="False positive rate", ylab="True positive rate",  
     lty=1,lwd=3,
     main=paste("AUC=",AUC))
abline(0, 1, lty=2, lwd=3)   
dev.off()

####time ROC####
setwd("timeROC")

install.packages("timeROC")
install.packages("survival")
library(timeROC)
library(survival)
library(tidyverse)
#install.packages("stats")


library("stats")           
inputFile="input.txt"      
outFile="pairDiff.pdf"     
ylab="ccRCC expression"      
setwd("D:\\biowolf\\bioR\\14.pairDiff")     

#Read input file
data=read.table(inputFile,sep="\t",header=T,check.names=F,row.names=1)
cond1=colnames(data)[1]
cond2=colnames(data)[2]

#Drawing
pdf(file=outFile, width=5, height=4)
ggpaired(data, cond1 = cond1, cond2 = cond2, fill = "condition", palette = "jco",
         legend.title="Condition",xlab="",ylab = ylab)+
  #stat_compare_means(paired = TRUE, label = "p.format", label.x = 1.35)
  stat_compare_means(paired = TRUE, symnum.args=list(cutpoints = c(0, 0.001, 0.01, 0.05, 1), symbols = c("***", "**", "*", "ns")),label = "p.signif",label.x = 1.35)
dev.off()

library(survival)
library(caret)
library(glmnet)
library(survminer)
library(timeROC)

coxPfilter=0.05        
setwd("C:\\biowolf\\cuproptosis\\14.model")   


rt=read.table("Raw data.txt", header=T, sep="\t", check.names=F, row.names=1)
rt$futime[rt$futime<=0]=1
rt$futime=rt$futime/365
rt[,3:ncol(rt)]=log2(rt[,3:ncol(rt)]+1)


bioForest=function(coxFile=null,forestFile=null,forestCol=null){

  rt <- read.table(coxFile,header=T,sep="\t",row.names=1,check.names=F)
  gene <- rownames(rt)
  hr <- sprintf("%.3f",rt$"HR")
  hrLow  <- sprintf("%.3f",rt$"HR.95L")
  hrHigh <- sprintf("%.3f",rt$"HR.95H")
  Hazard.ratio <- paste0(hr,"(",hrLow,"-",hrHigh,")")
  pVal <- ifelse(rt$pvalue<0.001, "<0.001", sprintf("%.3f", rt$pvalue))
  

  pdf(file=forestFile, width=7, height=6)
  n <- nrow(rt)
  nRow <- n+1
  ylim <- c(1,nRow)
  layout(matrix(c(1,2),nc=2),width=c(3,2.5))
  
 
  xlim = c(0,3)
  par(mar=c(4,2.5,2,1))
  plot(1,xlim=xlim,ylim=ylim,type="n",axes=F,xlab="",ylab="")
  text.cex=0.8
  text(0,n:1,gene,adj=0,cex=text.cex)
  text(1.5-0.5*0.2,n:1,pVal,adj=1,cex=text.cex);text(1.5-0.5*0.2,n+1,'pvalue',cex=text.cex,adj=1)
  text(3,n:1,Hazard.ratio,adj=1,cex=text.cex);text(3,n+1,'Hazard ratio',cex=text.cex,adj=1,)
  

  par(mar=c(4,1,2,1),mgp=c(2,0.5,0))
  LOGindex = 10 
  hrLow = log(as.numeric(hrLow),LOGindex)
  hrHigh = log(as.numeric(hrHigh),LOGindex)
  hr = log(as.numeric(hr),LOGindex)
  xlim = c(floor(min(hrLow,hrHigh)),ceiling(max(hrLow,hrHigh)))
  plot(1,xlim=xlim,ylim=ylim,type="n",axes=F,ylab="",xaxs="i",xlab="Hazard ratio")
  arrows(as.numeric(hrLow),n:1,as.numeric(hrHigh),n:1,angle=90,code=3,length=0.05,col="darkblue",lwd=2.5)
  abline(v=log(1,LOGindex),col="black",lty=2,lwd=2)
  boxcolor = ifelse(as.numeric(hr) > log(1,LOGindex), forestCol[1],forestCol[2])
  points(as.numeric(hr), n:1, pch = 15, col = boxcolor, cex=1.3)
  a1 = axis(1,labels=F,tick=F)
  axis(1,a1,10^a1)
  dev.off()
}



n=1      
for(i in 1:n){

  inTrain<-createDataPartition(y=rt[,2], p=0.5, list=F)
  train<-rt[inTrain,]
  test<-rt[-inTrain,]
  trainOut=cbind(id=row.names(train),train)
  testOut=cbind(id=row.names(test),test)
  

  outUniTab=data.frame()
  sigGenes=c("futime","fustat")
  for(i in colnames(train[,3:ncol(train)])){
    #cox
    cox <- coxph(Surv(futime, fustat) ~ train[,i], data = train)
    coxSummary = summary(cox)
    coxP=coxSummary$coefficients[,"Pr(>|z|)"]
    

    if(coxP<coxPfilter){
      sigGenes=c(sigGenes,i)
      outUniTab=rbind(outUniTab,
                      cbind(id=i,
                            HR=coxSummary$conf.int[,"exp(coef)"],
                            HR.95L=coxSummary$conf.int[,"lower .95"],
                            HR.95H=coxSummary$conf.int[,"upper .95"],
                            pvalue=coxSummary$coefficients[,"Pr(>|z|)"])
      )
    }
  }
  uniSigExp=train[,sigGenes]
  uniSigExpOut=cbind(id=row.names(uniSigExp),uniSigExp)
  if(ncol(uniSigExp)<6){next}