R生存分析代码,不知道是代码出现问题还是输入文件出现问题,报错!

R生存分析代码,不知道是代码出现问题还是输入文件出现问题,一直报错!自己检测感觉代码和输入文件都没有问题,麻烦各位伸出援手!

代码如下:

library(survival)
setwd("C:\\Users\\Administrator\\Desktop\\94geo\\07.survivalFilter")                 #工作目录(需修改)
pFilter=0.001                                                                      #显著性过滤条件
rt=read.table("expTime.txt",header=T,sep="\t",check.names=F,row.names=1)           #读取输入文件

outTab=data.frame()
sigGenes=c("futime","fustat")
for(gene in colnames(rt[,3:ncol(rt)])){
       if(sd(rt[,gene])<0.1){
         next}
       a=rt[,gene]<=median(rt[,gene]) 
       rt1=rt[a,]
       b=setdiff(rownames(rt),rownames(rt1))
       rt2=rt[b,]
       surTab1=summary(survfit(Surv(futime, fustat) ~ 1, data = rt1))
       surTab2=summary(survfit(Surv(futime, fustat) ~ 1, data = rt2))
       survivalTab1=cbind(time=surTab1$time, surv=surTab1$surv,lower=surTab1$lower,upper=surTab1$upper)
       survivalTab1=survivalTab1[survivalTab1[,"time"]<5,]
       if(class(survivalTab1)=="matrix"){
         survivalTab1=survivalTab1[nrow(survivalTab1),]
       }
       survivalTab2=cbind(time=surTab2$time, surv=surTab2$surv,lower=surTab2$lower,upper=surTab2$upper)
       survivalTab2=survivalTab2[survivalTab2[,"time"]<5,]
       if(class(survivalTab2)=="matrix"){
         survivalTab2=survivalTab2[nrow(survivalTab2),]
       }
       fiveYearsDiff=abs(survivalTab1["surv"]-survivalTab2["surv"])
       if(is.na(fiveYearsDiff)==TRUE){next}

       diff=survdiff(Surv(futime, fustat) ~a,data = rt)
       pValue=1-pchisq(diff$chisq,df=1)
       fit=survfit(Surv(futime, fustat) ~ a, data = rt)
       cox=coxph(Surv(futime, fustat) ~ rt[,gene], data = rt)
       coxSummary = summary(cox)
       coxP=coxSummary$coefficients[,"Pr(>|z|)"]
    
       if((pValue<pFilter) & (coxP<pFilter) & (fiveYearsDiff>0.1)){
           sigGenes=c(sigGenes,gene)
           outTab=rbind(outTab,
                       cbind(gene=gene,
                             KM=pValue,
                             HR=coxSummary$conf.int[,"exp(coef)"],
                             HR.95L=coxSummary$conf.int[,"lower .95"],
                             HR.95H=coxSummary$conf.int[,"upper .95"],
                                 coxPvalue=coxP) )
         }
}
write.table(outTab,file="survival.result.xls",sep="\t",row.names=F,quote=F)    #输出基因和p值表格文件
surSigExp=rt[,sigGenes]
surSigExp=cbind(id=row.names(surSigExp),surSigExp)
write.table(surSigExp,file="surSigExp.txt",sep="\t",row.names=F,quote=F)

报错图片如下

img

谢谢!

你好,我是有问必答小助手,非常抱歉,本次您提出的有问必答问题,技术专家团超时未为您做出解答


本次提问扣除的有问必答次数,将会以问答VIP体验卡(1次有问必答机会、商城购买实体图书享受95折优惠)的形式为您补发到账户。


因为有问必答VIP体验卡有效期仅有1天,您在需要使用的时候【私信】联系我,我会为您补发。