这是看《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】

（未完待续）