Songjinghe的个人博客分享 http://blog.sciencenet.cn/u/Songjinghe

博文

生物信息学中多重检验问题的一个例子(FDR、p值、q值)

已有 26372 次阅读 2016-4-5 17:57 |个人分类:统计|系统分类:科研笔记| 生物信息学, 多重检验, FDR, q-value

这是看《How does multiple testing correction work》一文的笔记和总结。


多重检验问题是搞生物信息学时几乎绕不开的问题,笔者不是生物信息学出身(学计算机的),被FDR,p-value,q-value搞得晕头转向,听人家说这篇文章讲的比较清楚,所以耐着性子看了一遍来给大家汇报。但是注意这只是一个比较直观的解释,如果想详细了解FDR的统计学意义请参考其他文章


首先它这个文章的思路有点怪,我想把它套在数理统计这个课程讲的思路上总是套不上,理了半天才理顺了。


这种东西肯定是要举个例子才能讲清楚的,作者举了一个例子是在人类21号染色体上检测CTCF基因(大家自己百度,好像是转录阻抑啥啥的一个基因吧,然而就算不理解也不太妨碍理解这篇文章)【To better understand this protein, you want to identify candidate CTCF binding sites in human chromosome 21】。方法是和某个论文提出的模型进行匹配,说白了DNA不就是一个只由ATGC四种字母组成的超级长的字符串嘛(这里的长度是6800万),然后从里面随机抽连续的20个字母拿出来输入一个函数,这个函数会输出这20个字符串和CTCF匹配程度的得分(-100到20不等的实数,得分大于17好像就算很好的匹配了——在论文里)。然后得分越高就说明越是这个基因。



论文里面讲的是通用的生物信息学手段,但是和一般的假设检验套不上,在本文中我会试着用假设检验的思路讲。



准备知识1:【null model】这个东西实际上就是H0为真时总体的分布,有两种方法可以得到:一是可以从原始数据出发构造【empirical null model】,方法有permutation test,bootstrap,文章里用的是把原始数据乱序排列后做抽样,可以认为是非参数方法。二是认为原始数据是按照一定规律“生成”的,这样可以用数学来表达,叫【analytical null model】,可以认为是参数法。


说明:下文中“观测”==“抽样”,虽然在统计上这两个意思好像不一样。



上帝视角:这是在看过整个例子后把一些东西放到前面讲,免得混了。就是在人类21号染色体上有68000000个字母(核苷酸),其中连续的20个字母得分超过17的有519个位置。如果把整个染色体的68000000个字母乱序排列,那么其中连续的20个字母得分超过17的有35个位置。先别管这些数是咋来的——上帝视角!后面会讲~



此时我们假设一个研究员,他目前只知道整个染色体的DNA上有68000000个字母,其他的都不知道。他要解决的问题是:人类21号染色体上有没有CTCF基因?如果有,(位置)在哪里?于是他做了一个假设检验。


H0:DNA上没有这个基因(整个DNA上没有得分超过17的片段)


H1:DNA上有这个基因(整个DNA上有得分超过17的片段)


然后他从整个DNA上随机抽取了一个位置(比如19390631),拿出连续的20个字母,一算得分——【标记A,待会还会回到这里】——26分!嗯,他感觉不错,看下显著性怎么样吧,于是他把整个DNA上的68000000个字母乱序排列(构造empirical null model),随机抽取连续的20个字母算得分,抽了1000次,没有得分超过17的。他很高兴,看来【由于随机的ATGC排序使计算得分大于17】是一个小概率事件,发生概率(即p值)最多也不超过1/1000=0.001,根据【小概率事件在一次观测中几乎不可能出现】的统计学原理,可以拒绝H0,认为DNA上有这个基因。于是他高兴的在实验报告上写到:在19390631处抽样一次,得分26,拒绝原假设,p<0.001。



好了,回到上帝视角,我们发现这真是一个幸运的研究员,简直可以去买彩票了!为什么呢?因为整个原始数据中连续的20个字母得分超过17的只有519个位置,在一次观测中就抽到这个位置的概率是519/68000000=0.000076。这个概率太小了,按理说至少要抽样1/0.000076=13157次才能抽到这么一个大于20的得分,而这个研究员居然一次抽样就抽到了!这让大家羡慕嫉妒恨啊,不行,我们得把他的运气调成和一般人一样,然后再看看会发生什么。于是回到【标记A】



——0分!嗯运气不是很好,继续抽样!——又抽样了1000次,居然没有一次得分是大于17的!研究员很郁闷——关键是无法得出任何结论,因为可能还真的就有那么几个位置得分大于17~但是,简直是大海捞针啊,这么抽样得抽到啥时候去!怒了!于是写了一个自动抽样的程序,从整个DNA头上开始抽样,比如第一次抽1~20个字母,第二次抽2~21个字母,依次类推,一直抽样到67999980~68000000,每次拿出来算一个分,看看到底有没有得分超过17的!这是把整个样本空间都遍历了一遍。哈哈~



过了很久,结果终于出来了,居然真的有519个位置的得分超过17!太激动了——赶紧写在实验报告上——等下,还没算显著性呢~仍然是构造empirical null model,随机抽取连续的20个字母算得分,抽了1000次,没有得分超过17的。他很高兴,看来【由于随机的ATGC排序使计算得分大于17】是一个小概率事件,发生概率(即p值)最多也不超过1/1000=0.001,根据【小概率事件在一次观测中几乎不可能出现】的统计学原理,可以拒绝H0,认为DNA上有这个基因——等下,好像不对劲,我们之前在原始数据中抽样了1000次也没有抽到得分大于17的啊?不对不对,关键是【小概率事件在一次观测中几乎不可能出现】这条统计学原理不能用在这里,因为我们得到519个位置是进行了67999980次抽样(而不是一次观测)的结果,这么大量的抽样,那就应该用另一个统计学原理【小概率事件在大量观测中几乎必然出现】,那么也就是说在519个位置中,肯定有【由于随机的ATGC排序使计算得分大于17】的情况。具体有多少呢?根据p值的定义,如果H0成立,p=0.001的意思是观测(抽样)1000次会产生一个小概率事件,那抽样67999980次就会产生约68000个小概率事件,即在所有得分大于17的位置中有68000个是【由于随机的ATGC排序使计算得分大于17】的!而我总共才观测到519个!——当然不用激动,我们的p值是小于0.001的,具体是多少还不知道。这里只是先从逻辑上推一下。


所以现在的情况是,虽然我找到了519个得分大于17的位置,但是仍然不能拒绝原假设(H0),因为不知道到底有多少个观测值是由于大量抽样产生的小概率事件【由于随机的ATGC排序使计算得分大于17】。故当前要紧的是算出p的准确值来。就是H0成立时【由于随机的ATGC排序使计算得分大于17】这个事件的概率。怎么算呢?硬算吧,把刚才的那个自动抽样程序拿过来然后在null model上跑一遍。结果发现有35个得分大于17的。嗯,看来519个得分超过17的位置里有35个是偶然因素导致的,那还有519-35=484个是随机因素不能解释的,所以我们可以拒绝原假设,认为DNA序列上有CTCF匹配,FDR=35/519=0.067,所以最终实验报告上的结论是抽样67999980次,得分大于17的519个,以FDR=0.067拒绝原假设。


另一个需要注意的地方是虽然知道有35个是False Positive/Discovery,但是并不知道具体是哪个。


理解了FDR,q-value就好理解了,由于FDR不具有单调性,所以q-value是“单调的FDR”。


其实看了这篇文章我仍然对于多重检验是否适用于我的问题比较晕(在我自己的那个实验中到底是不是要用FDR这个东西。我这里的实际情况是不同细菌之间算相关系数R,然后也有一个显著性水平)于是又看了另一篇文章《Tutorial in biostatistics:multiple hypothesis testing in genomics》,这个文章讲得比较深入。其中令我印象比较深的一点就是:多重检验问题其实广泛存在!只要你把多个检验放到一起,就有这个问题,哪怕这些检验之间完全不搭边(不同领域,不同问题,不同处理方法,不同样本量)只要是用假设检验来得出一个结论,就可以放在一起,就存在多重检验问题。顺着这个思路举个例子,也就是意味着在全世界所有使用假设检验的论文中(每个论文看做一次假设检验),肯定有一些得出的结论是有问题的!当然我们不知道是哪些论文有问题。因为【小概率事件在大量观测中几乎必然出现】。


(全文完)




【另一种思路,第一版草稿】


论文里面讲的是通用的生物信息学手段,但是和一般的假设检验套不上,我下面先用我自己的思路(假设检验)讲,最后再拐回论文上去:)



最初的目标是“我想知道这个染色体上有没有这一段基因”。那方法是从整个DNA序列上抽样检测(注意不是把所有的位置都检测一遍——这个是搞计算机的思路,从统计学的观点看就太简单粗暴了,我们先假设这么搞不行),抽一个样本(注意只有一个),然后做一次假设检验。这里H0是DNA上没有这个基因,H1是DNA上有这个基因。显著性水平设置成0.05(我随便设的,其实更小也可以下面会讲)。具体的操作方法就是从DNA序列中随机找一个位置,抽取连续的20个字母,输入那个打分函数看一下,如果大于17,那说明肯定有这段基因。为啥呢?



这就又扯到了【empirical null model】上了(本段内容与论文一致),这个东西实际上就是H0为真时总体(所有位置的得分情况)的分布。构造【empirical null model】的方法有很多中,什么permutation test啥的(当时我也是晕头转向,数理统计里面没讲这个呀!)bootstrap(不知道对不对),文章里用的是把原始数据乱序排列后做抽样(也是68000000次,差不多是全集了)找分布。总之我们得到了H0为真时总体的分布,然后一看里面只有35个得分大于17的,就是如果DNA上没有这个基因(H0成立),我观测(抽样)一次,计算得分大于17的概率是35/68000000=0.0000051(这个就是得分大于17的p值)——看来【由于随机的ATGC排序使计算得分大于17】是一个小概率事件。根据【小概率事件在一次观测中不可能出现】的统计学原理,拒绝H0,认为DNA上有这个基因。



但是事情没有这么简单:实际情况是一次抽样几乎肯定抽不到得分大于17的基因序列(上面那个抽到了是我虚构的)。为啥呢?这时候就要上那个简单粗暴的方法了,从头开始,把DNA所有位置都抽出来算一个值看看。于是进行68000000次抽样(全集),然后发现了519个得分大于17的位置。519/68000000=0.000076,也就是抽样一次就得到得分大于17的情况的概率是0.000076,用1/0.000076=13157即至少抽样13157次才可能得到一个得分大于17的位置!



所以说这个时候才真正回到论文的思路上,要解决的问题是:找到所有CTCF基因的位置,然后提取这些位置的基因做进一步的处理。此时的H0仍为DNA上没有这个基因,H1是DNA上有这个基因,但是由于这个基因其实很少,所以必须进行大量的抽样(理论上至少抽样13157次,实际进行了68000000次)才能找到这些基因。但是这样带来的问题是——就算找到了一些位置满足得分大于17,但仍然不能拒绝H0,因为上面拒绝原假设的统计学原理【小概率事件在一次观测中不可能出现】不成立了!不成立的原因是我们进行了大量(68000000次)而不是“一次”观测,统计学还有一条原理叫【小概率事件在大量观测中必然出现】所以找到的得分大于17的位置中肯定有一些是【由于随机的ATGC排序使计算得分大于17】


(未完待续)



https://wap.sciencenet.cn/blog-2827057-968037.html

上一篇:凝聚法层次聚类之ward linkage method
下一篇:茶道的形而上学:意念能让茶更好喝?
收藏 IP: 123.124.23.*| 热度|

2 谢平 xiyouxiyou

该博文允许注册用户评论 请点击登录 评论 (0 个评论)

数据加载中...

Archiver|手机版|科学网 ( 京ICP备07017567号-12 )

GMT+8, 2024-3-28 23:02

Powered by ScienceNet.cn

Copyright © 2007- 中国科学报社

返回顶部