ChIP-seq文章
8年前的ChIP-seq如何找peak
老板前几天丢了2010年的文章给我，让我去看下这篇文章的ChIP-seq分析结果中，他要的基因是不是也被调控了。于是我就看到了8年的ChIP-seq分析中是如何寻找Peak region。
数据来自于Illumina的单端测序，长度为42bp。分别2个ChIP-seq独立的生物学重复和2个对照生物学重复，但是一共有7份数据。因为第一组ChIP-seq还有3个技术重复，第一组对照还有2个技术重复，所以是3+1+2+1=7 。他们的ChIP-seq样本在测序之前用qPCR先进行了验证。
由于当时测序技术还一般，后面几个bp质量就急剧变差，于是先用了GenomeMapper直接裁剪到36bp，数据比对用的是SHORE，允许4个碱基错配，不允许gap。目前质量控制可用trimmomatic，比对可用bowtie1，接着用samtools转换格式和排序。这些都比较常规，不需要过多介绍，下面找到差异需要仔细讲讲。
作者的策略分为两步：第一步用slide window在ChIP-seq样本中找到潜在可能的富集区域。第二步对该区域进行定量，通过和对照组进行比较，从而找到统计学上显著富集的区域。
step1: 找到ChIP样本中潜在区域
为了实现第一步，首先需要用唯一匹配的read（uniquely mapping reads）去计算基因组的片段覆盖图（fragment coverage graph），也就是说在基因组上每个碱基要注释read覆盖深度(coverage depath)。
注意：作者说“如果你懂一点ChIP实验的话，那你就知道有一半的ChIP实验表明蛋白和DNA免疫沉淀的片段平均长度是200～300bp“，所以作者先把原来只有34bp的read覆盖在3‘方向拓展了130bp，假装自己测序是SE130. 不过现在基本都是PE100，PE150，你都不需要假装了。
由于我对ChIP这方面实验还不太懂，于是我去问了公众号《嘉因》的小编—小丫，她说如果是实验的打断效果是在300-700，如果是建库的size selection则是200bp。
得到全基因组每个碱基的read的覆盖深度后，作者使用2kb的window，以1bp为步长开始扫描基因组检测peak region。每一步都需要评估潜在的富集中心碱基（potential enrichment of the central base）。评估方法是单边泊松检验，其中lambda设为局部的平均深度（即所有碱基的覆盖深度除以碱基总数）。没有覆盖的地方就是ChIP-seq没有结合的区域，所以不计算哦。这一步得到的P<0.05的区域即为P1.
然后继续从P1里挑选出那些连续130个碱基的P值都低于0.05的部分作为P2.
P2区域继续和对照组进行比较，如果对照组的那些区域也是很高表达，那么这些区域继续被剔除，留下的就是P3. 标准是P2潜在的peak region对应的对照区域的平均覆盖度(mean coverage) 大于正处理的中位数+3个标准差（median coverage plus three deviations）
step2: 找到显著性富集的区域
这部分首先使用单边二项检验， 参数之一的N为ChIP和对照实验中比对到潜在区域到read总数。参数之二的成功概率‘r'，作者认为r = s/(s+1). 其中s为对照的比例因子(scaling factor).s的计算方式为: 将整条基因组切割为400bp大小的分箱(bin)，每个分箱分别记录处理样本和对照样本的比对reads数，s就是那个在这些分箱中让处理样本的read count和对照样本的read count的系数。
上面得到的一些列p值还需要用 Benjamini-Hochberg矫正方法转换为FDR。为了进一步提高peak区域的准确，作者还计算了per base excess，也就是正处理和对照相比每个碱基额外的覆盖度的均值chip reads – (s * control reads)]/ peak width
最后作者选择了两个独立重复中，FDR < 10^-10 和 per base excess > 0.25的peak region。至于为啥选择那样子的FDR，作者给出的原因是二项分布的模型不太好，所以要非常严格才行。
if the variance of the true distribution is larger. Alternatively, it may be that AP2 binds an unusually large number of regions in the genome at a lower affinity. In any case, we chose to err on the side of higher peak calling stringency.
