不看KM

不看KM

背景知识

讲课讲到批量的logrank test可以告诉我们两个组之间的生存率差别是否显著，收到提问说：

老师，刚才的 生存分析中，展示了所有基因的p值，有没有其他列可以判断预后是好还是差？

这是一个好问题呀，如果只看p值，只是知道差别是否显著，那么到底是基因表达量高的组预后好，还是低的组预后好？

搜索发现我们可以通过比较两组的中位生存时间来量化。但是这个不适用于曲线有交叉的情况。

没关系，我们就只看p＜0.05的基因就可以啦。

怎么得到两组的中位生存时间呢

一段最少的代码，用内置数据lung来做

代码语言：javascript代码运行次数：0运行复制

library(survminer)
library(survival)
fit <- survfit(Surv(time, status) ~ , data = lung)
summary(fit)$table[,'median']

代码语言：javascript代码运行次数：0运行复制

## =1 =2 
##   270   426

代码语言：javascript代码运行次数：0运行复制

plot(fit)

用表达矩阵里的基因来做

使用我的tinyarray包里的示例数据

代码语言：javascript代码运行次数：0运行复制

library(tinyarray)
str(exprSet_hub1)

代码语言：javascript代码运行次数：0运行复制

##  num [1:8, 1:177] 8.0 19.1 17.6 12.06 17.64 ...
##  - attr(*, "dimnames")=List of 2
##   ..$ : chr [1:8] "CXCL8" "F1" "COLA1" "ISG15" ...
##   ..$ : chr [1:177] "TCGA-A-A9IO-01A" "TCGA-US-A774-01A" "TCGA-HZ-A49H-01A" "TCGA-FB-A4P5-01A" ...

代码语言：javascript代码运行次数：0运行复制

str(meta1)

代码语言：javascript代码运行次数：0运行复制

## 'data.frame':    177 obs. of  4 variables:
##  $ sample   : chr  "TCGA-A-A9IO-01A" "TCGA-US-A774-01A" "TCGA-HZ-A49H-01A" "TCGA-FB-A4P5-01A" ...
##  $ event    : int  0 1 0 1 0 1 0 0 1 0 ...
##  $ X_PATIET: chr  "TCGA-A-A9IO" "TCGA-US-A774" "TCGA-HZ-A49H" "TCGA-FB-A4P5" ...
##  $ time     : int  1942 695 491 179 228 18 289 92 476 1854 ...

写了个函数，输入数据是整理好的表达矩阵和临床信息，表达矩阵的列名和临床信息表格的行名是一一对应的。

代码语言：javascript代码运行次数：0运行复制

which_better = function(exprSet_hub,meta,pvalue_cutoff = 1){
  x = surv_KM(exprSet_hub,meta,pvalue_cutoff = pvalue_cutoff)
  diffs = apply(exprSet_hub,1,function(g){
    #g = 'ICAM1'
    gene = ifelse(g>median(g),"high","low")
    fit <- survfit(Surv(time, event) ~ gene, data = meta)
    me = summary(fit)$table[,'median']
    diff =ifelse(( me[1]-me[2])>0,"high","low")
    return(diff)
  })
  re = data.frame(p = x,
             better_group = diffs[names(x)])
  return(re)
  }
which_better(exprSet_hub1,meta1)

代码语言：javascript代码运行次数：0运行复制

##                p better_group
## ICAM1  0.0218004         high
## F1    0.498822          low
## COLA1 0.54286         high
## MMP9   0.6977895          low
## CXCL10 0.71901          low
## COL1A2 0.741749         high
## ISG15  0.959846          low
## CXCL8  0.9599565         high

如果只想看p<0.05的基因那么就加上参数：

代码语言：javascript代码运行次数：0运行复制

which_better(exprSet_hub1,meta1,pvalue_cutoff = 0.05)

代码语言：javascript代码运行次数：0运行复制

##               p better_group
## ICAM1 0.0218004         high

等我有空，把这个函数再精装修一下放进tinyarray里面。

画图检验成果

代码语言：javascript代码运行次数：0运行复制

p = exp_surv(exprSet_hub1,meta1)
library(patchwork)
wrap_plots(p,nrow = 2)

感谢提问的小朋友，祝大家科研顺利O(∩_∩)O。

本文参与 腾讯云自媒体同步曝光计划，分享自。原始发表：2024-10-22，如有侵权请联系 cloudcommunity@tencent 删除plot函数量化数据搜索