生物信息学习的正确姿势

NGS系列文章包括NGS基础、高颜值在线绘图和分析、转录组分析 (Nature重磅综述|关于RNA-seq你想知道的全在这)、ChIP-seq分析 (ChIP-seq基本分析流程)、单细胞测序分析 (重磅综述:三万字长文读懂单细胞RNA测序分析的最佳实践教程)、DNA甲基化分析、重测序分析、GEO数据挖掘(典型医学设计实验GEO数据分析 (step-by-step))、批次效应处理等内容。

姊妹篇:总被审稿人提起的多重假设检验校正是什么?

前几天,Nature上一篇comment再度引发关于p-value如何使用和解释的文章:Scientists rise up against statistical significance,800多名科学家联合声明拒绝使用基于p-value置信区间贝叶斯因子等的二分法将研究结果分为统计显著和统计不显著两个部分,而是应该把置信区间改为兼容性区间, 描述区间所有值的实际含义,尤其是其所代表的的效果 (point estimate)或极值在哪。给定了统计假设,任何极值内的值与研究数据都是兼容的。基于此,作者可以更好的强调数据分析带来的期望值和不确定性,不再对结果过于自信或悲观。

不过一来统计界以后会怎么实施未知,二来签名也未发对p-value的正确使用。那么怎么理解P-value的含义?怎么算是正确使用P-value呢?怎么评估算出的P-value是否正常呢?就是我们下面要说的。基于传统,后面还是会继续使用显著性这一说法。

统计分析检验获取p-value是我们经常要做的一个工作,比如获得差异基因或富集分析等。通常计算后会得到数百、数千或数万个p-value。考虑到多重假设检验的问题,你可能会想着先做一个校正。

然而,你最先需要做的却是绘制一个直方图。怎么绘制?简单强大的在线绘图-第3版。

在做任何的多重假设检验校正、假阳性率控制或结果解释之前,先绘制这么一个p-value分布直方图,它可以告诉你在所有假设的p值分布,并帮您发现潜在的问题。

p-value分布直方图可能有下面6种可能,我们一一看来。

Anti-conservative p-value

如果p-value分布直方图如上图所示,左侧0值附近有个峰,右侧为近乎均匀分布,那么恭喜你,这是一个很好的分布。

0-1之间均匀分布的p-value代表原假设H0 (null hypothesis)的P值。为什么它们是均匀分布的呢?这是根据p-value的定义来的。在原假设下,p-value有5%的可能低于0.05, 10%的可能低于0.1,以此类推,就是一个均匀分布。

p-value接近于0值的峰代表的是备择假设H1 (alternative hypothesis) (也包含部分假阳性)。如果把原假设和备择假设分开,p-value的分布应该入下图所示:

首先可以看到在低p-value处也有一些原假设 (H0),因此不可以简单的说所有p-value<0.05的都是显著的,否则就会获得一些假阳性结果。而且一些备择假设 (H1)的p-value也比较高,这些就是不能通过本次统计检验方法获得的阳性结果,也称为假阴性结果。

多重假设检验校正就是确定显著性的合理阈值。

那么怎么判断多少假设是原假设,多少是可以拒绝原假设采用备择假设呢?可以从下面几张图有个直观认识,左侧Peak越高,越多的假设p-value趋近于0, 也就是显著的结果。右侧的柱子越高,更多原假设不能被拒绝。如果想获得定量的评估,可以使用qvalue包。

library(qvalue)
data(hedenfalk)
pvalues <- hedenfalk$p
qobj <- qvalue(p=pvalues)
summary(qobj)

输出不同p-value假设的累计数目

Call:
qvalue(p = pvalues)pi0:    0.669926Cumulative number of significant calls:<1e-04 <0.001 <0.01 <0.025 <0.05 <0.1   <1
p-value       15     76   265    424   605  868 3170
q-value        0      0     1     73   162  319 3170
local FDR      0      0     3     30    85  167 2241

估计原假设 (H0 null hypothesis)的整体比例 (π0),q-valuep-value的关系, qvalue即是定义某一个检验统计显著需要承受的最小假阳性率值。lfdr指在给定的p-value条件下,原假设 (H0)为真的后验概率值。

hist(obj)

均匀分布 Uniform p-value

假如,你的p-value是如下图所示,平平的均匀分布,怎么办呢?

看上去所有的假设都符合原假设,是不是意味着就没有办法拒绝原假设了?其实也不是:

  • 起码有一小部分的假设是备择假设,可以用过FDR校正方法如Benjamini-Hochber等鉴定出来。

  • 直接应用p-value<0.05是不合适的,假阳性率会很高。

双峰 Bimodal p-values

如前面所示在p-value=0处有一个峰,但在p-value=1处也有一个?怎么解释。

首先不要对这些p-value应用假阳性率控制。为什么呢?因为一部分FDR控制算法是基于P-value1附近是均匀分布的。如果不符合这个前提,计算出的显著性会很少。

下一步找出为什么p-value会有这个分布,针对性解决:

  • 是否使用的是单端检验 (one-tailed test) (如检验药物处理后基因表达上调)。如果是这样,p-value接近1的正好是相反的变化 (如基因表达下调)。如果您同时关注上下调,则采用双端检验 (two-sided test)。如果您不想包含另一种变化,则在检验前先过滤掉这些。(注:比如富集分析时只关注富集)

  • 是否pvalue接近1的情况都是病态情况,如基因差异表达分析中,一些软件会赋予在所有样品中都不表达的基因检验pvalue1,这样的情况直接过滤掉就好。(注:一般分析时是提前过滤。)

Conservative p-values

看到这个分布,不要鲁莽的下结论:没有任何统计显著的假设。如果真的没有统计显著性假设,p-value的分布应该是均匀的 Uniform, 这是因为p-value就是这么定义的:原假设下均匀分布。

如果p-value呈现这个分布,说明统计检验使用错了。其原因可能是数据的分布不符合统计检验的假设,比如统计检验适用于连续数据,而提供的是离散数据,或者统计检验适用于正态分布数据,而提供的数据严重不符合等。最好的解决办法是找一个友好的统计学家朋友帮助您。

我们一直强调可视化的是原始p-value的分布,如果使用的工具不小心提供的是校正后的p-value,比如使用Bonferroni correction,那么校正后的p-value可能是这个分布。

稀疏分布 Sparse p-values

如图所示,获得的p-value的值比较单一,假如做了10,000次统计检验,只获得很少的不同的检验p-value,可以使用下面的代码获取总共有多少不同的p-value

length(unique(mypvalues))

为什么会获得这样的p-value呢?

  • 自展或置换检验 (bootstrap or permutation test)的迭代次数太少。

  • 数据集小的时候运行了非参数检验 (如Wilcoxon rank-sum testSpearman correlation),尝试扩大样本量或数据转换为可以进行参数检验。

不要做假阳性率控制,因为p-value的分布不是连续的。

悟空庙宇P-value (“What the…?!?”)

像不像孙悟空变的一座庙,尾巴做旗杆?中间的P-value有个凸起,在1附近有个峰。

最好的方式是求助于统计学家,当然在这之前,看下数据的分布,了解下所用的统计方法,先有个直观认识。

所以p-value不是算出来就可以用了,观察其分布,可以帮助我们判断数据分布是否合适,选用的统计检验方法是否合适,后期如何进行处理,对结果解释增强可信度。

参考

  1. http://varianceexplained.org/statistics/interpreting-pvalue-histogram/

  2. http://www.bioconductor.org/packages/release/bioc/vignettes/qvalue/inst/doc/qvalue.pdf

  3. https://www.nature.com/articles/d41586-019-00857-9

  4. https://stats.stackexchange.com/questions/10613/why-are-p-values-uniformly-distributed-under-the-null-hypothesis#

GEO/TCGA数据

  • UCSC XENA - 集大成者(TCGA, ICGC)

  • ICGC数据库使用

  • TCGA数据库在线使用

  • BROAD开发的TCGA分析平台,强大的下载功能

  • cBioPortal功能强大的TCGA再分析平台

  • 这是数据更新最实时的TCGA网站,功能强大

  • 不懂R,如何进行GEO数据库表达谱的差异分析、富集分析、蛋白互作、可视化?

  • 典型医学设计实验GEO数据分析 (step-by-step) - Limma差异分析、火山图、功能富集

  • 典型医学设计实验GEO数据分析 (step-by-step) - 数据获取到标准化

往期精品(点击图片直达文字对应教程)

后台回复“生信宝典福利第一波”或点击阅读原文获取教程合集

什么,你算出的P-value看上去像齐天大圣变的庙?相关推荐

  1. c语言根据日期求星期不用公式,...迅速算出是星期几的方法给你任何一年看日期怎么能...-知道日期怎么算周几-数学-倪滩贝同学...

    概述:本道作业题是倪滩贝同学的课后练习,分享的知识点是知道日期怎么算周几,指导老师为冉老师,涉及到的知识点涵盖:...迅速算出是星期几的方法给你任何一年看日期怎么能...-知道日期怎么算周几-数学,下 ...

  2. AI算法 真的能算出人类的欲望吗?

    https://www.toutiao.com/a6686069891713221132/ 假期里,闲着也是闲着,我们聊点烧脑的话题. 以前,我给一家朋友的机器人公司写过一个<AI时代的爱情&g ...

  3. 拾谈“用最有效率的方法算出2乘以8等於几?”

    这是网上流传的"变态级JAVA程序员面试32问"的其中一题(二十八题),然后下面给出来的答案是 第二十八,编程题: 用最有效率的方法算出2乘以8等於几? 有C背景的程序员特别喜欢问 ...

  4. 2张图片就能「算出」一段视频,Reddit网友都惊呆了 | 旷视北大出品

    鱼羊 萧箫 发自 凹非寺 量子位 报道 | 公众号 QbitAI 只给AI两张图片,就能得到高帧率动态视频? 输入的两张图像,重叠后是这样的: 而算出来的视频,是酱婶的: 不错,这又是视频插帧算法的功 ...

  5. 30岁程序员吐槽:一分钟只能赚3.3元,混得太差!算出月薪后我服了

    现在很多年轻人都会在年龄还很小的时候,就开始思考着怎么去找一份比较好的工作.即使是进入职场工作,拿着一份整体水平比较低的薪资,也会像办法让自己的薪资待遇提高一些.毕竟谁不愿意让自己拿到更多的工资 在一 ...

  6. JAVA班级年龄平均值代码_java用list集合存储学生信息并算出成绩平均值操作

    需求 键盘输入五名学生信息并录入list集合; 输出每个学生的信息,计算并输出这五个学生Java语言成绩的平均值: 计算并输出他们Java语言成绩的最大值和最小值. 思路 用Scanner 键盘输入 ...

  7. java 学生信息 list_java用list集合存储学生信息并算出成绩平均值操作

    需求 键盘输入五名学生信息并录入list集合; 输出每个学生的信息,计算并输出这五个学生Java语言成绩的平均值: 计算并输出他们Java语言成绩的最大值和最小值. 思路 用Scanner 键盘输入 ...

  8. android double值排序,android根据Double类型数据经纬度算出距离再根据距离实现排序功能...

    前言 项目中用到全国的加油站数据加载 并根据经纬度算出距离 然后根据距离从小到大排序 主要是数据类型是Double 这里必须对数据进行封装 实现也不难 这里讲一下自己的实现方法和实现思路 效果图 先来 ...

  9. 如何用计算机算出男朋友的身高,【趣味物理】如何用物理方法测出男生的真实身高?...

    原标题:[趣味物理]如何用物理方法测出男生的真实身高? 毕导,本名毕啸天,清华大学化工系博士生,今日头条优质科普内容创作者. 如何测出男生的真实身高? 男生常常会隐瞒身高,173敢报178,175就敢 ...

最新文章

  1. 我用2年时间从财务到数据分析师!
  2. AI大牛周明打造的轻量“孟子模型”开源!靠10亿参数冲上CLUE榜第三,可用于新闻分类、文案生成...
  3. LODOP打印table表格宽度固定-超宽隐藏
  4. Javascript中数组去重的六种方法
  5. spring—事务控制
  6. 前端学习(3123):react-hello-react之props的基本使用
  7. 基于lucene语法的实时文本搜索与匹配--Tripod
  8. Eclipse 启动不了 Tomcat
  9. Kotlin从入门到放弃(三)——协程
  10. beautifulsoup_Py之Beautiful Soup:Beautiful Soup 4.2.0的简介、安装、使用方法
  11. 【java】信号量机制
  12. 树梅派应用27:通过USB蓝牙适配器连接BLE设备
  13. Android 点击键盘外 非输入框 关闭软键盘
  14. 计算基因上外显子碱基覆盖度(exon coverage depth):Samtool工具使用
  15. 修改服务器bond网口mode4,双25GE网卡做bond4测试,其中一个网口没有流量一个网口可以打满的问题分享★★★...
  16. redis单点故障问题
  17. 论文笔记| The Emergence, Advancement and Future of Textual Answer Triggering
  18. alias linux 执行命令,Linux系统alias命令编写实现命令别名方法介绍
  19. 甜品消消乐知识点总结
  20. 李梦娇口诀88条(视频+讲义)

热门文章

  1. 杰理之JLANC开发工具使用说明【篇】
  2. 水星路由器 Mercury MER1200G刷机教程,Archer C5V4刷回原厂固件,串口(TTL)刷机,需要一个USB转TTL工具
  3. win32 IOCTL学习总结
  4. 激活office python
  5. 《JavaScript函数式编程思想》——递归
  6. 张一鸣:10年面试2000人,我发现混的好的人,全都有这5种特质
  7. 如何把acer aspire5560G笔记本拆开清理灰尘
  8. 培训报名小程序实战开发
  9. 利用selenium保存静态网页
  10. mysql 的 infobright 数据库的 mediumblob 显示不了数据