背景知识
讲课讲到批量的logrank test可以告诉我们两个组之间的生存率差别是否显著，收到提问说：
老师，刚才的 生存分析中，展示了所有基因的p值，有没有其他列可以判断预后是好还是差？
这是一个好问题呀，如果只看p值，只是知道差别是否显著，那么到底是基因表达量高的组预后好，还是低的组预后好？
搜索发现我们可以通过比较两组的中位生存时间来量化。但是这个不适用于曲线有交叉的情况。
没关系，我们就只看p＜0.05的基因就可以啦。
怎么得到两组的中位生存时间呢
一段最少的代码，用内置数据lung来做
library(survminer)
library(survival)
fit <- survfit(Surv(time, status) ~ sex, data = lung)
summary(fit)$table[,'median']
## sex=1 sex=2
## 270 426
plot(fit)
用表达矩阵里的基因来做
使用我的tinyarray包里的示例数据
library(tinyarray)
str(exprSet_hub1)
## num [1:8, 1:177] 8.03 19.13 17.63 12.06 17.64 ...
## - attr(*, "dimnames")=List of 2
## ..$ : chr [1:8] "CXCL8" "FN1" "COL3A1" "ISG15" ...
## ..$ : chr [1:177] "TCGA-3A-A9IO-01A" "TCGA-US-A774-01A" "TCGA-HZ-A49H-01A" "TCGA-FB-A4P5-01A" ...
str(meta1)
## 'data.frame': 177 obs. of 4 variables:
## $ sample : chr "TCGA-3A-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_PATIENT: chr "TCGA-3A-A9IO" "TCGA-US-A774" "TCGA-HZ-A49H" "TCGA-FB-A4P5" ...
## $ time : int 1942 695 491 179 228 183 289 392 476 1854 ...
写了个函数，输入数据是整理好的表达矩阵和临床信息，表达矩阵的列名和临床信息表格的行名是一一对应的。
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(as.numeric( 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)
## p better_group
## ICAM1 0.0218004 high
## FN1 0.4983822 low
## COL3A1 0.5428363 high
## MMP9 0.6977895 low
## CXCL10 0.7190313 low
## COL1A2 0.7417349 high
## ISG15 0.9598463 low
## CXCL8 0.9599565 high
如果只想看p<0.05的基因那么就加上参数：
which_better(exprSet_hub1,meta1,pvalue_cutoff = 0.05)
## p better_group
## ICAM1 0.0218004 high
等我有空，把这个函数再精装修一下放进tinyarray里面。
画图检验成果
p = exp_surv(exprSet_hub1,meta1)
library(patchwork)
wrap_plots(p,nrow = 2)
感谢提问的小朋友，祝大家科研顺利O(∩_∩)O。
