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

博文

FitHiC V1算法解析(一)

已有 5349 次阅读 2020-3-14 14:06 |个人分类:Hi-C|系统分类:科研笔记| Fit-Hi-C, 显著互作, Hi-C, Loops

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

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

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

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

1. discrete binning approach

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

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

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

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

那么,对于某个locus pair,可以通过二项分布给出不多不少正好kcontacts的概率为:

于是,locus pair至少kcontacts的概率,为

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

 

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

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

(1) 如图,我们首先对染色体Chr1Chr2按照固定长度划分单元

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

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

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

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

上述这个过程服从二项分布。因此求loci pair (l1, l2)正好放了kPE 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-5kb5kb-10kb10kb-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表示starte表示end。例如如果bin 1对应20kb-25kb,则s120kbe125kb;如果bin 2对应25kb-30kb,则s225kbe230kb。那么对顺式互作重新定义M值和N值,如下:

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

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

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

 

有关p-valueFDR校正和q-value,我另外整理了一篇博文《p-value(P), FDR以及q-value(Q)》,值得一提的是该博文的框架是FitHiC通迅作者William S Noble2009年发表于NB上的一篇专门讲解p-value, FDRq-value的文章。建议研究课题中涉及p-valueq-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




https://blog.sciencenet.cn/blog-2970729-1223469.html

上一篇:遗传算法(genetic algorithm)简介
下一篇:FitHiC V1算法解析(二)
收藏 IP: 113.57.182.*| 热度|

0

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

数据加载中...

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

GMT+8, 2024-11-23 07:05

Powered by ScienceNet.cn

Copyright © 2007- 中国科学报社

返回顶部