2022-02-20

請(qǐng)問(wèn)大家,在運(yùn)行R語(yǔ)言的時(shí)候,出現(xiàn)了以下問(wèn)題,不知道是源文件的問(wèn)題,還是程序的問(wèn)題,求解

源文件如下

There were 50 or more warnings (use warnings() to see the first 50)

#install.packages("survival")

>?

> library(survival)? ? ? ? ? ? ?#引用包

> inputFile="expTime.txt"? ? ? ?#輸入文件

> pFilter=0.001? ? ? ? ? ? ? ? #顯著性過(guò)濾條件

> setwd("C:\\Users\\ASUS\\Desktop\\tmeimmune\\22.cox")? ? ? ?#設(shè)置共工作目錄

>?

> #讀取輸入文件

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

> rt$futime=rt$futime/365? ? ? ? ? ?#以年為單位,除以365

>?

> outTab=data.frame()

> for(gene in colnames(rt[,3:ncol(rt)])){

+? ? ?#單因素cox分析

+? ? ?cox=coxph(Surv(futime, fustat) ~ rt[,gene], data = rt)

+? ? ?coxSummary = summary(cox)

+? ? ?coxP=coxSummary$coefficients[,"Pr(>|z|)"]

+? ? ?

+? ? ?#KM檢驗(yàn),后續(xù)需要可視乎,要和生存相關(guān)

+? ? ?group=ifelse(rt[,gene]<=median(rt[,gene]),"Low","High")

+? ? ?diff=survdiff(Surv(futime, fustat) ~ group,data = rt)

+? ? ?pValue=1-pchisq(diff$chisq,df=1)

+? ? ?

+? ? ?#保存滿足條件的基因

+? ? ?if((pValue<pFilter) & (coxP<pFilter)){

+? ? ? ? ?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"],

+? ? ? ? ? ? ? ? ? ? ? ? ? ? pvalue=coxP) )

+? ? ?}

+ }

There were 50 or more warnings (use warnings() to see the first 50)

> #輸出基因和P值表格文件

> write.table(outTab,file="cox.result.txt",sep="\t",row.names=F,quote=F)

>?

> #繪制森林圖

> #讀取輸入文件

> rt <- read.table("cox.result.txt",header=T,sep="\t",row.names=1,check.names=F)

Error in read.table("cox.result.txt", header = T, sep = "\t", row.names = 1,? :?

? no lines available in input

> gene <- rownames(rt)

> hr <- sprintf("%.3f",rt$"HR")

> hrLow? <- sprintf("%.3f",rt$"HR.95L")

> hrHigh <- sprintf("%.3f",rt$"HR.95H")

> Hazard.ratio <- paste0(hr,"(",hrLow,"-",hrHigh,")")

> pVal <- ifelse(rt$pvalue<0.001, "<0.001", sprintf("%.3f", rt$pvalue))

>?

> #輸出圖形

> pdf(file="forest.pdf", width = 7,height = 6)

> n <- nrow(rt)

> nRow <- n+1

> ylim <- c(1,nRow)

> layout(matrix(c(1,2),nc=2),width=c(3,2))

>?

> #繪制森林圖左邊的基因信息

> xlim = c(0,3)

> par(mar=c(4,2.5,2,1))

> plot(1,xlim=xlim,ylim=ylim,type="n",axes=F,xlab="",ylab="")

> text.cex=0.8

> text(0,n:1,gene,adj=0,cex=text.cex)

> text(1.5-0.5*0.2,n:1,pVal,adj=1,cex=text.cex);text(1.5-0.5*0.2,n+1,'pvalue',cex=text.cex,font=2,adj=1)

Error in text.default(1.5 - 0.5 * 0.2, n:1, pVal, adj = 1, cex = text.cex) :?

? 'labels'長(zhǎng)度不能設(shè)成零

> text(3,n:1,Hazard.ratio,adj=1,cex=text.cex);text(3,n+1,'Hazard ratio',cex=text.cex,font=2,adj=1,)

>?

> #繪制森林圖

> par(mar=c(4,1,2,1),mgp=c(2,0.5,0))

> xlim = c(0,max(as.numeric(hrLow),as.numeric(hrHigh)))

Warning message:

In max(as.numeric(hrLow), as.numeric(hrHigh)) :

? no non-missing arguments to max; returning -Inf

> plot(1,xlim=xlim,ylim=ylim,type="n",axes=F,ylab="",xaxs="i",xlab="Hazard ratio")

Error in plot.window(...) : 'xlim'值不能是無(wú)限的

> arrows(as.numeric(hrLow),n:1,as.numeric(hrHigh),n:1,angle=90,code=3,length=0.05,col="darkblue",lwd=2.5)

Error in arrows(as.numeric(hrLow), n:1, as.numeric(hrHigh), n:1, angle = 90,? :?

? 不能將零長(zhǎng)度的座標(biāo)同其它長(zhǎng)度的座標(biāo)混合在一起

> abline(v=1,col="black",lty=2,lwd=2)

> boxcolor = ifelse(as.numeric(hr) > 1, 'red', 'green')

> points(as.numeric(hr), n:1, pch = 15, col = boxcolor, cex=1.3)

Error in xy.coords(x, y) : 'x' and 'y' lengths differ

> axis(1)

> dev.off()

null device?

? ? ? ? ? 1?

最后編輯于
?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請(qǐng)聯(lián)系作者
【社區(qū)內(nèi)容提示】社區(qū)部分內(nèi)容疑似由AI輔助生成,瀏覽時(shí)請(qǐng)結(jié)合常識(shí)與多方信息審慎甄別。
平臺(tái)聲明:文章內(nèi)容(如有圖片或視頻亦包括在內(nèi))由作者上傳并發(fā)布,文章內(nèi)容僅代表作者本人觀點(diǎn),簡(jiǎn)書系信息發(fā)布平臺(tái),僅提供信息存儲(chǔ)服務(wù)。

相關(guān)閱讀更多精彩內(nèi)容

友情鏈接更多精彩內(nèi)容