#install.packages("survival")
#install.packages("survminer")

library(survival)
library("survminer")

setwd("C:\\Users\\lexb4\\Desktop\\77CGGA\\10.singleGeneSurvival")                #工作目录（需修改）
rt=read.table("singleGeneSurData.txt",header=T,sep="\t",check.names=F,row.names=1)     #读取输入文件
gene="FCGBP"                                                                   #基因名字

a=ifelse(rt[,gene]<=median(rt[,gene]),"low","high")
diff=survdiff(Surv(futime, fustat) ~a,data = rt)
pValue=1-pchisq(diff$chisq,df=1)
fit=survfit(Surv(futime, fustat) ~ a, data = rt)
if(pValue<0.001){
		      pValue="<0.001"
		  }else{
		      pValue=paste0("=",round(pValue,3))
		  }
surPlot=ggsurvplot(fit, 
			           data=rt,
			           conf.int=TRUE,
			           pval=paste0("p",pValue),
			           pval.size=5,
			           risk.table=T,
			           legend.labs=c("high","low"),
			           legend.title=paste0(gene," level"),
			           xlab="Time(years)",
			           break.time.by = 1,
			           risk.table.title="",
			           palette=c("red", "blue"),
			           risk.table.height=.25)          
pdf(file=paste(gene,".survival.pdf",sep=""), width = 6, height = 5.5,onefile = FALSE)
print(surPlot)
dev.off()

summary(fit)                 #查看五年生存率

###

#install.packages('survival')

library(survival)
setwd("C:\\Users\\lexb4\\Desktop\\77CGGA\\11.uniIndep")                                   #设置工作目录
rt=read.table("singleGeneSurData.txt",header=T,sep="\t",check.names=F,row.names=1)        #读取输入文件
colnames(rt)=gsub("_status","",colnames(rt))

outTab=data.frame()
for(i in colnames(rt[,3:ncol(rt)])){
	 cox <- coxph(Surv(futime, fustat) ~ rt[,i], data = rt)
	 coxSummary = summary(cox)
	 coxP=coxSummary$coefficients[,"Pr(>|z|)"]
	 outTab=rbind(outTab,
	              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|)"])
	              )
}
write.table(outTab,file="uniCox.txt",sep="\t",row.names=F,quote=F)

######绘制森林图######
#读取输入文件
rt <- read.table("uniCox.txt",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="uniForest.pdf", width = 7,height = 4)
n <- nrow(rt)
nRow <- n+1
ylim <- c(1,nRow)
layout(matrix(c(1,2),nc=2),width=c(3,2.2))

#绘制森林图左边的临床信息
xlim = c(0,3)
par(mar=c(4,2,2,1))
plot(1,xlim=xlim,ylim=ylim,type="n",axes=F,xlab="",ylab="")
text.cex=0.75
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,font=2,adj=1)
text(3,n:1,Hazard.ratio,adj=1,cex=text.cex);text(3,n+1,'Hazard ratio',cex=text.cex,font=2,adj=1,)

#绘制森林图
par(mar=c(4,1,2,1),mgp=c(2,0.5,0))
xlim = c(0,max(as.numeric(hrLow),as.numeric(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=1,col="black",lty=2,lwd=2)
boxcolor = ifelse(as.numeric(hr) > 1, 'green', 'green')
points(as.numeric(hr), n:1, pch = 15, col = boxcolor, cex=1.3)
axis(1)
dev.off()

###

library(survival)
setwd("C:\\Users\\lexb4\\Desktop\\77CGGA\\12.multiIndep")                        #设置工作目录
rt=read.table("singleGeneSurData.txt",header=T,sep="\t",check.names=F,row.names=1)
colnames(rt)=gsub("_status","",colnames(rt))

multiCox=coxph(Surv(futime, fustat) ~ ., data = rt)
multiCoxSum=summary(multiCox)

outTab=data.frame()
outTab=cbind(
             HR=multiCoxSum$conf.int[,"exp(coef)"],
             HR.95L=multiCoxSum$conf.int[,"lower .95"],
             HR.95H=multiCoxSum$conf.int[,"upper .95"],
             pvalue=multiCoxSum$coefficients[,"Pr(>|z|)"])
outTab=cbind(id=row.names(outTab),outTab)
write.table(outTab,file="multiCox.txt",sep="\t",row.names=F,quote=F)

######绘制森林图######
#读取输入文件
rt <- read.table("multiCox.txt",header=T,sep="\t",row.names=1,check.names=F)
row.names(rt)=gsub("\\`","",row.names(rt))
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="multiForest.pdf", width = 7,height = 4)
n <- nrow(rt)
nRow <- n+1
ylim <- c(1,nRow)
layout(matrix(c(1,2),nc=2),width=c(3,2.2))

#绘制森林图左边的临床信息
xlim = c(0,3)
par(mar=c(4,2,2,1))
plot(1,xlim=xlim,ylim=ylim,type="n",axes=F,xlab="",ylab="")
text.cex=0.75
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,font=2,adj=1)
text(3,n:1,Hazard.ratio,adj=1,cex=text.cex);text(3,n+1,'Hazard ratio',cex=text.cex,font=2,adj=1,)

#绘制森林图
par(mar=c(4,1,2,1),mgp=c(2,0.5,0))
xlim = c(0,max(as.numeric(hrLow),as.numeric(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=1,col="black",lty=2,lwd=2)
boxcolor = ifelse(as.numeric(hr) > 1, 'red', 'red')
points(as.numeric(hr), n:1, pch = 15, col = boxcolor, cex=1.3)
axis(1)
dev.off()

###

library(survivalROC)
setwd("C:\\Users\\lexb4\\Desktop\\77CGGA\\13.ROC")      #设置工作目录
rt=read.table("singleGeneSurData.txt",header=T,sep="\t",check.names=F,row.names=1)    #读取cox回归风险文件
rocCol=c("red","green","blue")
aucText=c()

#绘制5年的ROC曲线
pdf(file="ROC.pdf",width=6,height=6)
par(oma=c(0.5,1,0,1),font.lab=1.5,font.axis=1.5)
roc=survivalROC(Stime=rt$futime, status=rt$fustat, marker = rt[,3], predict.time =5, method="KM")
plot(roc$FP, roc$TP, type="l", xlim=c(0,1), ylim=c(0,1),col=rocCol[1], 
  xlab="False positive rate", ylab="True positive rate",
  lwd = 2, cex.main=1.3, cex.lab=1.2, cex.axis=1.2, font=1.2)
aucText=c(aucText,paste0("five year"," (AUC=",sprintf("%.3f",roc$AUC),")"))
abline(0,1)

#绘制3年的ROC曲线
roc=survivalROC(Stime=rt$futime, status=rt$fustat, marker = rt[,3], predict.time =3, method="KM")
aucText=c(aucText,paste0("three year"," (AUC=",sprintf("%.3f",roc$AUC),")"))
lines(roc$FP, roc$TP, type="l", xlim=c(0,1), ylim=c(0,1),col=rocCol[2],lwd = 2)

#绘制1年的ROC曲线
roc=survivalROC(Stime=rt$futime, status=rt$fustat, marker = rt[,3], predict.time =1, method="KM")
aucText=c(aucText,paste0("one year"," (AUC=",sprintf("%.3f",roc$AUC),")"))
lines(roc$FP, roc$TP, type="l", xlim=c(0,1), ylim=c(0,1),col=rocCol[3],lwd = 2)

legend("bottomright", aucText,lwd=2,bty="n",col=rocCol)
dev.off()

###

setwd("C:\\Users\\lexb4\\Desktop\\77CGGA\\14.prepareClinicalCor")                            #修改工作目录
expFile="rocSigExp.txt"                                                            #表达数据文件
clinicalFile="clinical.txt"                                                       #临床数据文件
gene="FCGBP"

exp=read.table(expFile,sep="\t",header=T,check.names=F,row.names=1)                #读取表达数据文件
cli=read.table(clinicalFile,sep="\t",header=T,check.names=F,row.names=1)           #读取临床数据文件
samSample=intersect(row.names(exp),row.names(cli))
exp=exp[samSample,]
cli=cli[samSample,]
selectCol=c("futime","fustat",gene)
outTab=cbind(exp[,selectCol],cli)
outTab=cbind(id=row.names(outTab),outTab)
write.table(outTab,file="singleGeneCliData.txt",sep="\t",row.names=F,quote=F)

###

library(beeswarm)
setwd("G:\\77CGGA\\15.singleGeneClinical")                    #修改工作目录
file="singleGeneCliData.txt"                                                       #输入文件
rt=read.table(file,sep="\t",header=T,check.names=F,row.names=1)                    #读取表达数据文件
gene="FCGBP"                                                                   #基因名字

#临床相关性分析，输出图形结果
for(clinical in colnames(rt[,4:ncol(rt)])){
      #定义颜色
	    xlabel=vector()
			tab1=table(rt[,clinical])
			labelNum=length(tab1)
			dotCol=c(2,3)
			if(labelNum==3){
				dotCol=c(2,3,4)
			}
			if(labelNum==4){
				dotCol=c(2,3,4,5)
			}
			if(labelNum>4){
				dotCol=rainbow(labelNum)
			}
			for(i in 1:labelNum){
			  xlabel=c(xlabel,names(tab1[i]) )
			}
	    
	    #相关性检验
	    i=gene
		  rt1=rbind(expression=rt[,i],clinical=rt[,clinical])
		  rt1=as.matrix(t(rt1))
		  rt1 <- as.data.frame(rt1)
		  rt1[,1] <- as.numeric(rt1[,1])
		  if(labelNum==2){
		    rtTest<-wilcox.test(expression ~ clinical, data=rt1)
		  }else{
		    rtTest<-kruskal.test(expression ~ clinical, data = rt1)}
		  pValue=rtTest$p.value
		  pval=0
		  if(pValue<0.001){
			  pval="<0.001"
			}else{
			   pval=paste0("=",sprintf("%.03f",pValue))
		  }

      #可视化
      if(pValue<0.05){
					b = boxplot(expression ~ clinical, data = rt1,outline = FALSE, plot=F)
					yMin=min(b$stats)
					yMax = max(b$stats/5+b$stats)
					n = ncol(b$stats)
					outPdf=paste0(i,".",clinical,".pdf")
					width=ifelse(clinical=="Histology",14,7)
					pdf(file=outPdf,width = width,height = 5)
					par(mar = c(4.5,6,3,3))
					boxplot(expression ~ clinical, data = rt1,names=xlabel,
						     ylab = "Gene expression",main=paste0(i," (p",pval,")"),xlab=clinical,
						     cex.main=1.4, cex.lab=1.4, cex.axis=1.3,ylim=c(yMin,yMax),outline = FALSE)
				  beeswarm(expression ~ clinical, data = rt1, col =dotCol, lwd=0.1,
				         pch = 16, add = TRUE, corral="wrap")
				  dev.off()
		  }
}

###

library(plyr)
library(ggplot2)
library(grid)
library(gridExtra)

setwd("G:\\77CGGA\\18.multipleGSEA")                #设置工作目录
files=grep(".xls",dir(),value=T)                                         #获取目录下的所有xls文件
data = lapply(files,read.delim)                                          #读取每个文件
names(data) = files

dataSet = ldply(data, data.frame)
dataSet$pathway = gsub(".xls","",dataSet$.id)                            #将文件后缀删掉

gseaCol=c("#58CDD9","#7A142C","#5D90BA","#431A3D","#91612D","#6E568C","#E0367A","#D8D155","#64495D","#7CC767","#223D6C","#D20A13","#FFD121","#088247","#11AA4D")
pGsea=ggplot(dataSet,aes(x=RANK.IN.GENE.LIST,y=RUNNING.ES,colour=pathway,group=pathway))+
  geom_line(size = 1.5) + scale_color_manual(values = gseaCol[1:nrow(dataSet)]) +   
  labs(x = "", y = "Enrichment Score", title = "") + scale_x_continuous(expand = c(0, 0)) + 
  scale_y_continuous(expand = c(0, 0),limits = c(min(dataSet$RUNNING.ES - 0.02), max(dataSet$RUNNING.ES + 0.02))) +   
  theme_bw() + theme(panel.grid = element_blank()) + theme(panel.border = element_blank()) + theme(axis.line = element_line(colour = "black")) + theme(axis.line.x = element_blank(),axis.ticks.x = element_blank(),axis.text.x = element_blank()) + 
  geom_hline(yintercept = 0) +   theme(legend.position = c(0,0),legend.justification = c(0,0)) + #legend注释的位值
  guides(colour = guide_legend(title = NULL)) + theme(legend.background = element_blank()) + theme(legend.key = element_blank())+theme(legend.key.size=unit(0.5,'cm'))
pGene=ggplot(dataSet,aes(RANK.IN.GENE.LIST,pathway,colour=pathway))+geom_tile()+
  scale_color_manual(values = gseaCol[1:nrow(dataSet)]) + 
  labs(x = "high expression<----------->low expression", y = "", title = "") + 
  scale_x_discrete(expand = c(0, 0)) + scale_y_discrete(expand = c(0, 0)) +  
  theme_bw() + theme(panel.grid = element_blank()) + theme(panel.border = element_blank()) + theme(axis.line = element_line(colour = "black"))+
  theme(axis.line.y = element_blank(),axis.ticks.y = element_blank(),axis.text.y = element_blank())+ guides(color=FALSE)

gGsea = ggplot_gtable(ggplot_build(pGsea))
gGene = ggplot_gtable(ggplot_build(pGene))
maxWidth = grid::unit.pmax(gGsea$widths, gGene$widths)
gGsea$widths = as.list(maxWidth)
gGene$widths = as.list(maxWidth)
dev.off()

#将图形可视化，保存在"multipleGSEA.pdf"
pdf('multipleGSEA.pdf',      #输出图片的文件
     width=7,                #设置输出图片高度
     height=5.5)             #设置输出图片高度
par(mar=c(5,5,2,5))
grid.arrange(arrangeGrob(gGsea,gGene,nrow=2,heights=c(.8,.3)))
dev.off()

###

gene="FCGBP"                         #基因名称

library(limma)
library(pheatmap)

setwd("C:\\Users\\lexb4\\Desktop\\77CGGA\\20.pheatmap")        #设置工作目录

#读取输入文件
rt=read.table("normalize.txt",sep="\t",header=T,check.names=F)
rt=as.matrix(rt)
rownames(rt)=rt[,1]
exp=rt[,2:ncol(rt)]
dimnames=list(rownames(exp),colnames(exp))
data=matrix(as.numeric(as.matrix(exp)),nrow=nrow(exp),dimnames=dimnames)
data=avereps(data)
rt=data[rowMeans(data)>0.5,]
up=read.table("up.txt",header=F)
down=read.table("down.txt",header=F)

#数据整理
rt=rt[,order(rt[gene,])]
Gene=as.data.frame(t(rt)[,gene])
colnames(Gene)=gene
rt=rt[c(as.vector(up[,1]), as.vector(down[,1])),]
Type=c( rep("postive",nrow(up)), rep("negative",nrow(down)) )
names(Type)=rownames(rt)
Type=as.data.frame(Type)

#设置热图颜色
ann_colors = list(
  Gene = colorRampPalette(c("green", "black", "red"))(50),
  Type = c(postive = "#D95F02", negative = "#1B9E77") )
names(ann_colors)=c(gene,"Type")

#绘制热图
pdf(file="heatmap.pdf",width = 8,height = 5)
pheatmap(rt, 
         annotation_row=Type, 
         annotation_col=Gene, 
         annotation_colors = ann_colors,
         color = colorRampPalette(c("green", "black", "red"))(50),
         fontsize_row=7,
         fontsize_col=5,
         fontsize=7,
         cluster_cols = FALSE,
         cluster_rows = FALSE,
         show_colnames = F)
dev.off()

###

setwd("C:\\Users\\lexb4\\Desktop\\77CGGA\\21.circos")        #设置工作目录
inputFile="normalize.txt"                                  #输入文件
gene="FCGBP"                                           #基因名字

#读取输入文件，对数据进行处理
library(limma)
rt=read.table(inputFile,sep="\t",header=T,check.names=F)
rt=as.matrix(rt)
rownames(rt)=rt[,1]
exp=rt[,2:ncol(rt)]
dimnames=list(rownames(exp),colnames(exp))
data=matrix(as.numeric(as.matrix(exp)),nrow=nrow(exp),dimnames=dimnames)
data=avereps(data)
rt=data[rowMeans(data)>0.5,]

#提取相关基因，计算基因间相关系数
gene=read.table("gene.txt",header=F)
data=t(rt[as.vector(gene[,]),])
cor1=cor(data)

#引用圈图可视化R包
library(circlize)
library(corrplot)
options(stringsAsFactors=F)

#设置图形颜色
pdf("circos.pdf")
col = c(rgb(1,0,0,seq(1,0,length=32)),rgb(0,1,0,seq(0,1,length=32)))
cor1[cor1==1]=0
par(mar=c(2,2,2,4))
c1 = ifelse(c(cor1)>=0,rgb(1,0,0,abs(cor1)),rgb(0,1,0,abs(cor1)))
col1 = matrix(c1,nc=ncol(data))

#绘制圈图
circos.par(gap.degree =c(3,rep(2, nrow(cor1)-1)),start.degree = 180)
chordDiagram(cor1,grid.col=rainbow(ncol(data)),transparency = 0.5,col=col1,symmetric = T)
par(xpd=T)
colorlegend(col, vertical = T,labels=c(1,0,-1),xlim=c(1.1,1.3),ylim=c(-0.4,0.4))       #绘制图例
dev.off()
circos.clear()

###

library(survivalROC)
setwd("C:\\Users\\Administrator\\Desktop\\GEO\\geoSurvical\\ROC")      #设置工作目录
rt=read.table("singleGeneSurData.txt",header=T,sep="\t",check.names=F,row.names=1)    #读取cox回归风险文件
rocCol=c("red","green","blue")
aucText=c()

#绘制5年的ROC曲线
pdf(file="ROC.pdf",width=6,height=6)
par(oma=c(0.5,1,0,1),font.lab=1.5,font.axis=1.5)
roc=survivalROC(Stime=rt$futime, status=rt$fustat, marker = rt[,3], predict.time =5, method="KM")
plot(roc$FP, roc$TP, type="l", xlim=c(0,1), ylim=c(0,1),col=rocCol[1], 
  xlab="False positive rate", ylab="True positive rate",
  lwd = 2, cex.main=1.3, cex.lab=1.2, cex.axis=1.2, font=1.2)
aucText=c(aucText,paste0("five year"," (AUC=",sprintf("%.3f",roc$AUC),")"))
abline(0,1)

#绘制3年的ROC曲线
roc=survivalROC(Stime=rt$futime, status=rt$fustat, marker = rt[,3], predict.time =3, method="KM")
aucText=c(aucText,paste0("three year"," (AUC=",sprintf("%.3f",roc$AUC),")"))
lines(roc$FP, roc$TP, type="l", xlim=c(0,1), ylim=c(0,1),col=rocCol[2],lwd = 2)

#绘制1年的ROC曲线
roc=survivalROC(Stime=rt$futime, status=rt$fustat, marker = rt[,3], predict.time =1, method="KM")
aucText=c(aucText,paste0("one year"," (AUC=",sprintf("%.3f",roc$AUC),")"))
lines(roc$FP, roc$TP, type="l", xlim=c(0,1), ylim=c(0,1),col=rocCol[3],lwd = 2)

legend("bottomright", aucText,lwd=2,bty="n",col=rocCol)
dev.off()

###

library(survival)
library(survminer)

setwd("C:\\Users\\lexb4\\Desktop\\geoSurvical\\11.survivalPlot")                    #工作目录（需修改）
rt=read.table("singleGeneData.txt",header=T,sep="\t",check.names=F,row.names=1)     #读取输入文件
gene=colnames(rt)[3]                                                                #基因名字

a=ifelse(rt[,gene]<=median(rt[,gene]),"low","high")
diff=survdiff(Surv(futime, fustat) ~a,data = rt)
pValue=1-pchisq(diff$chisq,df=1)
fit=survfit(Surv(futime, fustat) ~ a, data = rt)
if(pValue<0.001){
		      pValue="<0.001"
		  }else{
		      pValue=paste0("=",round(pValue,3))
		  }
surPlot=ggsurvplot(fit, 
			       data=rt,
			       conf.int=TRUE,
			       pval=paste0("p",pValue),
			       pval.size=6,
			       risk.table=T,
			       legend.labs=c("high","low"),
			       legend.title=paste0(gene," level"),
			       xlab="Time(years)",
			       break.time.by = 1,
			       risk.table.title="",
			       palette=c("red", "blue"),
			       risk.table.height=.25)          
pdf(file=paste(gene,".survival.pdf",sep=""), width = 6.5, height = 5.5,onefile = FALSE)
print(surPlot)
dev.off()

summary(fit)                 #查看五年生存率

###

library(survival)
setwd("C:\\Users\\lexb4\\Desktop\\geoSurvical\\12.indep")                          #设置工作目录
rt=read.table("singleGeneData.txt",header=T,sep="\t",check.names=F,row.names=1)    #读取输入文件

#单因素独立预后分析
uniTab=data.frame()
for(i in colnames(rt[,3:ncol(rt)])){
	 cox <- coxph(Surv(futime, fustat) ~ rt[,i], data = rt)
	 coxSummary = summary(cox)
	 uniTab=rbind(uniTab,
	              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|)"])
	              )
}
write.table(uniTab,file="uniCox.txt",sep="\t",row.names=F,quote=F)

#多因素独立预后分析
multiCox=coxph(Surv(futime, fustat) ~ ., data = rt)
multiCoxSum=summary(multiCox)
multiTab=data.frame()
multiTab=cbind(
             HR=multiCoxSum$conf.int[,"exp(coef)"],
             HR.95L=multiCoxSum$conf.int[,"lower .95"],
             HR.95H=multiCoxSum$conf.int[,"upper .95"],
             pvalue=multiCoxSum$coefficients[,"Pr(>|z|)"])
multiTab=cbind(id=row.names(multiTab),multiTab)
write.table(multiTab,file="multiCox.txt",sep="\t",row.names=F,quote=F)


############绘制森林图函数############
bioForest=function(coxFile=null,forestCol=null,forestFile=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 = 6.5,height = 5)
		n <- nrow(rt)
		nRow <- n+1
		ylim <- c(1,nRow)
		layout(matrix(c(1,2),nc=2),width=c(3,2.2))
		
		#绘制森林图左边的临床信息
		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,font=2,adj=1)
		text(3,n:1,Hazard.ratio,adj=1,cex=text.cex);text(3,n+1,'Hazard ratio',cex=text.cex,font=2,adj=1,)
		
		#绘制森林图
		par(mar=c(4,1,2,1),mgp=c(2,0.5,0))
		xlim = c(0,max(as.numeric(hrLow),as.numeric(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=1,col="black",lty=2,lwd=2)
		boxcolor = ifelse(as.numeric(hr) > 1, forestCol, forestCol)
		points(as.numeric(hr), n:1, pch = 15, col = boxcolor, cex=1.3)
		axis(1)
		dev.off()
}
############绘制森林图函数############


bioForest(coxFile="uniCox.txt", forestCol="green", forestFile="uniForest.pdf")
bioForest(coxFile="multiCox.txt", forestCol="red", forestFile="multiForest.pdf")

###

library(ggpubr)
setwd("C:\\Users\\lexb4\\Desktop\\geoSurvical\\13.cliCor")                    #修改工作目录
file="singleGeneData.txt"                                                     #输入文件
exp=read.table(file,sep="\t",header=T,check.names=F,row.names=1)              #读取表达数据文件
cli=read.table("clinical.txt",sep="\t",header=T,check.names=F,row.names=1)    #读取临床数据文件

#合并数据
samSample=intersect(row.names(exp),row.names(cli))
exp=exp[samSample,]
cli=cli[samSample,]
rt=cbind(exp[,1:3],cli)
gene=colnames(rt)[3]       #获取基因名字

#临床相关性分析，输出图形结果
for(clinical in colnames(rt[,4:ncol(rt)])){
	data=rt[c(gene,clinical)]
	colnames(data)=c("gene","clinical")
	#设置比较组
	group=levels(factor(data$clinical))
	comp=combn(group,2)
	my_comparisons=list()
    for(i in 1:ncol(comp)){my_comparisons[[i]]<-comp[,i]}
	#绘制boxplot
	boxplot=ggboxplot(data, x="clinical", y="gene", color="clinical",
	          xlab=clinical,
	          ylab=paste(gene,"expression"),
	          legend.title=clinical,
	          add = "jitter")+ 
	stat_compare_means(comparisons = my_comparisons)
	pdf(file=paste0(clinical,".pdf"),width=5.5,height=5)
	print(boxplot)
	dev.off()
}

###

library(survival)
library(survminer)
setwd("D:\\biowolf\\RBP\\24.survival")              #设置工作目录

#绘制生存曲线函数
bioSurvival=function(inputFile=null,outFile=null){
	#读取输入文件
	rt=read.table(inputFile,header=T,sep="\t")
	#比较高低风险组生存差异，得到显著性p值
	diff=survdiff(Surv(futime, fustat) ~risk,data = rt)
	pValue=1-pchisq(diff$chisq,df=1)
	pValue=signif(pValue,4)
	pValue=format(pValue, scientific = TRUE)
	fit <- survfit(Surv(futime, fustat) ~ risk, data = rt)
		
	#绘制生存曲线
	surPlot=ggsurvplot(fit, 
		           data=rt,
		           conf.int=T,
		           pval=paste0("p=",pValue),
		           pval.size=4,
		           risk.table=TRUE,
		           legend.labs=c("High risk", "Low risk"),
		           legend.title="Risk",
		           xlab="Time(years)",
		           break.time.by = 1,
		           risk.table.title="",
		           palette=c("red", "blue"),
		           risk.table.height=.25)
	pdf(file=outFile,onefile = FALSE,width = 6.5,height =5.5)
	print(surPlot)
	dev.off()
}
bioSurvival(inputFile="trainRisk.txt",outFile="trainSurv.pdf")
bioSurvival(inputFile="testRisk.txt",outFile="testSurv.pdf")

###

library(survivalROC)
setwd("D:\\biowolf\\RBP\\25.ROC")             #设置工作目录

rocPlot=function(inputFile=null,outPdf=null){
	rt=read.table(inputFile,header=T,sep="\t",check.names=F,row.names=1)
	pdf(file=outPdf,width=5.5,height=5.5)
	par(oma=c(0.5,1,0,1),font.lab=1.5,font.axis=1.5)
	roc=survivalROC(Stime=rt$futime, status=rt$fustat, marker = rt$riskScore, 
		            predict.time =1, method="KM")
	plot(roc$FP, roc$TP, type="l", xlim=c(0,1), ylim=c(0,1),col='red', 
		 xlab="False positive rate", ylab="True positive rate",
		 main=paste("ROC curve (", "AUC = ",sprintf("%0.3f",roc$AUC),")"),
		 lwd = 2, cex.main=1.3, cex.lab=1.2, cex.axis=1.2, font=1.2)
	abline(0,1)
	dev.off()
}

#绘制train组ROC曲线
rocPlot(inputFile="trainRisk.txt",outPdf="trainROC.pdf")
#绘制test组ROC曲线
rocPlot(inputFile="testRisk.txt",outPdf="testROC.pdf")

###

library(pheatmap)
setwd("D:\\biowolf\\RBP\\26.riskPlot")             #设置工作目录

bioRiskPlot=function(inputFile=null,riskScoreFile=null,survStatFile=null,heatmapFile=null){
	rt=read.table(inputFile,sep="\t",header=T,row.names=1,check.names=F)   #读取输入文件
	rt=rt[order(rt$riskScore),]                                            #按照riskScore对样品排序
		
	#绘制风险曲线
	riskClass=rt[,"risk"]
	lowLength=length(riskClass[riskClass=="low"])
	highLength=length(riskClass[riskClass=="high"])
	line=rt[,"riskScore"]
	line[line>10]=10
	pdf(file=riskScoreFile,width = 10,height = 3.5)
	plot(line, type="p", pch=20,
		 xlab="Patients (increasing risk socre)", ylab="Risk score",
		 col=c(rep("green",lowLength),rep("red",highLength)) )
	abline(h=median(rt$riskScore),v=lowLength,lty=2)
	legend("topleft", c("High risk", "low Risk"),bty="n",pch=19,col=c("red","green"),cex=1.2)
	dev.off()
		
	#绘制生存状态图
	color=as.vector(rt$fustat)
	color[color==1]="red"
	color[color==0]="green"
	pdf(file=survStatFile,width = 10,height = 3.5)
	plot(rt$futime, pch=19,
		 xlab="Patients (increasing risk socre)", ylab="Survival time (years)",
		 col=color)
	legend("topleft", c("Dead", "Alive"),bty="n",pch=19,col=c("red","green"),cex=1.2)
	abline(v=lowLength,lty=2)
	dev.off()
		
	#绘制风险热图
	rt1=rt[c(3:(ncol(rt)-2))]
	rt1=log2(rt1+1)
	rt1=t(rt1)
	annotation=data.frame(type=rt[,ncol(rt)])
	rownames(annotation)=rownames(rt)
	pdf(file=heatmapFile,width = 10,height = 3.5)
	pheatmap(rt1, 
		     annotation=annotation, 
		     cluster_cols = FALSE,
		     fontsize_row=11,
		     show_colnames = F,
		     fontsize_col=3,
		     color = colorRampPalette(c("green", "black", "red"))(50) )
	dev.off()
}
bioRiskPlot(inputFile="trainRisk.txt",riskScoreFile="train.riskScore.pdf",survStatFile="train.survStat.pdf",heatmapFile="train.heatmap.pdf")
bioRiskPlot(inputFile="testRisk.txt",riskScoreFile="test.riskScore.pdf",survStatFile="test.survStat.pdf",heatmapFile="test.heatmap.pdf")

###

library(rms)
setwd("G:\\课件\\29.nomogram")                         #设置工作目录

#TCGA列线图绘制
riskFile="trainRisk.txt"
outFile="train.Nomogram.pdf"
risk=read.table(riskFile,header=T,sep="\t",check.names=F,row.names=1)        #读取风险文件
rt=risk[,1:(ncol(risk)-2)]

#数据打包
dd <- datadist(rt)
options(datadist="dd")
#生成函数
f <- cph(Surv(futime, fustat) ~ FCGBP	+ PRS_type	+ Grade	+ Age	+ Chemo_status	+ IDH_mutation_status	+ p19q_codeletion_status, x=T, y=T, surv=T, data=rt, time.inc=1)
surv <- Survival(f)
#建立nomogram
nom <- nomogram(f, fun=list(function(x) surv(1, x), function(x) surv(3, x), function(x) surv(5, x)), 
    lp=F, funlabel=c("1-year survival", "3-year survival", "5-year survival"), 
    maxscale=100, 
    fun.at=c(0.99, 0.9, 0.8, 0.7, 0.5, 0.3,0.1,0.01))  
#nomogram可视化
pdf(file=outFile,height=6,width=9)
plot(nom)
dev.off()

###

library(survival)
library(survminer)

inputFile="scoreTime.txt"        #输入文件名称
setwd("C:\\Users\\lexb4\\Desktop\\TMEimmune\\09.survival")    #工作目录（需修改）
rt=read.table(inputFile,header=T,sep="\t",check.names=F)      #读取输入文件
rt$futime=rt$futime/365                                       #除以365，单位改成年

outTab=data.frame()
for(score in colnames(rt[,4:ncol(rt)])){
	#根据score中位值对样品分组
	a=ifelse(rt[,score]<=median(rt[,score]),"Low","High")
	#高低组生存差异
	diff=survdiff(Surv(futime, fustat) ~a,data = rt)
	pValue=1-pchisq(diff$chisq,df=1)
	#计算每个时间节点病人数目
	fit=survfit(Surv(futime, fustat) ~ a, data = rt)
    #绘制生存曲线
    scoreKm=cbind(score=score,KM=pValue)
	outTab=rbind(outTab,scoreKm)
	if(pValue<0.001){
		pValue="p<0.001"
	}else{
		pValue=paste0("p=",sprintf("%.03f",pValue))
	}
	titleName=score
	surPlot=ggsurvplot(fit, 
						data=rt,
						conf.int=TRUE,
						pval=pValue,
						pval.size=6,
						risk.table=T,
						#ncensor.plot = TRUE,
						legend.labs=c("high","low"),
						legend.title=titleName,
						xlab="Time(years)",
						break.time.by = 1,
						risk.table.title="",
						palette=c("red", "blue"),
						risk.table.height=.25)          
	pdf(file=paste0("sur.",score,".pdf"), width = 6.5, height = 5.5,onefile = FALSE)
	print(surPlot)
	dev.off()
}
#输出基因和p值表格文件
write.table(outTab,file="survResult.xls",sep="\t",row.names=F,quote=F)

###

options(stringsAsFactors=F)
library(limma)
library(ggpubr)

scoreFile="scores.txt"            #score文件
cliFile="clinical.txt"            #临床数据文件
setwd("C:\\Users\\lexb4\\Desktop\\TMEimmune\\10.cliCor")     #修改工作目录

#读取score文件，并对输入文件整理
rt=read.table(scoreFile,sep="\t",header=T,check.names=F,row.names=1)
data=as.matrix(rt)
rownames(data)=gsub("(.*?)\\-(.*?)\\-(.*?)\\-(.*?)\\-.*","\\1\\-\\2\\-\\3",rownames(data))
data=avereps(data)

#读取临床数据文件
cli=read.table(cliFile,sep="\t",header=T,check.names=F,row.names=1)

#合并数据
samSample=intersect(row.names(data),row.names(cli))
data=data[samSample,]
cli=cli[samSample,]
rt=cbind(data,cli)

#临床相关性分析，输出图形结果
for(clinical in colnames(rt[,4:ncol(rt)])){
	for(score in colnames(rt[,1:3])){
		data=rt[c(score,clinical)]
		colnames(data)=c("score","clinical")
		data=data[(data[,"clinical"]!="unknow"),]
		#设置比较组
		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
		boxplot=ggboxplot(data, x="clinical", y="score", color="clinical",
		          xlab=clinical,
		          ylab=score,
		          legend.title=clinical,
		          add = "jitter")+ 
			stat_compare_means(comparisons = my_comparisons)
		#输出图片
		pdf(file=paste0(score,".",clinical,".pdf"),width=5.5,height=5)
		print(boxplot)
		dev.off()
	}
}

###

#引用包
library("limma")
library("pheatmap")

fdrFilter=0.05                    #fdr临界值
logFCfilter=1                     #logFC临界值
scoreType="StromalScore"          #按StromalScore分组
inputFile="symbol.txt"            #输入文件
scoreFile="scores.txt"            #score文件
setwd("C:\\Users\\lexb4\\Desktop\\TMEimmune\\11.StromalDiff")    #设置工作目录

#读取表达输入文件
rt=read.table(inputFile,sep="\t",header=T,check.names=F)
rt=as.matrix(rt)
rownames(rt)=rt[,1]
exp=rt[,2:ncol(rt)]
dimnames=list(rownames(exp),colnames(exp))
data=matrix(as.numeric(as.matrix(exp)),nrow=nrow(exp),dimnames=dimnames)
data=avereps(data)
data=data[rowMeans(data)>0.1,]


#读取score文件,根据score中位值对样品分组
score=read.table(scoreFile,sep="\t",header=T,check.names=F)
med=median(score[,scoreType])
conTab=score[score[,scoreType]<=med,]
treatTab=score[score[,scoreType]>med,]
con=as.vector(conTab[,1])
treat=as.vector(treatTab[,1])
conNum=length(con)
treatNum=length(treat)
data=cbind(data[,con],data[,treat])

#差异分析
outTab=data.frame()
Type=c(rep(1,conNum),rep(2,treatNum))
for(i in row.names(data)){
	geneName=unlist(strsplit(i,"\\|",))[1]
	geneName=gsub("\\/", "_", geneName)
	rt=rbind(expression=data[i,],Type=Type)
	rt=as.matrix(t(rt))
	wilcoxTest<-wilcox.test(expression ~ Type, data=rt)
	conGeneMeans=mean(data[i,1:conNum])
	treatGeneMeans=mean(data[i,(conNum+1):ncol(data)])
	logFC=log2(treatGeneMeans)-log2(conGeneMeans)  
	pvalue=wilcoxTest$p.value
	conMed=median(data[i,1:conNum])
	treatMed=median(data[i,(conNum+1):ncol(data)])
	diffMed=treatMed-conMed
	if( ((logFC>0) & (diffMed>0)) | ((logFC<0) & (diffMed<0)) ){  
		  outTab=rbind(outTab,cbind(gene=i,conMean=conGeneMeans,treatMean=treatGeneMeans,logFC=logFC,pValue=pvalue))
	 }
}
pValue=outTab[,"pValue"]
fdr=p.adjust(as.numeric(as.vector(pValue)),method="fdr")
outTab=cbind(outTab,fdr=fdr)

#输出所有基因的差异情况
write.table(outTab,file="Stromal.all.xls",sep="\t",row.names=F,quote=F)

#输出差异表格
diffSig=outTab[( abs(as.numeric(as.vector(outTab$logFC)))>logFCfilter & as.numeric(as.vector(outTab$fdr))<fdrFilter),]
write.table(diffSig,file="Stromal.Diff.xls",sep="\t",row.names=F,quote=F)
write.table(diffSig,file="Stromal.Diff.txt",sep="\t",row.names=F,quote=F)

#输出差异基因的表达数据
diffExp=rbind(ID=colnames(data[as.vector(diffSig[,1]),]),data[as.vector(diffSig[,1]),])
write.table(diffExp,file="Stromal.diffExp.txt",sep="\t",col.names=F,quote=F)

#绘制差异基因热图
geneNum=50      #上调和下调绘制热图的基因数目
diffSig=diffSig[order(as.numeric(as.vector(diffSig$logFC))),]
diffGeneName=as.vector(diffSig[,1])
diffLength=length(diffGeneName)
hmGene=c()
if(diffLength>(2*geneNum)){
    hmGene=diffGeneName[c(1:geneNum,(diffLength-geneNum+1):diffLength)]
}else{
    hmGene=diffGeneName
}
hmExp=data[hmGene,]
hmExp=log2(hmExp+0.001)
Type=c(rep("Low",conNum),rep("High",treatNum))
names(Type)=colnames(data)
Type=as.data.frame(Type)
pdf(file="Stromal.heatmap.pdf",height=8,width=10)
pheatmap(hmExp, 
         annotation=Type, 
         color = colorRampPalette(c("blue", "white", "red"))(50),
         cluster_cols =F,
         show_colnames = F,
         fontsize = 8,
         fontsize_row=6,
         fontsize_col=8)
dev.off()

###

#引用包
library("clusterProfiler")
library("org.Hs.eg.db")
library("enrichplot")
library("ggplot2")

pvalueFilter=0.05         #p值过滤条件
qvalueFilter=0.05         #矫正后的p值过滤条件

setwd("C:\\Users\\lexb4\\Desktop\\TMEimmune\\15.GO")         #设置工作目录
rt=read.table("id.txt",sep="\t",header=T,check.names=F)      #读取id.txt文件
rt=rt[is.na(rt[,"entrezID"])==F,]                            #去除基因id为NA的基因
gene=rt$entrezID
geneFC=2^rt$logFC
names(geneFC)=gene

#定义颜色类型
colorSel="qvalue"
if(qvalueFilter>0.05){
	colorSel="pvalue"
}

#GO富集分析
kk=enrichGO(gene = gene,OrgDb = org.Hs.eg.db, pvalueCutoff =1, qvalueCutoff = 1, ont="all", readable =T)
GO=as.data.frame(kk)
GO=GO[(GO$pvalue<pvalueFilter & GO$qvalue<qvalueFilter),]
#保存富集结果
write.table(GO,file="GO.txt",sep="\t",quote=F,row.names = F)

#定义显示Term数目
showNum=10
if(nrow(GO)<30){
	showNum=nrow(GO)
}

#柱状图
pdf(file="barplot.pdf",width = 9,height = 7)
bar=barplot(kk, drop = TRUE, showCategory =showNum,split="ONTOLOGY",color = colorSel) + facet_grid(ONTOLOGY~., scale='free')
print(bar)
dev.off()
		
#气泡图
pdf(file="bubble.pdf",width = 9,height = 7)
bub=dotplot(kk,showCategory = showNum, orderBy = "GeneRatio",split="ONTOLOGY", color = colorSel) + facet_grid(ONTOLOGY~., scale='free')
print(bub)
dev.off()
		
#圈图
pdf(file="circos.pdf",width = 9,height = 6.5)
cnet=cnetplot(kk, foldChange=geneFC, showCategory = 5, circular = TRUE, colorEdge = TRUE)
print(cnet)
dev.off()

###

#引用包
library("clusterProfiler")
library("org.Hs.eg.db")
library("enrichplot")
library("ggplot2")

pvalueFilter=0.05        #p值过滤条件
qvalueFilter=0.05        #矫正后的p值过滤条件

setwd("C:\\Users\\lexb4\\Desktop\\TMEimmune\\16.KEGG")        #设置工作目录
rt=read.table("id.txt",sep="\t",header=T,check.names=F)       #读取id.txt文件
rt=rt[is.na(rt[,"entrezID"])==F,]                             #去除基因id为NA的基因
gene=rt$entrezID
geneFC=2^rt$logFC
names(geneFC)=gene

#定义颜色类型
colorSel="qvalue"
if(qvalueFilter>0.05){
	colorSel="pvalue"
}

#kegg富集分析
kk <- enrichKEGG(gene = gene, organism = "hsa", pvalueCutoff =1, qvalueCutoff =1)
KEGG=as.data.frame(kk)
KEGG$geneID=as.character(sapply(KEGG$geneID,function(x)paste(rt$Gene[match(strsplit(x,"/")[[1]],as.character(rt$entrezID))],collapse="/")))
KEGG=KEGG[(KEGG$pvalue<pvalueFilter & KEGG$qvalue<qvalueFilter),]
#保存富集结果
write.table(KEGG,file="KEGG.txt",sep="\t",quote=F,row.names = F)

#定义显示Term数目
showNum=30
if(nrow(KEGG)<showNum){
	showNum=nrow(KEGG)
}

#柱状图
pdf(file="barplot.pdf",width = 10,height = 7)
barplot(kk, drop = TRUE, showCategory = showNum, color = colorSel)
dev.off()

#气泡图
pdf(file="bubble.pdf",width = 10,height = 7)
dotplot(kk, showCategory = showNum, orderBy = "GeneRatio",color = colorSel)
dev.off()

#圈图
pdf(file="circos.pdf",width = 11,height = 7)
kkx=setReadable(kk, 'org.Hs.eg.db', 'ENTREZID')
cnetplot(kkx, foldChange=geneFC,showCategory = 5, circular = TRUE, colorEdge = TRUE,node_label="all")
dev.off()