文献阅读 | TopDom: an efficient and

2021-07-27 本文已影响0人阿狸的窝

原文：
Shin H, Shi Y, Dai C, Tjong H, Gong K, Alber F, Zhou XJ. TopDom: an efficient and deterministic method for identifying topological domains in genomes. Nucleic Acids Res. 2016 Apr 20;44(7):e70. doi: 10.1093/nar/gkv1505. Epub 2015 Dec 23. PMID: 26704975; PMCID: PMC4838359.

软件地址：
http://zhoulab.usc.edu/TopDom/.

概要

文章的作者团队开发了一种从Hi-C contact matrices中鉴定TAD的算法TopDom。与之前的算法相比，TopDom不需要复杂的参数调节（parameter tunning），用户仅需要提供window size一个参数。并且与DI相比，TopDom鉴定的TAD大小更小，数目更多。
作者将TopDom应用于2人2鼠4种不同细胞系，再次证明了TAD在不同细胞类型间保守，并且作者观察到管家基因更倾向于分布在这些在细胞类型间保守的TAD边界附近。

背景

Toplogical domain 特点：

同一个TD内的染色质性质相似（组蛋白修饰、活性基因密度、复制时间等）
在不同物种、同一物种的不同细胞类型间保守

已有的domain识别算法

发表年份	作者	算法原理	参考文献
2012	Sexton et al.	distance-scaling factor	[1]
2012	Dixon et al.	DI + HMM	[2]
2012	Hou et al.	Bayesian probability model	[3]
2014	Filippova et al.	Armatus	[4]
2014	Levi-Leduc et al.	HicSeg	[5]
2014	Rao et al.	Arrowhead	[6]

已有算法主要存在的问题

计算耗时
需要调节的参数过多，且参数的选择对结果影响大
不同方法间计算得到的TAD一致性差

TopDom 算法原理

概述

TopDom分为3个步骤：

计算binSignal Value
binSignal的局部最小值作为TAD boundary
统计学检验过滤false detection

计算binSignal Value

给定一个window size $w$ （用户指定参数）
对于每一个bin $i$ ，计算位于其上游的 $w$ 个bin和位于其下游的 $w$ 个bin之间两两的互作频率的平均值
$binSignal(i) = {1\over w^2} \sum_{l=1}^w\sum_{m=1}^w C_{U_i(l), D_i(m)}$
计算线性基因组上每条染色体上的所有bin的binSignal值，可以得到binSignal曲线

截屏2021-07-27 下午3.30.14.png

寻找局部最小值（local minima)

因为TAD内的互作频率应当远高于TAD间互作，所以作者预计TAD中心所在bin的binSignal应当更高，而TD boundary所在bin的binSigna较低。根据这一假设，作者提出，binSignal曲线的局部最低点应当是TD的boundary所在。

但是一个重要的问题是，某些局部最低点可能体现的是数据的噪音。作者希望对binSignal曲线进行平滑，并使用一系列的折线段对来替代原曲线。

在此，作者采用了Kumar Ray et al.的Piecewise linear curve fitting算法[7-8]。

伪代码：

F = F_last = 0
bin_start = start
j = bin_start + 2
while j <= end
    L = Length(bin_start, bin_j)
    E = FittingError() 
    # Fitting error: the sum of distances from the points to the fitted line segments. 
    F = L - E
    if F < F_last:
       # Use bin (j-1) as the end of the last segment and also the start of the next segment
        SaveTurningPoint(j-1)
        bin_start = j - 1
        j = bin_start  + 2 
        F_last = 0
    else:
        # Try if a longer segment will be better
        F_last = F
        j = j + 1

由此方法找出的转折点中的local minima被视作 TD boundary-like bins

false positive 过滤

对于一个bin $i$ ，
定义位于 $i$ 上游的 $w$ 个bin与位于 $i$ 下游的 $w$ 个bin之间的互作为 between interaction
定义位于 $i$ 上游的 $w$ 个bin内部的interaction和位于 $i$ 下游的 $w$ 个bin之间的互作为 within interactions

对于染色体上的每一个bin，使用Wilcox Rank Sum test 检验该bin的 between interactions的normalized contact frequency（校正距离) 是否显著小于 within interactions的normalized contact frequency

使用 Wilcox p-value < 0.05 对 local minima 进行过滤。

local minumum → TD

对于由两个相邻的local minima确定的一段region：

如果region内所有bin都满足 P<0.05，则认为该region为boundary region（或TD-free chromatin region）
否则认为 region 为 TD

注：作者指出，绝大多数local minima都满足 p<0.05，因此在此步仅有很少的bin被过滤

如何选择 window size / TD质量评估

随着window size $w$ 的增大，鉴定得到的TD的size增大而数目减少。

Figure 2 | TD 大小与数目随window size变化的关系

为了选择一个恰当的window size，首先需要一个找到一个指标用于评估计算出来的TD的合理性。

作者考虑到TD的一个重要性质是：位于同一个TD内部的bins的 contact frequency profiles 有更高的一致性。因此，作者计算每个TD内的任意2个bin之间Pearson相关系数（Pearson correlation coefficient)，认为PCC越高的window size越合理。

比较TopDom和已有方法

作者将TopDom与2种已有的方法（DI[2]， HiCseg [5]）进行比较（基于数据：mESC, mCortex, hESC, IMR90）

TopDom 鉴定到的TD数目更多，size更小
基于3种质量指标（intra-TD PCC, intra-TD wPCC, Intra-inter Diff），TAD表现均更优（IMR90除外）

![Table 1 | 比较3种方法鉴定的TD的大小与数目] (https://img.haomeiwen.com/i21550152/a06e48a22bbcbd22.png?imageMogr2/auto-orient/strip%7CimageView2/2/w/1240)

Figure 5 | 比较3种方法鉴定的TD的质量

参考文献

[1] Sexton,T., Yaffe,E., Kenigsberg,E., Bantignies,F., Leblanc,B., Hoichman,M., Parrinello,H., Tanay,A. and Cavalli,G. (2012) Three-dimensional folding and functional organization principles of the Drosophila genome. Cell, 148, 458–472
[2] Dixon,J.R., Selvaraj,S., Yue,F., Kim,A., Li,Y., Shen,Y., Hu,M., Liu,J.S. and Ren,B. (2012) Topological domains in mammalian genomes identified by analysis of chromatin interactions. Nature, 485, 376–380.
[3] Hou,C., Li,L., Qin,Z.S. and Corces,V.G. (2012) Gene density, transcription, and insulators contribute to the partition of the Drosophila genome into physical domains. Mol. Cell, 48, 471–484.
[4] Filippova,D., Patro,R., Duggal,G. and Kingsford,C. (2014) Identification of alternative topological domains in chromatin. Algorithms Mol. Biol., 9, 14.
[5] Levy-Leduc,C., Delattre,M., Mary-Huard,T. and Robin,S. (2014) Two-dimensional segmentation for analyzing Hi-C data. Bioinformatics, 30, I386–I392.
[6] Rao,S.S., Huntley,M.H., Durand,N.C., Stamenova,E.K., Bochkov,I.D., Robinson,J.T., Sanborn,A.L., Machol,I., Omer,A.D. and Lander,E.S. (2014) A 3D map of the human genome at kilobase resolution reveals principles of chromatin looping. Cell, 159,
1665–1680.
[7] Ray,B.K. and Ray,K.S. (1993) Determination of optimal polygon from digital curve using L 1 norm. Pattern Recognit., 26, 505–509.
[8] Ray,B.K. and Ray,K.S. (1994) A non-parametric sequential method for polygonal approximation of digital curves. Pattern Recognit. Lett., 15, 161–167.