c语言酶切算法,科学网—FitHiC V1算法解析(一) - 卢锐的博文

FitHiC V1主要用于识别中程顺式互作

处理Hi-C数据，最自然的分辨率划分方法是基于限制性内切酶切出来的酶切片段，即一个酶切片段为一个最小单位。但是，因为测序深度和基因组上感兴起的size不同，例如有的研究关注TAD结构，有的研究关注Loops结构，所以目前要么按固定大小对基因组划分bin，要么将多个连续的酶切片段(metafragment)合到一起划分bin。FitHiC中将这些分辨率单元称为"locus"。一对locus之间的连接数称为"contact count"，那么基因组上所有的locus pair的contact count构建的矩阵称为"contact map"。FitHiC V1关注的是中程连接("mid-range" contact,例如对酵母而言是10-250kb，对其它复杂真核生物则为50kb-10Mb)，且为顺式互作(intra-chromosomal)。

事实上，因为可能出现随机成环(random looping)，即便是中程顺式互作的检测也是一项非常大的挑战。另外，很多enhancers与promoters之间的连接，或者多个promoters之间的连接都是中程顺式互作。因此，FitHiC V1(发表于2014年1月)采用了循序渐进的策略，先解决相对较容易的中程顺式互作，再在此基础上进行改进和优化，六年后完成了FitHiC V2版(发表于2020年1月)，实现了反式互作以及任意范围顺式互作的连接检测工作(这部分后续会另写博文介绍)。

对中程顺式互作的准确检查，除了要控制随机引物成环效应(random ploymer looping effect)，还得考虑来自实验和技术方面的偏好性，例如GC含量(GC content)、可比对性(mappability)，交联偏好(cross-linking preference)、酶切片段长度(fragment length)，以及环化长度(circularization length)等。

1. discrete binning approach

这个方法是Zhijun Duan于2010提出来的[1]。顺带一提，Zhijun Duan也是DNase Hi-C技术的发明者。下面以反式互作(trans-interaction)为例，讲解其算法思路。

我们将Hi-C进行比对、过滤和校正之后，统计locus pair之间的contact count，仅考虑有contact的locus pair。先定义两个值：

[1] M值。它表示染色体间locus pair的数目。假设随机取一对locus pair，取到的概率是等可能的，即概率为p=1/M。

[2] N值。它表示染色体间locus pair观测到的contact counts的总数。

那么，对于某个locus pair，可以通过二项分布给出不多不少正好k个contacts的概率为：

于是，locus pair至少有k个contacts的概率，为

这个累积概率即p-value值。

可能有的人无法理解，我们这边就用示意图帮助大家认识一下：

注：图中一条虚线表示一个contact

(1)如图，我们首先对染色体Chr1和Chr2按照固定长度划分单元

注意：为了和后面equal-occupancy binning方法中再次划分的bin在名称上进行区分，FitHiC特意将开始划分的单元(分辨率)称为loci，后续所有操作的最小单元是loci。另外，loci的复数形式为locus。

为了简化模型，例子中假设这个物种只有两条染色体，即图中的Chr1和Chr2。Chr1和Chr2上contact count大于0的locus pair有M个，而且这M个locus pair概率相等，其概率为p=1/M。即从Chr1上随机挑选1个loci，从Chr2上随机挑选另一个loci，这两个locus构成locus pair，那么在这些locus pair中contact count大于0的有M对，那么从这M对中再随机挑选1对的概率为1/M

(2)如果总的contact count为N，即我们总的有N对PE reads，这些PE reads一条在Chr1上，一条在Chr2上。且PE reads是随机地连接locus pair，我们可以想像一下，第一对PE reads，R1随机地放到Chr1上一个loci，R2随机地放到Chr2上一个loci；然后是第二对PE reads，R1放到Chr1上随机的一个loci，R2放到Chr2上随机的一个loci。这样将所有的N对PE reads都放到locus pairs中。

(3)在上步中，如果我们只关注某一对locus pair，如上图中loci pair (l1,l2)。那么每一次放一对PE reads，只会出现两种情况，要么放到这对locus pair里，要么不放在这一对里。放到这一对locus pair的概率是p=1/M。

上述这个过程服从二项分布。因此求loci pair (l1,l2)正好放了k对PE reads的概率，即公式(1)。回想一下，大学课本里的例子，抛5次正常的硬币，求正好有3次正面朝上的概率。

理解了上面的过程是二项分布，接着算loci pair (l1,l2)放的PE reads数至少为k的概率，即k对的概率加上k+1对的概率，一直加到N对的概率，为N时表示所有的PE reads都放到loci pair (l1,l2)的概率。总的累积概率即p-value值。

刚才提到的是反式互作的情况，顺式互作会更加复杂一些，因为顺式互作会存在距离效应(distance effect)——同一条染色体上的两个点，在一维基因组上距离越远，在三维构象中其互作越少，且这种越少呈现幂律分布(power law distribution)。考虑到距离效应，Zhijun Duan[1]团队想出来的解决办法是再次划分bin，即locus pair之间的距离按照0-5kb，5kb-10kb，10kb-15kb等等，如下图。因为离得比较近的位点间随机连接比较多，所以0-20kb这种近程互作直接清理掉。可以认为locus pair中距离相等(例如相距45kb-50kb)的这些locus pair，它们的概率是相等的，因此也是服从二项分布。

注：

(1)为了绘图方便，这里假设loci的长度为1.66kb(即5kb/3)，在实际问题中loci的长度会有不同。

(2)图中一条虚线表示一个contact

显然，相距较近的locus pair，其contact更多，即图中红色的连接要比该loci蓝色的连接多；但是相距相差不大的locus pair其连接数相差不大，例如图中两个由红色虚线连接的locus pair之间的contact count相差不大。

因此，可以将相距20-25kb的所有locus pair作为一组单独分析；将相距25-30kb的所有locus pair作为另一组单独分析。

我们将重新划分出来的距离作为bin，这里的bin i表示的是第i个距离范围[si, ei)，其中s表示start，e表示end。例如如果bin 1对应20kb-25kb，则s1为20kb，e1为25kb；如果bin 2对应25kb-30kb，则s2为25kb，e2为30kb。那么对顺式互作重新定义M值和N值，如下：

[1] Mi值，表示距离为[si, ei)范围的locus pair总数

[2] Ni值，表示距离为[si, ei)范围的locus pair观测到的contact counts的总数。

代入公式(1)和(2)也可以计算出p-value值，随后为了获取到更显著的结果，很自然地采用了q-value。

有关p-value、FDR校正和q-value，我另外整理了一篇博文《p-value(P值), FDR以及q-value(Q值)》，值得一提的是该博文的框架是FitHiC通迅作者William S Noble于2009年发表于NB上的一篇专门讲解p-value, FDR和q-value的文章。建议研究课题中涉及p-value和q-value的读者仔细阅读一下文献[3]。

参考材料：

[1] Zhijun Duan, Mirela Andronescu, Kevin Schutz, Sean McIlwain, Yoo Jung Kim, Choli Lee, Jay Shendure, Stanley Fields, C. Anthony Blau, William S. Noble. A three-dimensional model of the yeast genome. 2010

[2]Ferhat Ay, Timothy L. Bailey, William Stafford Noble.Statistical confidence estimation for Hi-C data reveals regulatory chromatin contacts. 2014

[3] William S Noble. How does multiple testing correction work? 2009

转载本文请联系原作者获取授权，同时请注明本文来自卢锐科学网博客。

链接地址：http://blog.sciencenet.cn/blog-2970729-1223469.html

上一篇：遗传算法(genetic algorithm)简介

下一篇：FitHiC V1算法解析(二)