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)
报错图片如下
谢谢!
你好,我是有问必答小助手,非常抱歉,本次您提出的有问必答问题,技术专家团超时未为您做出解答
本次提问扣除的有问必答次数,将会以问答VIP体验卡(1次有问必答机会、商城购买实体图书享受95折优惠)的形式为您补发到账户。
因为有问必答VIP体验卡有效期仅有1天,您在需要使用的时候【私信】联系我,我会为您补发。