生存分析凭什么不需要矫正P值

生存分析是大数据时代,筛选目标基因的超级有效策略。各种数据挖掘文章本质上都是要把目标基因集缩小,比如表达量矩阵通常是2万多个蛋白编码基因,不管是表达芯片还是RNA-seq测序的,采用何种程度的差异分析,最后都还有成百上千个目标基因。如果是临床队列,通常是会跟生存分析进行交集,或者多个数据集差异结果的交集,比如:多个数据集整合神器-RobustRankAggreg包 ,这样的基因集就是100个以内的数量了,但是仍然有缩小的空间,比如lasso等统计学算法,最后搞成10个左右的基因组成signature即可顺利发表。

虽然生存分析如此重要而且如此常见,但是仍然有一些未解之谜,不同数据库来源,病人的不同时期的记录信息,以及不同的阈值分组,拿到的结果居然是可以不一样的!虽然大家都倾向于做各种花式分析,然后挑选具有统计学显著意义的生存分析结果。

生存分析最重要的是病人分组

我在生信技能树多次分享过生存分析的细节;

可以看到,有基因表达量高低分组,基因突变与否分组,多个基因表达量和突变联合分组,甲基化高低分组,gsea和gsva等基因集得分进行分组,五花八门,其中200块的代码我的学徒免费送给你,GSVA和生存分析 视频值得看!

TCGA数据库的RNA-seq表达矩阵全部基因高低表达分组后批量生存分析

虽然说可以各种花式分析,然后挑选具有统计学显著意义的生存分析结果,但是最开始基本上都是对全部基因分组后批量生存分析,可以是表达量高低,包括mRNA,lncRNA,miRNA的表达量,以及甲基化信号值高低等等,一个基因可以把病人分组,只要是有分组,就可以进行一次生存分析。

比如我们可以下载TCGA数据库的RNA-seq表达矩阵,读入到R里面构建成为 expr 这个数据变量,然后整理好临床表型,构建成为phe这个变量,接下来就可以使用下面的代码对RNA-seq表达矩阵全部基因高低表达分组后批量生存分析:

## 批量生存分析 使用 logrank test 方法
mySurv=with(phe,Surv(time, event))
log_rank_p <- apply(expr , 1 , function(gene){
 # gene=as.numeric(expr[1,])
 phe$group=ifelse(gene>median(gene),'high','low') 
 data.survdiff=survdiff(mySurv~group,data=phe)
 p.val = 1 - pchisq(data.survdiff$chisq, length(data.survdiff$n) - 1)
 #cat(p.val) 
 return(p.val)
})
log_rank_p=sort(log_rank_p)
head(log_rank_p)
boxplot(log_rank_p)

但是最近挑选具有统计学显著意义的生存分析结果的基因时候,发现很多基因都是表达量相关的,也就是说,它们尽管说是不同的基因在不同病人表达量不同,但是它们对病人的分组效果其实是类似的。

那么,我们一下子对几万个基因进行批量生存分析,每一次每一个基因的生存分析都是独立的P值,为什么我们没有对这样的P值进行矫正呢?

大家耳熟能详的矫正P值有,adjust.p,q值,以及FDR,他们的作用都是把P值放大,这样之前那些小于0.05或者0.01的具有统计学显著的基因就不再显著啦,就是把筛选标准严格一点而已。

生存分析凭什么不需要矫正P值?

难道就是因为我们希望统计学显著的生存结果,就选择性展示它吗?

文末友情推荐

要想真正入门生物信息学建议务必购买全套书籍,一点一滴攻克计算机基础知识,书单在:什么,生信入门全套书籍仅需160
如果大家没有时间自行慢慢摸索着学习,可以考虑我们生信技能树官方举办的学习班:

如果你没有服务器的话,做NGS数据分析实战可能会有点勉强,建议考虑:每天不足一块钱,定制生信云送给你

Comments are closed.