有一组数据的表达量全是0,删掉就可以了
> library(survival)
> rt=read.table("mutSurvInput.txt",header=T,sep="\t",check.names=F)
> rt$futime=rt$futime/365
> outTab=data.frame()
>
> for(gene in colnames(rt[,4:ncol(rt)])){
+ a=rt[,gene]
+ diff=survdiff(Surv(futime, fustat) ~a,data = rt)
+ pValue=1-pchisq(diff$chisq,df=1)
+ outTab=rbind(outTab,cbind(gene=gene,pvalue=pValue))
+ if(pValue<0.001){
+ pValue=format(pValue, scientific = TRUE)
+ }else{
+ pValue=sprintf("%.3f",pValue)
+ }
+
+ fit <- survfit(Surv(futime, fustat) ~ a, data = rt)
+ summary(fit)
+
+ pdf(file=paste(gene,".survival.pdf",sep=""),
+ width = 5.5,
+ height =5,
+ )
+ plot(fit,
+ lwd=2,
+ col=c("red","blue"),
+ xlab="Time (year)",
+ mark.time=T,
+ ylim=c(0,1.05),
+ ylab="Survival rate",
+ main=paste(gene,"(p=", pValue ,")",sep="") )
+ legend("topright",
+ c("Mutation","Wild"),
+ lwd=2,
+ col=c("red","blue"))
+ dev.off()
+ }
Error in survdiff.fit(y, groups, strata.keep, rho) :
There is only 1 group