双条件生存曲线

rt.png

library(limma)
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)
colnames(data)=gsub("(.*?)\\.(.*?)\\.(.*?)\\.(.*?)\\..*", "\\1\\-\\2\\-\\3", colnames(data))
data=t(data)

data.png

risk=read.table("risk.txt", header=T, sep="\t", check.names=F, row.names=1)

risk.txt.png

sameSample=intersect(row.names(risk), row.names(data))
risk=risk[sameSample, "risk",drop=F]
data=data[sameSample, , drop=F]
rt=cbind(risk, data)

rt.png

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

time.txt.png

#数据合并并输出结果
sameSample=intersect(row.names(rt),row.names(cli))
data=rt[sameSample,]
cli=cli[sameSample,]
out=cbind(cli,data)
out=cbind(id=row.names(out),out)
write.table(out,file="input.txt",sep="\t",row.names=F,quote=F)

input.txt.png

library(survival)
library(survminer)
inputFile="input.txt"         
outFile="survival.pdf"        
var1="AXCC"                   #用于生存分析的变量

rt=read.table(inputFile, header=T, sep="\t", check.names=F)

rt.png

#根据基因表达，对数据分组
rt$futime<-rt$futime/365
rt$risk<-paste("模型",rt$risk)
a=rt$risk
b=ifelse(rt[,var1]<median(rt[,var1]),paste0(var1," low"),paste0(var1," high"))
Type=paste(a,"+",b)
rt=cbind(rt,Type)

rt.png

#生存差异统计
length=length(levels(factor(Type)))
diff=survdiff(Surv(futime, fustat) ~Type,data = rt)
pValue=1-pchisq(diff$chisq,df=length-1)
if(pValue<0.001){
    pValue="p<0.001"
}else{
    pValue=paste0("p=",sprintf("%.03f",pValue))
}
fit <- survfit(Surv(futime, fustat) ~ Type, data = rt)
        
#绘制生存曲线
surPlot=ggsurvplot(fit, 
                   data=rt,
                   conf.int=F,  #不含置信区间
                   pval=pValue,
                   pval.size=6,
                   legend.labs=levels(factor(rt[,"Type"])),
                   legend.title="Type",
                   #font.legend=6,
                   legend = "right",
                   xlab="Time(years)",
                   break.time.by = 5, #x坐标轴间隔多少年
                   risk.table.title="",
                   risk.table=F,     #风险表格
                   risk.table.height=.25)
pdf(file=outFile,onefile = FALSE,width = 7.5,height =5)
print(surPlot)
dev.off()