for(i in links){
w = which(link==i)
x1 = xnumsample[w[1]]
y1 = data[w[1]]
x2 = xnumsample[w[2]]
y2 = data[w[2]]
segments(x1,y1,x2,y2)
}
par(xpd=T)
arrows(1,max(data)*1.1,2,max(data)*1.1,angle=90,code=3,length=0.1)
text(1.5,max(data)*1.05,paste("p =",pval),pos=3,cex=1)
dev.off()
#}
#data = df[,1]
#cell = colnames(df)[1]
#cell = gsub(' ','_',cell)
#outpdf=paste(geneName,"-diff.Paired.",plotType,"-h",height,"w",width,"-",pvalue.Pairdot,".pdf",sep="")
#Plot(data,group,outpdf)
#if (!requireNamespace("BiocManager", quietly = TRUE))
#    install.packages("BiocManager")
#BiocManager::install("limma", version = "3.12")
geneName="FTL"                                                           #
#【A】#############################
###  ---different expression  ---###
####################################---#-----------------------------------------------------------------
library(limma)
filename.singlegeneName.txt <- paste("The data for Figure S1A and B.txt",sep = "")
rt=read.table(filename.singlegeneName.txt,sep="\t",header=T,row.names=1,check.names=F)#row.names=F,
#duplicated(rt$ID)
colnames(rt)=c("expression","Type")
rt$expression=log2(rt$expression+1)
ylab = paste(geneName," log2 (expression)",sep="")
##ggplot------
library(plyr)
library(ggplot2)
library(ggpubr)
options(stringsAsFactors=F)
#p-value
clinical= colnames(rt)[2]
gene=colnames(rt)[1]
xlabel=vector()
tab1=table(rt[,clinical])
labelNum=length(tab1)
dotCol=MyPal1
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]) )
}
outTab=data.frame()
rt1=rbind(expression=rt[,gene],clinical=rt[,clinical])
rt1=data.frame(t(rt1))
rt1[,1]=as.numeric(as.character(rt1[,1]))
class(rt1[,1])
if(labelNum==2){
wilcoxTest<-wilcox.test(expression ~ clinical, data=rt1)
}else{
wilcoxTest<-kruskal.test(expression ~ clinical, data = rt1)}
pValue=wilcoxTest$p.value
outTab=rbind(outTab,cbind(gene=gene,pVal=pValue))
pval=0
if(pValue<0.001){
pval=signif(pValue,4)
pval=format(pval, scientific = TRUE)
}else{
pval=round(pValue,3)
}
plotType="dotbPlot"
clinical= colnames(rt)[2]
gene=colnames(rt)[1]
data=rt[c(gene,clinical)]
data=data[(data[,gene] > -4),]
colnames(data)=c("gene","clinical")
data=data[(data[,"clinical"]!="unknow"),]
data=na.omit(data)
#group
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
plotcli<-ggboxplot(data, x="clinical", y="gene", color="clinical",
xlab=clinical,
ylab=ylab,
legend.title=clinical,
add = "jitter")+
stat_compare_means(comparisons = my_comparisons)
width=5
height=6
outFile=paste("Figure S1A.pdf",sep="")
pdf(file=outFile,width=width,height=height)
print(plotcli)
dev.off()
#【B】#############################
###---    paired expression   ---###
####################################---#-------------------------------------------------------------------
MyPal <- c("#E41A1C","#377EB8","#4DAF4A","#FF7F00","#7CC767","#223D6C","#D20A13","#FFD121","#088247","#11AA4D")#（"#0571b0"-蓝,"#ca0020"-红）
rt.singleGene=read.table(filename.singlegeneName.txt,sep="\t",header=T,check.names=F)             #读取文件
colnames(rt.singleGene)=c("ID","expression","Type")
rt.singleGene$expression=log2(rt.singleGene$expression+1)
colnames(rt.singleGene)=c("ID",geneName,"Type")
ylab = paste(geneName," log2 (expression)",sep="")
add.ID=rt.singleGene[,1]
#ID
rt.singleGene.addID=cbind(rt.singleGene,add.ID)
rt.singleGene.addID[,1]=substr(rt.singleGene.addID[,1],1,12)
#normal-tumor
rt.singleGene.addID2 = rt.singleGene.addID
rt.singleGene.Normal=rt.singleGene.addID[grep(pattern="Normal",rt.singleGene.addID[,3]),]
rt.singleGene.Tumor=rt.singleGene.addID[grep(pattern="Tumor",rt.singleGene.addID[,3]),]
rt.singleGene.Normal0=rt.singleGene.Normal
rt.singleGene.Tumor0=rt.singleGene.Tumor
rownames(rt.singleGene.Normal0)=rt.singleGene.Normal0[,1]
rownames(rt.singleGene.Tumor0)=rt.singleGene.Tumor0[,1]
#same name
rt.singleGene.Normal=rt.singleGene.addID[grep(pattern="Normal",rt.singleGene.addID[,3]),]
rownames(rt.singleGene.Normal)=rt.singleGene.Normal[,1]
selename.normal = rownames(rt.singleGene.Normal)
rownames(rt.singleGene.Tumor)=rt.singleGene.Tumor[,1]
rt.singleGene.Tumor=rt.singleGene.Tumor[selename.normal,]
rt.singleGene.Tumor=na.omit(rt.singleGene.Tumor)
rt.singleGene.Pairprepare=merge(rt.singleGene.Tumor,rt.singleGene.Tumor,all=FALSE)
selename.Pair= rownames(rt.singleGene.Pairprepare)=rt.singleGene.Pairprepare[,1]
rt.singleGene.Normal.Pair=rt.singleGene.Normal0[selename.Pair,]
rt.singleGene.Tumor.Pair=rt.singleGene.Tumor0[selename.Pair,]
rownames(rt.singleGene.Normal.Pair)<-NULL
rownames(rt.singleGene.Tumor.Pair)<-NULL
rt.singleGene.NormalTumor.Pair=rbind(rt.singleGene.Normal.Pair,rt.singleGene.Tumor.Pair)
rt.singleGene.NormalTumor.Pair=rt.singleGene.NormalTumor.Pair[,c(4,2,3,1)]
colnames(rt.singleGene.NormalTumor.Pair)[1] ="ID"
rt.singleGene.geneExp=rt.singleGene.NormalTumor.Pair[,-4]
rt.singleGene.Group=rt.singleGene.NormalTumor.Pair[,-2]
#filename.paired.geneName.txt <- paste(geneName,"-","Exp.Paired.singleGene.txt",sep = "")
#write.table(rt.singleGene.geneExp,file=filename.paired.geneName.txt,sep="\t",quote=F,row.names=F)
#filename.paired.Group.txt <- paste(geneName,"-","Paired.samplegroup.txt",sep = "")
#write.table(rt.singleGene.Group,file=filename.paired.Group.txt,sep="\t",quote=F,row.names=F,col.names=F)
#---paired
group = rt.singleGene.Group#
#group = read.table(filename.paired.Group.txt,header=F,sep="\t")                               #读取分组信息文件
df = rt.singleGene.geneExp#
#df = read.table(filename.paired.geneName.txt,row.names=1,header=T,sep="\t",check.names=F)        #读取输入文件
row.names(group) = group[,1]
row.names(df) = df[,1]
df=df[,-1]
m = match(group[,1],rownames(df))
df = df[m,]
#dot----
plotType="dot"
data = df[,1]
cell = colnames(df)[1]
cell = gsub(' ','_',cell)
#Plot = function(data,group,outpdf){
xfactors = as.factor(group[,2])
xnumsample = as.numeric(xfactors)
xaxis = levels(xfactors)
link = group[,3]
links = unique(group[,3])
x1data = data[xnumsample==1]
x2data = data[xnumsample==2]
wilcoxP=wilcox.test(x1data,x2data,exact=FALSE)$p.value#)$p.value
pvalue.Pairdot=signif(wilcoxP,3)
if(pvalue.Pairdot<0.001){
pval=signif(pvalue.Pairdot,4)
pval=format(pvalue.Pairdot, scientific = TRUE)
}else{
pval=round(pvalue.Pairdot,3)
}
width=5
height=6
outpdf=paste("Figure S1B.pdf",sep="")
pdf(outpdf,width=width,height=height)
par(las=1)
yMin=min(data)-max(data)/10                     #纵轴最小值                      #纵轴最小值
yMax=max(data)*1.1                     #纵轴最大值
plot(1,xlim=c(0.5,2.5),ylim=c(yMin,yMax),type="n",xlab="",ylab = ylab,xaxt="n",cex.main=1.5, cex.lab=1.3, cex.axis=1.2)
points(rep(1,length(x1data)),x1data,pch=16,cex=2,col=MyPal[1])
points(rep(2,length(x2data)),x2data,pch=15,cex=2,col=MyPal[2])
axis(1,1:2,xaxis)
for(i in links){
w = which(link==i)
x1 = xnumsample[w[1]]
y1 = data[w[1]]
x2 = xnumsample[w[2]]
y2 = data[w[2]]
segments(x1,y1,x2,y2)
}
par(xpd=T)
arrows(1,max(data)*1.1,2,max(data)*1.1,angle=90,code=3,length=0.1)
text(1.5,max(data)*1.05,paste("p =",pval),pos=3,cex=1)
dev.off()
#}
#data = df[,1]
#cell = colnames(df)[1]
#cell = gsub(' ','_',cell)
#outpdf=paste(geneName,"-diff.Paired.",plotType,"-h",height,"w",width,"-",pvalue.Pairdot,".pdf",sep="")
#Plot(data,group,outpdf)
#if (!requireNamespace("BiocManager", quietly = TRUE))
#    install.packages("BiocManager")
#BiocManager::install("limma", version = "3.8")
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
#---OS---#-------------------------
##---
filename.OSDFS.exp=paste("The data for Figure S1C and D.txt",sep = "")
rt.OSDFS=read.table(filename.OSDFS.exp,header=T, row.names=1,sep="\t")                   #读取输入文件
rt=rt.OSDFS[,-(3:4)]
res.cut=surv_cutpoint(rt, time = "Follow_days", event = "Follow_status", variables = gene)
cutoff=as.numeric(res.cut$cutpoint[1])
print(cutoff)
Median=median(rt[,gene])
print(Median)
Type1=ifelse(rt[,gene]<=Median, "Low", "High")
Type2=ifelse(rt[,gene]<=cutoff, "Low", "High")
rt$risk=Type1
rt$riskCutoff=Type2
outTaball=rbind(ID=c("futime","fustat","riskscore","riskMedian","riskCutoff"), rt)
colnames(outTaball)=outTaball[1,]
outTaball=outTaball[-1,]
outTaball$futime[outTaball$futime <0 ] = 0
outTaball$futime=as.numeric(outTaball$futime)
outTaball$futime=outTaball$futime/365
#Median
outTabMedian=outTaball[,c("futime","fustat","riskscore","riskMedian")]
outTabMedian=rbind(ID=c("futime","fustat","riskscore","risk"), outTabMedian)
filename.OS.riskMedian<- paste(gene,"-OS.firebrowse.down-Median.txt",sep = "")
write.table(outTabMedian, file=filename.OS.riskMedian, sep="\t", quote=F, col.names=F)
##--Median------
#
data=read.table(filename.OS.riskMedian,header=T, row.names=1,sep="\t")                   #读取输入文件
#10years
data$fustat[data$futime>10 ] <- 0
data$futime[data$futime>10 ] <- 10
cli="OS"
surDatarisk = data
#---
diffrisk=survdiff(Surv(futime, fustat) ~ risk, data = surDatarisk)
lengthrisk=length(levels(factor(surDatarisk[,"risk"])))
pValuerisk=1-pchisq(diffrisk$chisq, df=1)
#if(pValuerisk<0.01){
#  pValuerisk="p<0.01"
#}else{
#  pValuerisk=paste0("p=",sprintf("%.03f",pValuerisk))
#}
pValuerisk=signif(pValuerisk,4)
pValuerisk=format(pValuerisk, scientific = TRUE)
pValuerisk.OS.Median=pValuerisk
fit.risk <- survfit(Surv(futime, fustat) ~ risk, data = surDatarisk)
#print(surv_median(fit))
#????????????
bioCol=MyPal
bioColrisk=bioCol[1:length(levels(factor(surDatarisk[,"risk"])))]
surPlotrisk=ggsurvplot(fit.risk,
data=surDatarisk,
conf.int=F,
pval=pValuerisk,
pval.size=6,
legend.labs=levels(factor(surDatarisk[,"risk"])),
legend = c(0.8, 0.8),
legend.title=paste(gene),
font.legend=12,
xlab="Time(years)",
break.time.by = 1,
palette = bioCol,
surv.median.line = "hv",
risk.table=T,
cumevents=F,
risk.table.height=.25)
width=7
height=5.5
filename.tcga.OS.RiskMedian.pdf <- paste("OS.pdf",sep = "")
pdf(file=filename.tcga.OS.RiskMedian.pdf,onefile = FALSE, width=width, height=height)
print(surPlotrisk)
dev.off()
#---DFS---#-------------------------
filename.OSDFS.exp=paste("The data for Figure S1C and D.txt",sep = "")
rt.OSDFS=read.table(filename.OSDFS.exp,header=T, row.names=1,sep="\t")                   #读取输入文件
rt=rt.OSDFS[,-(1:2)]
res.cut=surv_cutpoint(rt, time = "DFS_days", event = "DFS_status", variables = gene)
cutoff=as.numeric(res.cut$cutpoint[1])
print(cutoff)
Median=median(rt[,gene])
print(Median)
Type1=ifelse(rt[,gene]<=Median, "Low", "High")
Type2=ifelse(rt[,gene]<=cutoff, "Low", "High")
rt$risk=Type1
rt$riskCutoff=Type2
outTaball=rbind(ID=c("futime","fustat","riskscore","riskMedian","riskCutoff"), rt)
colnames(outTaball)=outTaball[1,]
outTaball=outTaball[-1,]
outTaball$futime[outTaball$futime <0 ] = 0
outTaball$futime=as.numeric(outTaball$futime)
outTaball$futime=outTaball$futime/365
#--输出Median的分组
outTabMedian=outTaball[,c("futime","fustat","riskscore","riskMedian")]
outTabMedian=rbind(ID=c("futime","fustat","riskscore","risk"), outTabMedian)
filename.DFS.riskMedian<- paste(geneName,"-DFS.firebrowse.down-Median.txt",sep = "")
write.table(outTabMedian, file=filename.DFS.riskMedian, sep="\t", quote=F, col.names=F)
##--Median的分组
#比较高低风险组生存差异，得到显著性p值TCGA取10年，GEO取5年
data=read.table(filename.DFS.riskMedian,header=T, row.names=1,sep="\t")                   #读取输入文件
#对时间选择10年
data$fustat[data$futime>10 ] <- 0
data$futime[data$futime>10 ] <- 10
cli="DFS"
surDatarisk = data
#---生存分析
diffrisk=survdiff(Surv(futime, fustat) ~ risk, data = surDatarisk)
lengthrisk=length(levels(factor(surDatarisk[,"risk"])))
pValuerisk=1-pchisq(diffrisk$chisq, df=1)
#if(pValuerisk<0.01){
#  pValuerisk="p<0.01"
#}else{
#  pValuerisk=paste0("p=",sprintf("%.03f",pValuerisk))
#}
pValuerisk=signif(pValuerisk,4)
pValuerisk=format(pValuerisk, scientific = TRUE)
pValuerisk.DFS.Median=pValuerisk
fit.risk <- survfit(Surv(futime, fustat) ~ risk, data = surDatarisk)
#print(surv_median(fit))
#????????????
bioCol=MyPal
bioColrisk=bioCol[1:length(levels(factor(surDatarisk[,"risk"])))]
surPlotrisk=ggsurvplot(fit.risk,
data=surDatarisk,
conf.int=F,
pval=pValuerisk,
pval.size=6,
legend.labs=levels(factor(surDatarisk[,"risk"])),
legend = c(0.8, 0.8),
legend.title=paste(geneName),
font.legend=12,
xlab="Time(years)",
break.time.by = 1,
palette = bioCol,
surv.median.line = "hv",
risk.table=T,
cumevents=F,
risk.table.height=.25)
width=7
height=5.5
filename.tcga.DFS.RiskMedian.pdf <- paste("DFS.pdf",sep = "")
pdf(file=filename.tcga.DFS.RiskMedian.pdf,onefile = FALSE, width=width, height=height)
print(surPlotrisk)
dev.off()
#if (!requireNamespace("BiocManager", quietly = TRUE))
#    install.packages("BiocManager")
#BiocManager::install("limma", version = "3.8")
library("limma")
gene=geneName="FTL"                                                           #需要提取的基因
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
#---OS---#-------------------------
##---
filename.OSDFS.exp=paste("The data for Figure S1C and D.txt",sep = "")
rt.OSDFS=read.table(filename.OSDFS.exp,header=T, row.names=1,sep="\t")                   #读取输入文件
rt=rt.OSDFS[,-(3:4)]
res.cut=surv_cutpoint(rt, time = "Follow_days", event = "Follow_status", variables = gene)
cutoff=as.numeric(res.cut$cutpoint[1])
print(cutoff)
Median=median(rt[,gene])
print(Median)
Type1=ifelse(rt[,gene]<=Median, "Low", "High")
Type2=ifelse(rt[,gene]<=cutoff, "Low", "High")
rt$risk=Type1
rt$riskCutoff=Type2
outTaball=rbind(ID=c("futime","fustat","riskscore","riskMedian","riskCutoff"), rt)
colnames(outTaball)=outTaball[1,]
outTaball=outTaball[-1,]
outTaball$futime[outTaball$futime <0 ] = 0
outTaball$futime=as.numeric(outTaball$futime)
outTaball$futime=outTaball$futime/365
#Median
outTabMedian=outTaball[,c("futime","fustat","riskscore","riskMedian")]
outTabMedian=rbind(ID=c("futime","fustat","riskscore","risk"), outTabMedian)
filename.OS.riskMedian<- paste(gene,"-OS.firebrowse.down-Median.txt",sep = "")
write.table(outTabMedian, file=filename.OS.riskMedian, sep="\t", quote=F, col.names=F)
##--Median------
#
data=read.table(filename.OS.riskMedian,header=T, row.names=1,sep="\t")                   #读取输入文件
#10years
data$fustat[data$futime>10 ] <- 0
data$futime[data$futime>10 ] <- 10
cli="OS"
surDatarisk = data
#---
diffrisk=survdiff(Surv(futime, fustat) ~ risk, data = surDatarisk)
lengthrisk=length(levels(factor(surDatarisk[,"risk"])))
pValuerisk=1-pchisq(diffrisk$chisq, df=1)
#if(pValuerisk<0.01){
#  pValuerisk="p<0.01"
#}else{
#  pValuerisk=paste0("p=",sprintf("%.03f",pValuerisk))
#}
pValuerisk=signif(pValuerisk,4)
pValuerisk=format(pValuerisk, scientific = TRUE)
pValuerisk.OS.Median=pValuerisk
fit.risk <- survfit(Surv(futime, fustat) ~ risk, data = surDatarisk)
#print(surv_median(fit))
#????????????
bioCol=MyPal
bioColrisk=bioCol[1:length(levels(factor(surDatarisk[,"risk"])))]
surPlotrisk=ggsurvplot(fit.risk,
data=surDatarisk,
conf.int=F,
pval=pValuerisk,
pval.size=6,
legend.labs=levels(factor(surDatarisk[,"risk"])),
legend = c(0.8, 0.8),
legend.title=paste(gene),
font.legend=12,
xlab="Time(years)",
break.time.by = 1,
palette = bioCol,
surv.median.line = "hv",
risk.table=T,
cumevents=F,
risk.table.height=.25)
width=7
height=5.5
filename.tcga.OS.RiskMedian.pdf <- paste("OS.pdf",sep = "")
pdf(file=filename.tcga.OS.RiskMedian.pdf,onefile = FALSE, width=width, height=height)
print(surPlotrisk)
dev.off()
#---DFS---#-------------------------
filename.OSDFS.exp=paste("The data for Figure S1C and D.txt",sep = "")
rt.OSDFS=read.table(filename.OSDFS.exp,header=T, row.names=1,sep="\t")                   #读取输入文件
rt=rt.OSDFS[,-(1:2)]
res.cut=surv_cutpoint(rt, time = "DFS_days", event = "DFS_status", variables = gene)
cutoff=as.numeric(res.cut$cutpoint[1])
print(cutoff)
Median=median(rt[,gene])
print(Median)
Type1=ifelse(rt[,gene]<=Median, "Low", "High")
Type2=ifelse(rt[,gene]<=cutoff, "Low", "High")
rt$risk=Type1
rt$riskCutoff=Type2
outTaball=rbind(ID=c("futime","fustat","riskscore","riskMedian","riskCutoff"), rt)
colnames(outTaball)=outTaball[1,]
outTaball=outTaball[-1,]
outTaball$futime[outTaball$futime <0 ] = 0
outTaball$futime=as.numeric(outTaball$futime)
outTaball$futime=outTaball$futime/365
#--输出Median的分组
outTabMedian=outTaball[,c("futime","fustat","riskscore","riskMedian")]
outTabMedian=rbind(ID=c("futime","fustat","riskscore","risk"), outTabMedian)
filename.DFS.riskMedian<- paste(geneName,"-DFS.firebrowse.down-Median.txt",sep = "")
write.table(outTabMedian, file=filename.DFS.riskMedian, sep="\t", quote=F, col.names=F)
##--Median的分组
#比较高低风险组生存差异，得到显著性p值TCGA取10年，GEO取5年
data=read.table(filename.DFS.riskMedian,header=T, row.names=1,sep="\t")                   #读取输入文件
#对时间选择10年
data$fustat[data$futime>10 ] <- 0
data$futime[data$futime>10 ] <- 10
cli="DFS"
surDatarisk = data
#---生存分析
diffrisk=survdiff(Surv(futime, fustat) ~ risk, data = surDatarisk)
lengthrisk=length(levels(factor(surDatarisk[,"risk"])))
pValuerisk=1-pchisq(diffrisk$chisq, df=1)
#if(pValuerisk<0.01){
#  pValuerisk="p<0.01"
#}else{
#  pValuerisk=paste0("p=",sprintf("%.03f",pValuerisk))
#}
pValuerisk=signif(pValuerisk,4)
pValuerisk=format(pValuerisk, scientific = TRUE)
pValuerisk.DFS.Median=pValuerisk
fit.risk <- survfit(Surv(futime, fustat) ~ risk, data = surDatarisk)
#print(surv_median(fit))
#????????????
bioCol=MyPal
bioColrisk=bioCol[1:length(levels(factor(surDatarisk[,"risk"])))]
surPlotrisk=ggsurvplot(fit.risk,
data=surDatarisk,
conf.int=F,
pval=pValuerisk,
pval.size=6,
legend.labs=levels(factor(surDatarisk[,"risk"])),
legend = c(0.8, 0.8),
legend.title=paste(geneName),
font.legend=12,
xlab="Time(years)",
break.time.by = 1,
palette = bioCol,
surv.median.line = "hv",
risk.table=T,
cumevents=F,
risk.table.height=.25)
width=7
height=5.5
filename.tcga.DFS.RiskMedian.pdf <- paste("DFS.pdf",sep = "")
pdf(file=filename.tcga.DFS.RiskMedian.pdf,onefile = FALSE, width=width, height=height)
print(surPlotrisk)
dev.off()
