boxplot(expression ~ clinical, data = rt1,names=xlabel,
ylab = paste0("Risk score"),xlab = paste0(clinical),col=dotCol,
cex.main=1.6, cex.lab=1.4, cex.axis=1.3,ylim=c(yMin,yMax),outline = FALSE)
segments(1,ySeg, n,ySeg);
segments(1,ySeg, 1,ySeg2)
segments(n,ySeg, n,ySeg2)
text((1+n)/2,ySeg,labels=paste0("p=",pval),cex=1.5,pos=3)
dev.off()
}
###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="riskscore"
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,".",gene,".noX-ggboxplot.w4h6.pdf"),width=4,height=6)
print(boxplot)
dev.off()
}
for(cli in colnames(rtmerge.order[,2:(ncol(rtmerge.order)-1)])){
#for(cli in colnames(rtmerge.order[2])){
#cli = colnames(rtmerge.order[3])
clinical=cli                                                       #定义临床类型
gene="riskscore"                                                            #输入基因名称
xlabel=vector()
tab1=table(rt[,clinical])
labelNum=length(tab1)
dotCol=c("#377EB8", "#E41A1C")
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.cliName=paste0("pValue.",cli,sep = "")
pValue=pValue.cliSub=wilcoxTest$p.value
outTab=rbind(outTab,cbind(gene=gene,pVal=pValue))
pval=0
if(pValue.cliSub<0.001){
pval=signif(pValue.cliSub,4)
pval=format(pval, scientific = TRUE)
}else{
pval=round(pValue.cliSub,3)
}
b = boxplot(expression ~ clinical, data = rt1,outline = FALSE, plot=F)
yMin=min(b$stats)
yMax = max(b$stats/5+b$stats)
ySeg = max(b$stats/10+b$stats)
ySeg2 = max(b$stats/12+b$stats)
n = ncol(b$stats)
pValue.cliSub=round(pValue.cliSub,3)
width=5
height=6
outPdf=paste0(geneName,"-",clinical,".","firebrowse.down-Median-Clin.AJCC-boxplot.h",height,"w",width,"-",pValue.cliSub,".pdf")
pdf(file=outPdf,
width=width,height=height)
par(mar = c(4,7,3,3))
boxplot(expression ~ clinical, data = rt1, names=xlabel,
ylab = ylab,xlab = paste0(clinical),col=dotCol,
cex.main=1.6, cex.lab=1.4, cex.axis=1.3,ylim=c(yMin,yMax),outline = FALSE)
segments(1,ySeg, n,ySeg);
segments(1,ySeg, 1,ySeg2)
segments(n,ySeg, n,ySeg2)
text((1+n)/2,ySeg,labels=paste0("p=",pval),cex=1.5,pos=3)
dev.off()
#}
###用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="riskscore"
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=ylab,
legend.title=clinical,
palette=bioCol
)+
stat_compare_means(comparisons=my_comparisons)
width=5
height=6
outPdf=paste0(geneName,"-",clinical,".","firebrowse.down-Median-Clin.AJCC-ggboxplot.h",height,"w",width,"-",pValue.cliSub,".pdf")
pdf(file=outPdf,
width=width,height=height)
print(boxplot)
dev.off()
#}
##用ggplot输出
#for(cli in colnames(rtmerge.order[,2:(ncol(rtmerge.order)-1)])){
#cli = colnames(rt)[11]
clinical=cli
gene="riskscore"
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=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
outPdf=paste0(geneName,"-",clinical,".","firebrowse.down-Median-Clin.AJCC-gboxp.h",height,"w",width,"-",pValue.cliSub,".pdf")
pdf(file=outPdf,
width=width,height=height)
print(plotcli)
dev.off()
}
#cli = colnames(rt)[11]
clinical=cli
gene="riskscore"
data=rt[c(gene,clinical)]
bioCol=c("#377EB8","#E41A1C","#FF9900","#6E568C","#7CC767","#223D6C","#D20A13","#FFD121","#088247","#11AA4D")
bioCol=bioCol[1:length(unique(data[,clinical]))]
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,".",gene,".noX-ggboxplot.w4h6.pdf"),width=4,height=6)
print(boxplot)
dev.off()
filename.Tcga.Risk.Acli.merge <- paste("The data for Figure 4-clinical information-riskscoreorder.txt",sep = "")
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="riskscore"
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()
}
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$riskscore=log2(rtrisk$riskscore+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
riskScoreorderltoh = order(rtmerge[,"riskscore"],decreasing=F)
rtmerge.order=rtmerge[riskScoreorderltoh,]
show(rtmerge.order$riskScore)
filename.Tcga.Risk.Acli.merge <- paste("The data for Figure 4-clinical information-riskscoreorder.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="riskscore"
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()
}
for(cli in colnames(rtmerge.order[,2:(ncol(rtmerge.order)-1)])){
#for(cli in colnames(rtmerge.order[2])){
#cli = colnames(rtmerge.order[3])
clinical=cli                                                       #定义临床类型
gene="riskscore"                                                            #输入基因名称
xlabel=vector()
tab1=table(rt[,clinical])
labelNum=length(tab1)
dotCol=c("#377EB8", "#E41A1C")
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.cliName=paste0("pValue.",cli,sep = "")
pValue=pValue.cliSub=wilcoxTest$p.value
outTab=rbind(outTab,cbind(gene=gene,pVal=pValue))
pval=0
if(pValue.cliSub<0.001){
pval=signif(pValue.cliSub,4)
pval=format(pval, scientific = TRUE)
}else{
pval=round(pValue.cliSub,3)
}
b = boxplot(expression ~ clinical, data = rt1,outline = FALSE, plot=F)
yMin=min(b$stats)
yMax = max(b$stats/5+b$stats)
ySeg = max(b$stats/10+b$stats)
ySeg2 = max(b$stats/12+b$stats)
n = ncol(b$stats)
pValue.cliSub=round(pValue.cliSub,3)
width=5
height=6
outPdf=paste0(geneName,"-",clinical,".","firebrowse.down-Median-Clin.AJCC-boxplot.h",height,"w",width,"-",pValue.cliSub,".pdf")
pdf(file=outPdf,
width=width,height=height)
par(mar = c(4,7,3,3))
boxplot(expression ~ clinical, data = rt1, names=xlabel,
ylab = ylab,xlab = paste0(clinical),col=dotCol,
cex.main=1.6, cex.lab=1.4, cex.axis=1.3,ylim=c(yMin,yMax),outline = FALSE)
segments(1,ySeg, n,ySeg);
segments(1,ySeg, 1,ySeg2)
segments(n,ySeg, n,ySeg2)
text((1+n)/2,ySeg,labels=paste0("p=",pval),cex=1.5,pos=3)
dev.off()
#}
###用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="riskscore"
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=ylab,
legend.title=clinical,
palette=bioCol
)+
stat_compare_means(comparisons=my_comparisons)
width=5
height=6
outPdf=paste0(geneName,"-",clinical,".","firebrowse.down-Median-Clin.AJCC-ggboxplot.h",height,"w",width,"-",pValue.cliSub,".pdf")
pdf(file=outPdf,
width=width,height=height)
print(boxplot)
dev.off()
#}
##用ggplot输出
#for(cli in colnames(rtmerge.order[,2:(ncol(rtmerge.order)-1)])){
#cli = colnames(rt)[11]
clinical=cli
gene="riskscore"
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=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
outPdf=paste0(geneName,"-",clinical,".","firebrowse.down-Median-Clin.AJCC-gboxp.h",height,"w",width,"-",pValue.cliSub,".pdf")
pdf(file=outPdf,
width=width,height=height)
print(plotcli)
dev.off()
}
#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$riskscore=log2(rtrisk$riskscore+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
riskScoreorderltoh = order(rtmerge[,"riskscore"],decreasing=F)
rtmerge.order=rtmerge[riskScoreorderltoh,]
show(rtmerge.order$riskScore)
filename.Tcga.Risk.Acli.merge <- paste("The data for Figure 4-clinical information-riskscoreorder.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="riskscore"
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()
}
