生存分析凭什么不需要矫正P值
生存分析是大数据时代筛选目标基因的超级有效策略。各种数据挖掘文章本质上都是要把目标基因集缩小,比如表达量矩阵通常是2万多个蛋白编码基因,不管是表达芯片还是RNA-seq测序的,采用何种程度的差异分析,最后都还有成百上千个目标基因。如果是临床队列,通常是会跟生存分析进行交集,或者多个数据集差异结果的交集,比如:多个数据集整合神器-RobustRankAggreg包 ,这样的基因集就是100个以内的数量了,但是仍然有缩小的空间,比如lasso等统计学算法,最后搞成10个左右的基因组成signature即可顺利发表。
虽然生存分析如此重要而且如此常见,但是仍然有一些未解之谜,不同数据库来源,病人的不同时期的记录信息,以及不同的阈值分组,拿到的结果居然是可以不一样的!虽然大家都倾向于做各种花式分析,然后挑选具有统计学显著意义的生存分析结果。
生存分析最重要的是病人分组
我在生信技能树多次分享过生存分析的细节;
基因表达量高低分组的cox和连续变量cox回归计算的HR值差异太大? 学徒作业-两个基因突变联合看生存效应 TCGA数据库里面你的基因生存分析不显著那就TMA吧 对“不同数据来源的生存分析比较”的补充说明 批量cox生存分析结果也可以火山图可视化 既然可以看感兴趣基因的生存情况,当然就可以批量做完全部基因的生存分析 多测试几个数据集生存效应应该是可以找到统计学显著的! 我不相信kmplot这个网页工具的结果(生存分析免费做) 为什么不用TCGA数据库来看感兴趣基因的生存情况 200块的代码我的学徒免费送给你,GSVA和生存分析 集思广益-生存分析可以随心所欲根据表达量分组吗 生存分析时间点问题 寻找生存分析的最佳基因表达分组阈值 apply家族函数和for循环还是有区别的(批量生存分析出图bug) TCGA数据库生存分析的网页工具哪家强 KM生存曲线经logRNA检验后也可以计算HR值
可以看到,有基因表达量高低分组,基因突变与否分组,多个基因表达量和突变联合分组,甲基化高低分组,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值?
难道就是因为我们希望统计学显著的生存结果,就选择性展示它吗?