######Video source: https://shop119322454.taobao.com #install.packages("survival") setwd("D:\\tcgagene\\survival") #工作目录(需修改) library(survival) rt=read.table("clinicalExp.txt",header=T,sep="\t",check.names=F) rt$futime=rt$futime/365 #如果以月为单位,除以30;以年为单位,除以365 outTab=data.frame() for(gene in colnames(rt[,4:ncol(rt)])){ a=rt[,gene]<median(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)) #pValue=round(pValue,3) pValue=signif(pValue,4) pValue=format(pValue, scientific = TRUE) fit <- survfit(Surv(futime, fustat) ~ a, data = rt) summary(fit) tiff(file=paste(gene,".survival.tiff",sep=""), width = 14, #图片的宽度 height =14, #图片的高度 units ="cm", compression="lzw", bg="white", res=600) plot(fit, lwd=2, col=c("red","blue"), xlab="Time (month)", mark.time=T, ylab="Survival rate", main=paste("Survival curve (p=", pValue ,")",sep="")) legend("topright", c(paste(gene," high expression",sep=""), paste(gene," low expression",sep="") ), lwd=2, col=c("red","blue")) dev.off() } write.table(outTab,file="survival.xls",sep="\t",row.names=F,quote=F) ###Video source: http://study.163.com/u/biowolf ######Video source: https://shop119322454.taobao.com ######速科生物: http://www.biowolf.cn/ ######作者邮箱:2740881706@qq.com ######作者微信: seqBio ######QQ群: 259208034
时间: 2024-03-07 18:53:54 浏览: 110
这段代码是用R语言进行生存分析的,主要是对TCGA(The Cancer Genome Atlas)中的基因表达数据进行分析。具体分析过程如下:
1. 读取生存数据文件,将生存时间(futime)的单位转换为年;
2. 对每个基因进行分析,将样本按基因表达水平的中位数分为高表达组和低表达组;
3. 计算高表达组和低表达组之间的生存曲线差异,得到p值;
4. 绘制生存曲线图,并在图中标注p值;
5. 将每个基因的p值和生存曲线图输出到文件;
6. 最后将所有基因的p值输出到Excel文件中。
其中,生存分析是一种用于研究事件发生时间和事件相关因素之间关系的统计方法,常用于癌症等疾病的临床研究中。这段代码中使用了R语言中的survival包来实现生存分析。
阅读全文