TCGA数据批量运行Coxph函数

news2025/2/25 21:35:19

df数据框形如:
在这里插入图片描述

djs.coxph <- function(df,genelist){
  library(survival)
  library(survminer)
  
  dir.create("./survival")
  setwd("./survival")
  
  # 准备好的生存分析数据框,变量中包括OS.time,OS以及values of gene expression 
  df <- as.data.frame(df)
  genelist <- genelist
  
  
  # 生成文件头,用于保存cox分析结果
  colname<-c("gene","beta", "HR (95% CI for HR)", "wald.test", "p.value")
  write.table(t(colname),file="./summary_HR.csv",
              sep=",",append=T,col.names=F,row.names=F)
  
  # 对每一个gene运行coxph
  
  lapply(genelist, function(x){
    univ_formulas <- paste('Surv(OS.time, OS)~', x)
    univ_models <- coxph(x, data = df)
    
    x <- summary(univ_models)
    p.value<-signif(x$wald["pvalue"], digits=2)
    wald.test<-signif(x$wald["test"], digits=2)
    beta<-signif(x$coef[1], digits=2);#coeficient beta
    HR <-signif(x$coef[2], digits=2);#exp(beta)
    HR.confint.lower <- signif(x$conf.int[,"lower .95"], 2)
    HR.confint.upper <- signif(x$conf.int[,"upper .95"],2)
    HR <- paste0(HR, " (", 
                 HR.confint.lower, "-", HR.confint.upper, ")")
    res<-c(i,beta, HR, wald.test, p.value)
    
    # 写入结果
    write.table(t(res),file="./summary_HR.csv",
                sep=",",append=T,col.names=F,row.names=F)
  })
}
djs.KMplot <- function(df,genelist,group){
  library(survival)
  library(survminer)
  
  dir.create("./survival")
  setwd("./survival")
  
  # 准备好的生存分析数据框,变量中包括OS.time,OS以及values of gene expression 
  df <- as.data.frame(df)
  genelist <- genelist
  group <- group
  
  # 判断使用那种分组方法
  if(group == "median"){
    lapply(genelist, function(x){
      df$group <- ifelse(df[,x] >= median(df[,x]),"high","low")
      # KMplot
      fit <- survfit(Surv(OS.time, OS) ~ group,data = df)
      a <- ggsurvplot(fit,
                      pval = TRUE,
                      conf.int=TRUE,
                      pval.size=5,
                      xlab=i,
                      palette=c("red", "blue"),
                      legend.labs=c("High", "Low"),
                      risk.table=T,
                      risk.table.height=.25)
      
      b <- surv_pvalue(fit)
      # 输出pvalue of logrank test
      write.table(t(c(x,b$pval)),file="./pvalue.of.survivaldata.csv",
                  sep=",",append=T,col.names=F,row.names=F)
      # 输出 风险事件表
      write.table(x,file="./risktable.of.survivaldata.csv",
                  sep=",",append=T,col.names=F,row.names=F)
      write.table(a$data.survtable,file="./risktable.of.survivaldata.csv",
                  sep=",",append=T,col.names=T,row.names=F)
      # 输出KMplot
      png(paste("./",x,"_survival.png",sep = ""))
      print(a)
      dev.off() 
    })
  }
  
  if(group == "mean"){
    lapply(genelist, function(x){
      df$group <- ifelse(df[,x] >= mean(df[,x]),"high","low")
      # KMplot
      fit <- survfit(Surv(OS.time, OS) ~ group,data = df)
      a <- ggsurvplot(fit,
                      pval = TRUE,
                      conf.int=TRUE,
                      pval.size=5,
                      xlab=i,
                      palette=c("red", "blue"),
                      legend.labs=c("High", "Low"),
                      risk.table=T,
                      risk.table.height=.25)
      
      b <- surv_pvalue(fit)
      # 输出pvalue of logrank test
      write.table(t(c(x,b$pval)),file="./pvalue.of.survivaldata.csv",
                  sep=",",append=T,col.names=F,row.names=F)
      # 输出 风险事件表
      write.table(x,file="./risktable.of.survivaldata.csv",
                  sep=",",append=T,col.names=F,row.names=F)
      write.table(a$data.survtable,file="./risktable.of.survivaldata.csv",
                  sep=",",append=T,col.names=T,row.names=F)
      # 输出KMplot
      png(paste("./",x,"_survival.png",sep = ""))
      print(a)
      dev.off() 
    })
  }
  
  if(group == "quantile"){
      lapply(genelist, function(x){
        df$group <- ifelse(df[,x] >= quantile(df[,x])[[4]],"high",
                           ifelse(df[,x] <= quantile(df[,x])[[2]],"low","undetermine"))
        # KMplot
        fit <- survfit(Surv(OS.time, OS) ~ group,data = df[df$group != "undetermine",])
        a <- ggsurvplot(fit,
                        pval = TRUE,
                        conf.int=TRUE,
                        pval.size=5,
                        xlab=i,
                        palette=c("red", "blue"),
                        legend.labs=c("High", "Low"),
                        risk.table=T,
                        risk.table.height=.25)
        
        b <- surv_pvalue(fit)
        # 输出pvalue of logrank test
        write.table(t(c(x,b$pval)),file="./pvalue.of.survivaldata.csv",
                    sep=",",append=T,col.names=F,row.names=F)
        # 输出 风险事件表
        write.table(x,file="./risktable.of.survivaldata.csv",
                    sep=",",append=T,col.names=F,row.names=F)
        write.table(a$data.survtable,file="./risktable.of.survivaldata.csv",
                    sep=",",append=T,col.names=T,row.names=F)
        # 输出KMplot
        png(paste("./",x,"_survival.png",sep = ""))
        print(a)
        dev.off() 
      })
  }
}

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处:http://www.coloradmin.cn/o/820580.html

如若内容造成侵权/违法违规/事实不符,请联系多彩编程网进行投诉反馈,一经查实,立即删除!

相关文章

论文阅读 - Few-shot Network Anomaly Detection via Cross-network Meta-learning

论文链接&#xff1a;https://arxiv.org/pdf/2102.11165.pdf 目录 摘要&#xff1a; 引言 问题定义 方法 Graph Deviation Networks Cross-network Meta-learning 摘要&#xff1a; 网络异常检测旨在找到与绝大多数行为显着不同的网络元素&#xff08;例如节点、边、子图…

河北沃克仓储解决方案最新布局|HEGERLS四向穿梭车在高标仓和楼层仓中的应用

河北沃克金属制品有限公司是业内十分稀缺可提供整体物流仓储解决方案并落地的企业&#xff0c;既拥有自主研发和生产等一整套核心软硬件的能力&#xff0c;又具备丰富的整体方案规划与实施经验。经过多年积累和开拓&#xff0c;河北沃克金属制品有限公司业务已拓展覆盖近众多行…

如何列出phpMyAdmin左侧菜单中的所有表格 - 不分页 - 显示数据库所有的表

效果图 原来 优化后 步骤 点击logo&#xff0c;回到首页 点击设置 》 导航面板 》 导航树 》 节点中最大项数 》 应用 ok

Windows下RocketMQ的启动

下载地址&#xff1a;下载 | RocketMQ 解压后 一、修改runbroker.cmd 修改 bin目录下的runbroker.cmd set "JAVA_OPT%JAVA_OPT% -server -Xms2g -Xmx2g" set "JAVA_OPT%JAVA_OPT% -XX:MaxDirectMemorySize15g" set "JAVA_OPT%JAVA_OPT% -cp %CLASSP…

jmeter之接口测试(http接口测试)

基础知识储备 一、了解jmeter接口测试请求接口的原理 客户端--发送一个请求动作--服务器响应--返回客户端 客户端--发送一个请求动作--jmeter代理服务器---服务器--jmeter代理服务器--服务器 二、了解基础接口知识&#xff1a; 1、什么是接口&#xff1a;前端与后台之间的…

MySQL 的 Join 查询及 Hash Join 优化 | StoneDB 技术分享会 #3

StoneDB开源地址 https://github.com/stoneatom/stonedb 设计&#xff1a;小艾 审核&#xff1a;丁奇、宇亭 编辑&#xff1a;宇亭 作者一&#xff1a;徐鑫强&#xff08;花名&#xff1a;无花果&#xff09; 电子科技大学-计算机技术-在读硕士、StoneDB 内核研发实习生 作…

BES 平台 SDK之代码架构讲解二

本文章是基于BES2700 芯片&#xff0c;其他BESxxx 芯片可做参考&#xff0c;如有不当之处&#xff0c;欢迎评论区留言指出。 BES 平台 SDK之代码架构讲解一_谢文浩的博客-CSDN博客 上篇文章粗略的对整个SDK 目录下的文件进行了说明&#xff0c;接下来会对SDK 比较详细的介绍。…

C语言实用调试详解

目录 什么是bug? 调试是什么?有多重要? 调试是什么? 调试的基本步骤 Debug和Release的介绍 Windows环境调试介绍 调试环境的准备 学会快捷键 调试的时候查看程序当前信息 查看临时变量的值 查看内存信息 查看调用堆栈 查看汇编信息 查看寄存器信息 一些调试…

HCIP 重发布+路由策略总结

重发布 在同一个网络拓结构中&#xff0c;如果存在多种不同的路由协议&#xff0c;由于不同路由协议的机制各有不同&#xff0c;对路由的处理也不相同&#xff0c;这就在网络中造成了路由信息的隔离&#xff0c;在路由器的边界路由器上&#xff0c;将某种路由协议的路由信息引…

[网络工程师]-网络规划与设计-网络故障分析与处理

网络环境越复杂,发生故障的可能性越大,引发故障的原因也就越难确定。网络故障往往具有特定的故障现象。这些现象可能比较笼统,也可能比较特殊。利用特定的故障排查工具及技巧,在具体的网络环境下观察故障现象,细致分析,最终必然可以查出一个或多个引发故障的原因。一旦能…

gitlab上传代码

输入 git clone https地址&#xff0c;此地址可以在&#xff0c;gitlab项目上拷贝到本地&#xff0c;看本地电脑会出现在gitlab上新建的项目&#xff0c;并进入该目录下 将要上传的代码拷贝到该目录 依次输入一下代码 git init &#xff08;用于在目录中创建新的 Git 仓库。…

打造独一无二的花店小程序,轻松搭建步骤详解

随着移动互联网的快速发展&#xff0c;花店也开始意识到拥有一个专属的小程序能够提升用户体验、增加销售额。那么&#xff0c;如何快速搭建一个漂亮、实用的花店小程序呢&#xff1f;下面就为大家介绍一下具体的步骤。 第一步&#xff0c;使用第三方制作平台。如乔拓云网是一个…

setEagerlyType字段理解

官方文档介绍&#xff1a;V5.0.4版本开始一对一关联预载入支持两种方式&#xff1a;JOIN方式&#xff08;一次查询&#xff09;和IN方式&#xff08;两次查询&#xff09;&#xff0c;如果要使用IN方式关联预载入&#xff0c;在关联定义方法中添加。 这句话的意思是jion方式关联…

阿里云服务器免费试用及搭建WordPress网站

文章目录 前言一、免费试用1、选择使用产品2、进行产品配置3、远程连接阿里云服务器①、重置实例密码②、SecureCRT 远程链接③、Workbench 远程链接二、搭建 WordPress 网站1、开放搭建 WordPress 需要的端口2、搭建 LAMP 环境①、Linux 系统升级和更新源②、安装 Apache2③、…

【Excel】记录Match和Index函数的用法

最近一直用到的两个处理EXCEL表格数据的函数向大家介绍一下&#xff0c;写这篇博文的目的也是为了记录免得自己忘记了&#xff0c;嘻嘻。 先上百度的链接 Match函数的用法介绍&#xff1a;https://jingyan.baidu.com/article/2fb0ba40b4933941f3ec5f71.html 小结&#xff1a;…

Java从入门到精通(二)· 基本语法

Java从入门到精通&#xff08;二&#xff09; 基本语法 一 变量 1.字面量 计算机是用来处理数据的&#xff0c;字面量就是告诉程序员&#xff1a;数据在程序中的书写格式。 特殊的字符&#xff1a; \n 表示换行&#xff0c; \t 表示一个制表符&#xff0c;即一个tab 2.变量…

【JVM】什么是双亲委派机制

文章目录 1、类加载机制2、双亲委派模型2.1、介绍2.2、为什么需要双亲委派2.3、源码解析 3、破坏双亲委派3.1、介绍3.2、破坏实现3.3、破坏双亲委派的例子 4、线程上下文类加载器 1、类加载机制 类加载阶段分为加载、连接、初始化三个阶段&#xff0c;而加载阶段需要通过类的全…

【Ansible】Ansible自动化运维工具之playbook剧本搭建LNMP架构

LNMP 一、playbooks 分布式部署 LNMP1. 环境配置2. 安装 ansble3. 安装 nginx3.1 准备 nginx 相关文件3.2 编写 lnmp.yaml 的 nginx 部分3.3 测试 nginx4. 安装 mysql4.1 准备 mysql 相关文件4.2 编写 lnmp.yaml 的 mysql 部分4.3 测试 mysql5. 安装 php5.1 编写 lnmp.yaml 的 …

健身时戴什么耳机比较好、健身用的运动耳机推荐

运动健身已经成为一种潮流&#xff0c;有的人为了追求马甲线和八大块腹肌&#xff0c;还有的人为了缓解学习和工作的压力。但你在运动健身的时候难免会有烦躁、疲惫的时候&#xff0c;如果这时有音乐的加入那就完美了&#xff0c;因为美妙的歌声能冲淡这种疲惫感&#xff0c;让…

面向金融科技方向选手!一级学会背书,AI选股与可视分析大赛来啦

金融量化领域邂逅人工智能&#xff0c;将会迸发出怎样的火花&#xff1f; 在深度学习、强化学习和自然语言处理等技术取得不断突破和创新的今天&#xff0c;AI如何赋能量化投资领域&#xff0c;助力开发者打造表现优异&#xff0c;更加安全可靠的量化模型&#xff1f; 第四届CS…