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
filename.OS.riskMedian<- paste("The data for Figure 5 and 6-HNSC.OS and DFS.txt",sep = "")
library("limma")
gene=geneName="FTH1"                                                           #需要提取的基因
library("ggplot2")
library(survival)
library(survminer)
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
filename.OS.riskMedian<- paste("The data for Figure 5 and 6-HNSC.OS and DFS.txt",sep = "")
scoreFile=filename.OS.riskMedian     #score
cliFile="The data for Figure 5 and 6-HNSC.clinical.txt"             #临床数据文件
cli=read.table(cliFile, header=T, sep="\t", check.names=F, row.names=1)
sameSample=intersect(row.names(cli), row.names(score))
View(cli)
filename.OS.riskMedian<- paste("The data for Figure 5 and 6-HNSC.OS and DFS.txt",sep = "")
scoreFile=filename.OS.riskMedian     #score
cliFile="The data for Figure 5 and 6-HNSC.clinical.txt"             #临床数据文件
score=read.table(scoreFile, header=T, sep="\t",check.names=F, row.names=1)
cli=read.table(cliFile, header=T, sep="\t", check.names=F, row.names=1)
sameSample=intersect(row.names(cli), row.names(score))
score=score[sameSample,]
cli=cli[sameSample,]
data=cbind(futime=score[,"futime"], fustat=score[,"fustat"], cli, group=score[,"risk"])
head(data)
#提取临床数据
#trait="Age.2"
for(trait in cli.chageName){
rt=data[,c("futime","fustat",trait,"group")]
#rt=rt[(rt[,trait]!="unknow"),]
rt<-na.omit(rt)
colnames(rt)=c("futime", "fustat", "clinical", "group")
#对时间选择10年
rt$fustat[rt$futime>10 ] <- 0
rt$futime[rt$futime>10 ] <- 10
tab=table(rt[,"clinical"])
tab=tab[tab!=0]
#对每个临床信息里面的分类进行循环
for(j in names(tab)){
rt1=rt[(rt[,"clinical"]==j),]
tab1=table(rt1[,"group"])
tab1=tab1[tab1!=0]
labels=names(tab1)
if(length(labels)==2){
titleName=j
if((trait=="age") | (trait=="Age") | (trait=="AGE")){
titleName=paste0("age", j)
}
diff=survdiff(Surv(futime, fustat) ~group,data = rt1)
#设置Pvalue的名字
pValue.subOSName=paste0("pvalue.",trait,".",j,sep = "")
pValue.OS=1-pchisq(diff$chisq, df=1)
if(pValue.OS<0.001){
pval=signif(pValue.OS,4)
pval=format(pval, scientific = TRUE)
}else{
pVal=paste0("p=",sprintf("%.03f", pValue.OS))
}
fit <- survfit(Surv(futime, fustat) ~ group, data = rt1)
bioCol=bioCol[1:length(levels(factor(rt1[,"group"])))]
pValue.OS=round(pValue.OS,3)
#绘制生存曲线
surPlot=ggsurvplot(fit,
data=rt1,
conf.int=F,
pval=pVal,
pval.size=6,
title=paste0("Patients with ",titleName),
legend.title=gene,
legend.labs=labels,
font.legend=12,
xlab="Time(years)",
break.time.by = 1,
palette=bioCol,
risk.table=TRUE,
risk.table.title="",
risk.table.col = "strata",
risk.table.height=.25)
#输出图片
j=gsub(">=","ge",j);j=gsub("<=","le",j);j=gsub(">","gt",j);j=gsub("<","lt",j)
pdf(file=paste0(geneName,"-OS.",trait,".",j,"-Median-Clin.AJCC","-",pValue.OS,".pdf"), onefile = FALSE, width = 5.5, height =5)
print(surPlot)
dev.off()
}
}
}
cli.chageName=c("Age.2","Gender.2","Alcohol.2","Smoking.2","Grade.2nGX",
"T_A.2nTX",
"N_A.2nNX",
"M_A.2nMX",
"Stage_c.2nX",
"Hpv_p16")
#提取临床数据
#trait="Age.2"
for(trait in cli.chageName){
rt=data[,c("futime","fustat",trait,"group")]
#rt=rt[(rt[,trait]!="unknow"),]
rt<-na.omit(rt)
colnames(rt)=c("futime", "fustat", "clinical", "group")
#对时间选择10年
rt$fustat[rt$futime>10 ] <- 0
rt$futime[rt$futime>10 ] <- 10
tab=table(rt[,"clinical"])
tab=tab[tab!=0]
#对每个临床信息里面的分类进行循环
for(j in names(tab)){
rt1=rt[(rt[,"clinical"]==j),]
tab1=table(rt1[,"group"])
tab1=tab1[tab1!=0]
labels=names(tab1)
if(length(labels)==2){
titleName=j
if((trait=="age") | (trait=="Age") | (trait=="AGE")){
titleName=paste0("age", j)
}
diff=survdiff(Surv(futime, fustat) ~group,data = rt1)
#设置Pvalue的名字
pValue.subOSName=paste0("pvalue.",trait,".",j,sep = "")
pValue.OS=1-pchisq(diff$chisq, df=1)
if(pValue.OS<0.001){
pval=signif(pValue.OS,4)
pval=format(pval, scientific = TRUE)
}else{
pVal=paste0("p=",sprintf("%.03f", pValue.OS))
}
fit <- survfit(Surv(futime, fustat) ~ group, data = rt1)
bioCol=bioCol[1:length(levels(factor(rt1[,"group"])))]
pValue.OS=round(pValue.OS,3)
#绘制生存曲线
surPlot=ggsurvplot(fit,
data=rt1,
conf.int=F,
pval=pVal,
pval.size=6,
title=paste0("Patients with ",titleName),
legend.title=gene,
legend.labs=labels,
font.legend=12,
xlab="Time(years)",
break.time.by = 1,
palette=bioCol,
risk.table=TRUE,
risk.table.title="",
risk.table.col = "strata",
risk.table.height=.25)
#输出图片
j=gsub(">=","ge",j);j=gsub("<=","le",j);j=gsub(">","gt",j);j=gsub("<","lt",j)
pdf(file=paste0(geneName,"-OS.",trait,".",j,"-Median-Clin.AJCC","-",pValue.OS,".pdf"), onefile = FALSE, width = 5.5, height =5)
print(surPlot)
dev.off()
}
}
}
scoreFile=filename.OS.riskCutoff     #score分组文件
library("limma")
gene=geneName="FTH1"                                                           #需要提取的基因
library("ggplot2")
library(survival)
library(survminer)
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
filename.OS.riskMedian<- paste("The data for Figure 5 and 6-HNSC.OS and DFS.txt",sep = "")
scoreFile=filename.OS.riskMedian     #score
cliFile="The data for Figure 5 and 6-HNSC.clinical.txt"
score=read.table(scoreFile, header=T, sep="\t",check.names=F, row.names=1)
cli=read.table(cliFile, header=T, sep="\t", check.names=F, row.names=1)
sameSample=intersect(row.names(cli), row.names(score))
score=score[sameSample,]
filename.OS.riskMedian<- paste("The data for Figure 5 and 6-HNSC.OS.txt",sep = "")
scoreFile=filename.OS.riskMedian     #score
cliFile="The data for Figure 5 and 6-HNSC.clinical.txt"
score=read.table(scoreFile, header=T, sep="\t",check.names=F, row.names=1)
cli=read.table(cliFile, header=T, sep="\t", check.names=F, row.names=1)
sameSample=intersect(row.names(cli), row.names(score))
score=score[sameSample,]
cli=cli[sameSample,]
data=cbind(futime=score[,"futime"], fustat=score[,"fustat"], cli, group=score[,"risk"])
head(data)
#AJCC
cli.chageName=c("Age","Gender","Alcohol","Smoking","Grade",
"T",
"N",
"M",
"Stage",
"Hpv")
#提取临床数据
#trait="Age.2"
for(trait in cli.chageName){
rt=data[,c("futime","fustat",trait,"group")]
#rt=rt[(rt[,trait]!="unknow"),]
rt<-na.omit(rt)
colnames(rt)=c("futime", "fustat", "clinical", "group")
#对时间选择10年
rt$fustat[rt$futime>10 ] <- 0
rt$futime[rt$futime>10 ] <- 10
tab=table(rt[,"clinical"])
tab=tab[tab!=0]
#对每个临床信息里面的分类进行循环
for(j in names(tab)){
rt1=rt[(rt[,"clinical"]==j),]
tab1=table(rt1[,"group"])
tab1=tab1[tab1!=0]
labels=names(tab1)
if(length(labels)==2){
titleName=j
if((trait=="age") | (trait=="Age") | (trait=="AGE")){
titleName=paste0("age", j)
}
diff=survdiff(Surv(futime, fustat) ~group,data = rt1)
#设置Pvalue的名字
pValue.subOSName=paste0("pvalue.",trait,".",j,sep = "")
pValue.OS=1-pchisq(diff$chisq, df=1)
if(pValue.OS<0.001){
pval=signif(pValue.OS,4)
pval=format(pval, scientific = TRUE)
}else{
pVal=paste0("p=",sprintf("%.03f", pValue.OS))
}
fit <- survfit(Surv(futime, fustat) ~ group, data = rt1)
bioCol=bioCol[1:length(levels(factor(rt1[,"group"])))]
pValue.OS=round(pValue.OS,3)
#绘制生存曲线
surPlot=ggsurvplot(fit,
data=rt1,
conf.int=F,
pval=pVal,
pval.size=6,
title=paste0("Patients with ",titleName),
legend.title=gene,
legend.labs=labels,
font.legend=12,
xlab="Time(years)",
break.time.by = 1,
palette=bioCol,
risk.table=TRUE,
risk.table.title="",
risk.table.col = "strata",
risk.table.height=.25)
#输出图片
j=gsub(">=","ge",j);j=gsub("<=","le",j);j=gsub(">","gt",j);j=gsub("<","lt",j)
pdf(file=paste0(geneName,"-OS.",trait,".",j,".pdf"), onefile = FALSE, width = 5.5, height =5)
print(surPlot)
dev.off()
}
}
}
#读取输入文件
scoreFile=filename.DFS.riskMedian     #score分组文件
View(cli)
scoreFile="The data for Figure 5 and 6-HNSC.DFS.txt"     #score
score=read.table(scoreFile, header=T, sep="\t",check.names=F, row.names=1)
sameSample=intersect(row.names(cli), row.names(score))
score=score[sameSample,]
cli=cli[sameSample,]
data=cbind(futime=score[,"DFS_futime"], fustat=score[,"DFS_fustat"], cli, group=score[,"risk"])
head(data)
#提取临床数据
#trait="Age.2"
for(trait in cli.chageName){
rt=data[,c("futime","fustat",trait,"group")]
#rt=rt[(rt[,trait]!="unknow"),]
rt<-na.omit(rt)
colnames(rt)=c("futime", "fustat", "clinical", "group")
#对时间选择10年
rt$fustat[rt$futime>10 ] <- 0
rt$futime[rt$futime>10 ] <- 10
tab=table(rt[,"clinical"])
tab=tab[tab!=0]
c_dirSubOSDFS = paste('G:\\生信JM\\2021-All.selected-onegene.HNSCC.assembly\\',geneName,'\\Figure 5 ',geneName,' Subanalyses OSDFS-riskMedian-Clin.AJCC',sep = "")
if (file.exists(c_dirSubOSDFS) == F){
dir.create(c_dirSubOSDFS)
}else { }
setwd(c_dirSubOSDFS)             #设置路径
#对每个临床信息里面的分类进行循环
for(j in names(tab)){
rt1=rt[(rt[,"clinical"]==j),]
tab1=table(rt1[,"group"])
tab1=tab1[tab1!=0]
labels=names(tab1)
if(length(labels)==2){
titleName=j
if((trait=="age") | (trait=="Age") | (trait=="AGE")){
titleName=paste0("age", j)
}
diff=survdiff(Surv(futime, fustat) ~group,data = rt1)
pvalue.subDFSName=paste0("pvalue.",trait,".",j,sep = "")
pvalue.DFS=1-pchisq(diff$chisq, df=1)
if(pvalue.DFS<0.001){
pval=signif(pvalue.DFS,4)
pval=format(pval, scientific = TRUE)
}else{
pval=paste0("p=",sprintf("%.03f", pvalue.DFS))
}
fit <- survfit(Surv(futime, fustat) ~ group, data = rt1)
bioCol=bioCol[1:length(levels(factor(rt1[,"group"])))]
#pvalue.DFS=format(pvalue.DFS, scientific = TRUE)
pvalue.DFS=round(pvalue.DFS,3)
#绘制生存曲线
surPlot=ggsurvplot(fit,
data=rt1,
conf.int=F,
pval=pval,
pval.size=6,
title=paste0("Patients with ",titleName),
legend.title=gene,
legend.labs=labels,
font.legend=12,
xlab="Time(years)",
ylab="DFS probability",
break.time.by = 1,
palette=bioCol,
risk.table=TRUE,
risk.table.title="",
risk.table.col = "strata",
risk.table.height=.25)
#输出图片
j=gsub(">=","ge",j);j=gsub("<=","le",j);j=gsub(">","gt",j);j=gsub("<","lt",j)
pdf(file=paste0(geneName,"-DFS.",trait,".",j,"-Median-Clin.AJCC","-",pvalue.DFS,".pdf"), onefile = FALSE, width = 5.5, height =5)
print(surPlot)
dev.off()
}
}
}
scoreFile="The data for Figure 5 and 6-HNSC.DFS.txt"     #score
score=read.table(scoreFile, header=T, sep="\t",check.names=F, row.names=1)
sameSample=intersect(row.names(cli), row.names(score))
score=score[sameSample,]
cli=cli[sameSample,]
data=cbind(futime=score[,"DFS_futime"], fustat=score[,"DFS_fustat"], cli, group=score[,"risk"])
head(data)
#提取临床数据
#trait="Age.2"
for(trait in cli.chageName){
rt=data[,c("futime","fustat",trait,"group")]
#rt=rt[(rt[,trait]!="unknow"),]
rt<-na.omit(rt)
colnames(rt)=c("futime", "fustat", "clinical", "group")
#对时间选择10年
rt$fustat[rt$futime>10 ] <- 0
rt$futime[rt$futime>10 ] <- 10
tab=table(rt[,"clinical"])
tab=tab[tab!=0]
c_dirSubOSDFS = paste('G:\\生信JM\\2021-All.selected-onegene.HNSCC.assembly\\',geneName,'\\Figure 5 ',geneName,' Subanalyses OSDFS-riskMedian-Clin.AJCC',sep = "")
if (file.exists(c_dirSubOSDFS) == F){
dir.create(c_dirSubOSDFS)
}else { }
setwd(c_dirSubOSDFS)             #设置路径
#对每个临床信息里面的分类进行循环
for(j in names(tab)){
rt1=rt[(rt[,"clinical"]==j),]
tab1=table(rt1[,"group"])
tab1=tab1[tab1!=0]
labels=names(tab1)
if(length(labels)==2){
titleName=j
if((trait=="age") | (trait=="Age") | (trait=="AGE")){
titleName=paste0("age", j)
}
diff=survdiff(Surv(futime, fustat) ~group,data = rt1)
pvalue.subDFSName=paste0("pvalue.",trait,".",j,sep = "")
pvalue.DFS=1-pchisq(diff$chisq, df=1)
if(pvalue.DFS<0.001){
pval=signif(pvalue.DFS,4)
pval=format(pval, scientific = TRUE)
}else{
pval=paste0("p=",sprintf("%.03f", pvalue.DFS))
}
fit <- survfit(Surv(futime, fustat) ~ group, data = rt1)
bioCol=bioCol[1:length(levels(factor(rt1[,"group"])))]
#pvalue.DFS=format(pvalue.DFS, scientific = TRUE)
pvalue.DFS=round(pvalue.DFS,3)
#绘制生存曲线
surPlot=ggsurvplot(fit,
data=rt1,
conf.int=F,
pval=pval,
pval.size=6,
title=paste0("Patients with ",titleName),
legend.title=gene,
legend.labs=labels,
font.legend=12,
xlab="Time(years)",
ylab="DFS probability",
break.time.by = 1,
palette=bioCol,
risk.table=TRUE,
risk.table.title="",
risk.table.col = "strata",
risk.table.height=.25)
#输出图片
j=gsub(">=","ge",j);j=gsub("<=","le",j);j=gsub(">","gt",j);j=gsub("<","lt",j)
pdf(file=paste0(geneName,"-DFS.",trait,".",j,pvalue.DFS,".pdf"), onefile = FALSE, width = 5.5, height =5)
print(surPlot)
dev.off()
}
}
}
#对每个临床信息里面的分类进行循环
for(j in names(tab)){
rt1=rt[(rt[,"clinical"]==j),]
tab1=table(rt1[,"group"])
tab1=tab1[tab1!=0]
labels=names(tab1)
if(length(labels)==2){
titleName=j
if((trait=="age") | (trait=="Age") | (trait=="AGE")){
titleName=paste0("age", j)
}
diff=survdiff(Surv(futime, fustat) ~group,data = rt1)
pvalue.subDFSName=paste0("pvalue.",trait,".",j,sep = "")
pvalue.DFS=1-pchisq(diff$chisq, df=1)
if(pvalue.DFS<0.001){
pval=signif(pvalue.DFS,4)
pval=format(pval, scientific = TRUE)
}else{
pval=paste0("p=",sprintf("%.03f", pvalue.DFS))
}
fit <- survfit(Surv(futime, fustat) ~ group, data = rt1)
bioCol=bioCol[1:length(levels(factor(rt1[,"group"])))]
#pvalue.DFS=format(pvalue.DFS, scientific = TRUE)
pvalue.DFS=round(pvalue.DFS,3)
#绘制生存曲线
surPlot=ggsurvplot(fit,
data=rt1,
conf.int=F,
pval=pval,
pval.size=6,
title=paste0("Patients with ",titleName),
legend.title=gene,
legend.labs=labels,
font.legend=12,
xlab="Time(years)",
ylab="DFS probability",
break.time.by = 1,
palette=bioCol,
risk.table=TRUE,
risk.table.title="",
risk.table.col = "strata",
risk.table.height=.25)
#输出图片
j=gsub(">=","ge",j);j=gsub("<=","le",j);j=gsub(">","gt",j);j=gsub("<","lt",j)
pdf(file=paste0(geneName,"-DFS.",trait,".",j,pvalue.DFS,".pdf"), onefile = FALSE, width = 5.5, height =5)
print(surPlot)
dev.off()
}
}
#提取临床数据
#trait="Age.2"
for(trait in cli.chageName){
rt=data[,c("futime","fustat",trait,"group")]
#rt=rt[(rt[,trait]!="unknow"),]
rt<-na.omit(rt)
colnames(rt)=c("futime", "fustat", "clinical", "group")
#对时间选择10年
rt$fustat[rt$futime>10 ] <- 0
rt$futime[rt$futime>10 ] <- 10
tab=table(rt[,"clinical"])
tab=tab[tab!=0]
#对每个临床信息里面的分类进行循环
for(j in names(tab)){
rt1=rt[(rt[,"clinical"]==j),]
tab1=table(rt1[,"group"])
tab1=tab1[tab1!=0]
labels=names(tab1)
if(length(labels)==2){
titleName=j
if((trait=="age") | (trait=="Age") | (trait=="AGE")){
titleName=paste0("age", j)
}
diff=survdiff(Surv(futime, fustat) ~group,data = rt1)
pvalue.subDFSName=paste0("pvalue.",trait,".",j,sep = "")
pvalue.DFS=1-pchisq(diff$chisq, df=1)
if(pvalue.DFS<0.001){
pval=signif(pvalue.DFS,4)
pval=format(pval, scientific = TRUE)
}else{
pval=paste0("p=",sprintf("%.03f", pvalue.DFS))
}
fit <- survfit(Surv(futime, fustat) ~ group, data = rt1)
bioCol=bioCol[1:length(levels(factor(rt1[,"group"])))]
#pvalue.DFS=format(pvalue.DFS, scientific = TRUE)
pvalue.DFS=round(pvalue.DFS,3)
#绘制生存曲线
surPlot=ggsurvplot(fit,
data=rt1,
conf.int=F,
pval=pval,
pval.size=6,
title=paste0("Patients with ",titleName),
legend.title=gene,
legend.labs=labels,
font.legend=12,
xlab="Time(years)",
ylab="DFS probability",
break.time.by = 1,
palette=bioCol,
risk.table=TRUE,
risk.table.title="",
risk.table.col = "strata",
risk.table.height=.25)
#输出图片
j=gsub(">=","ge",j);j=gsub("<=","le",j);j=gsub(">","gt",j);j=gsub("<","lt",j)
pdf(file=paste0(geneName,"-DFS.",trait,".",j,pvalue.DFS,".pdf"), onefile = FALSE, width = 5.5, height =5)
print(surPlot)
dev.off()
}
}
}
