Insulation score（IS），俗称绝缘分数，用于计算识别三维基因组中的拓扑关联结构域TAD。

首次提出是在：

1，概念

为染色体上的基因组区间分配‘绝缘评分’的方法。该评分用于衡量跨越每个区间的所有相互作用的总和。绝缘曲线的最低点代表高绝缘区域，我们将其定义为TAD边界。

关键点来了，很多人一开始看这个定义的时候很容易忽略一个词——across，

指的是跨过这个bin的互作的总和，什么意思呢？

我可以滑动窗口（在hic矩阵的对角线上），我窗口有大有小不是关键，包含某个bin，但是我计算的却是across/或者说是跳过这个bin的互作总和。

是不是很难理解？

我举个例子就很好懂了，

如果这个bin是一个边界（boundary），那么互作能不能跨过这个边界传递给远方，理论上是不能的，所以across跨越边界的互作总和是很小的；

反之，如果某个bin不是边界，就是正常的在TAD内部的一个区域，那么跨越这个bin的互作之和，和正常的TAD内的互作之和有什么区别吗，其实差别不会很大。

所以这个值，insulation score，顾名思义，我们把它当作是一个类似于CNN中局部local receptor field的话，那么它关注的其实就是局部的互作信息，而且这个局部感受野还挖掉了中心的eye。

这个信息只有在遇到边界的时候才能够体现出互作的显著差异，

即边界所分隔开的区域，互作之和，肯定显著低于没有被隔开的、同样大小的区域的互作之和。

绝缘分数本质——量化基因组位点对交互的阻碍强度。

2，识别算法

首先我先举一个常见的理解错误，也不能说是错误，应该是误区（实际上insulation score在2015年这篇nature文献发表之后有很多改版和变种，所以也不乏有些文献可能就是我下面的这种想法）；

Warning：下面这个是我刚开始的想法以及伪代码实现，并不是我按照2015年这篇算法开发文献的原汁原味的理解（要看对应的伪代码的直接跳过这里，看我后面第2块伪代码）！！！！！！！！！！！！！！

  TAD边界检测算法 (Hi-C数据分析)                         
───────────────────────────────────────────────────────输入:                                                 • hic_matrix: ICE校正矩阵(10kb/50kb bins)            • window_size=500kb (绝缘分析窗口)                   输出:                                                 • boundaries: [[start,end],...] (边界区域列表)        
───────────────────────────────────────────────────────1. 初始化参数                                          bin_size ← 10kb                                    window_bins ← window_size/bin_size                 insulation_scores ← array[n]                      2. 计算绝缘分数 (∀ i ∈ [0,n-1]):                       if i∉[window_bins, n-window_bins-1]:                scores[i] ← NaN                                else:                                             ▸ 方块区域: rows/cols ∈ [i±window_bins]         ▸ score = (∑方块 - M[i,i]) / (有效bins数)        3. 归一化:                                           avg ← mean(scores[有效区域])                        scores ← log₂(scores/avg)                         
───────────────────────────────────────────────────────4. 边界检测:                                           a) 计算delta导数(100kb窗口):                        δ[i] = mean(left_diff) - mean(right_diff)       b) 候选边界:                                        ∀i where δ[i]>0 ∧ δ[i+1]<0                      c) 强度过滤:                                        strength = (left_peak - score) + (score - right_valley)keep if strength ≥ 0.1                          
───────────────────────────────────────────────────────                                                        
5. 输出边界区域:                                       ∀boundary → [boundary±3 bins] (±30kb)

这里我们要区分几个概念

bin_size：这个很简单，就是hic分析中的分辨率，简单来说就是将一条染色体分成以多少kb为单位的区域，比如说10kb为单位的区域，划分出来了上万个bin，那么bin_size就是10kb，不熟悉的可以查看我之前hic分析pipeline系列的文章。
matrix_size：其实严格来说没有这个概念，只不过是我在这里为了分析透彻才提出的，其实就是hic互作矩阵的大小。举例来说，chr1我们有1个互作矩阵图谱，假设是1Mb（假设），那么1Mb就是这个矩阵的size，我们所有要识别的TAD等分析手段，都是在这个大背景下的matrix中实现的。我们的矩阵是1Mb x 1Mb，如果考虑分成bin的话，比如说bin_size还是10kb，那么这个matrix就是100bin x 100bin。
window_size：这个没什么好说的了，就是窗口的大小，但是我上面这里写的window其实不是我们常规意义理解中的滑动窗口，按照我的理解，这里改成flank_region其实更合理，就是侧翼区域，其实就是将每个bin为中心，上下游各自延伸出来window_size，比如说window_size是5个bin的话，那么就是以当前bin为中心，在hic矩阵的对角线上上下各自延伸5个bin，那么这个square其实就是5x2+1=11个bins了。
square（diamond window）：square其实才是我们在宏观角度上看来的，真实在hic矩阵对角线上滑动的那个窗口。

（1）绘制insulation score-bin index曲线

首先我们肯定有一个全长范围的matrix，比如说chr1，我们假设是1Mb x 1Mb的窗口大小，这个是我们要研究的全景范围。

然后比如说我们的hic矩阵的分辨率也就是bin_size是10kb大小，

那么这个全景窗口也就是matrix，我们实际上就有100bin x 100bin划分的网格区域——》这个时候只是有了分bin的大小，也就是这个网格划分出了大小；

然后我们在有了这个窗口也就是田字格之后，我们再引入观察的窗口window（这里观察的窗口实际上是为了我们上面提到的square所服务的），所以我们还是直接提square为主，这里的square，比如说我们使用1个长度为500kb的square的矩形窗口去滑动，那么实际上这个square就是一个50bin x 50bin大小的窗口。

也就是说1个50bin x 50bin大小的正方形，在1个100bin x 100bin大小的背景窗口下沿着对角线滑动，那么可以获得多少个滑动子序列呢？ 100-50+1个，其实对角线上就是1D序列，我们还是称这些子序列为square。

总之每一个square，因为里面都是50bin x 50bin，就算是1bin x 1bin，我们也可以计算出或者是统计出contact数来，所以我们可以获取该square中的contact总和，每一个square都可以做到。

然后对于每一个square，因为是有50bin x 50bin，理论上是可以计算出对每一个1bin x 1bin的均值（也就是将这个square的contact总和/单位bin分箱的总个数），也就是10kb的1bin，问题是这个均值，到底给square里的谁呢？毕竟这个square里有50个bin的对角线，难道这50个对角线上的bin都分配到这个square里平均的contact吗？

答案是否定的，这个平均下来的contact实际上就是score，score是与bin一一对应的，所以这个score实际上是分配给对应的bin？

那这个bin到底是谁呢？

从直接上讲，我们肯定是下意识地觉得这个bin就是中心化的bin，这是符合直觉的。

并且，method里提到了在matrix首尾超过500kb的bin，实际上就不分配score，这不是正好意味着我们的score分配原则是依据square定的吗！

这里的square，如果沿着hic矩阵的对角线滑动，正好在左上角以及右下角划出了square size的范围之外，实际上就相当于是没有这个square，那么square就没有对应的bin了，那自然就也没有对应的score了。

所以，从理论上讲，其实这里的score是和square对应的，也就是和bin对应的，

所以逻辑上bin-square-score一一对应，可以绘制出score-bin index曲线，那么如何做到这一点，对于一个就有很多个bin在对角线上的square，符合直觉的方式也就只有分配给中心的bin了。

边界处理：对于距离矩阵起点或终点500 kb以内的分箱，不分配绝缘分数，因为500 kb × 500 kb的窗口会超出矩阵范围。

这个其实就很好理解了，大背景窗口window，在对角线上左上角以及右下角，分别划出来1个square作为边界，我们的square滑动的范围界限就是这两个边界square之间；

所以我们的问题，其实是要解决每个对角线上的bin如何计算出insulation score的问题————》

一旦计算出来了，每个bin都有1个对应的insulation score，那么我们就可以做归一化，就是对角线上所有的bin的score求和，然后分子是某个bin的score，分母是所有bin求和的score，然后比值可以log2对数化。

这个时候其实对于每一个bin，如果x轴是bin index，y轴是score的话，其实就已经可以画出来1个score曲线了；

然后这个score曲线的谷底，也就是min，可以看作是hic互作急剧下降的位点，这些bin其实就被解释为TAD边界。

所以这里的理解关键就是，如何将滑动的square与bin一一对应上？

还是使用文献中的方法，其实这里的描述已经很清晰了：

across each bin (up to a specified distance up- and downstream of the bin)

我们的bin（被across的），是如何被组织成需要聚合的square的contact，其实就是在这个bin上下游划分或者说是各自延伸一定distance，这里的distance其实就是我上面提到的window_size的概念，

然后我们每一个bin，有了上下游的window_size，也就是bin±window_size，才组织成了square这个概念。

而且这里其实也说明了，我们的window_size也就是侧翼flank的区域大小，其实是自定义的（相当于是1个超参数）。

当然，实际上这个算法有很多变种版本，主要是关于其中的square的定义部分，

不同的变种版本，具体设置的[Xstart，Xend]x[Ystart，Yend]实际上是不一样的；

在我上面理解的版本中，由[Xstart，Xend]x[Ystart，Yend]构成的square（也就是diamond window），实际上采用的是是在中心化的bin的左右两侧都flank一个window出来；

有的版本中，比如说下面的：

这也和我们前面的伪代码中减去中心的这个bin的contact的某种思想是对应的。

好了，以上就是我提到的第一个误区，也就是我一开始的想法！！！！！！！！

那么，真实原文实现的细节是怎么样的呢？

参考具体原始实现细节：

https://github.com/dekkerlab/crane-nature-2015

原始脚本是使用perl写的，详细细节可以参考，

我们需要重点关注的细节是这个square方块如何选取的，原始使用的脚本中对于这个square的坐标的设计：

如果我们还是使用500kb的square_size来举例的话，那么这一个square选取的其实不是一个对称矩阵：

其实结合这个图形，就更容易理解了，我们这会正经看一下这个图：

所谓的across其实和“do not include the diagonal”是一致的，

我们仔细看这个热图，发现了什么没有（一般越是CNS主刊的图越是不可能随意画的，所以这就算是一个简单的算法示意图，还是有细节可以考究的）。

我们注意看这个热图中滑动的square的起点以及终点，

我们可以发现，左上角的第一个square，实际上只和标记mark为3的那个bin有接触

从这个图里，我们就很容易发现并理解，为什么bin1和2就没有对应的square，因为这个square再接着按照图中的轨迹向左上方移动就会超出边界线。

再仔细看一下，很明显这个square的window_size就有2个bin大小，当然这个不是问题的关键，可能有的人会有疑问，如果我把这个square放在对角线上，不也同样是能够获取mark 2的bin的信息的吗？

那按理来说我不应该是从index为2的bin开始就有1个score吗？

但实际上这个示意图中只有从index 3开始的score信息！

所以，这其实就是理解算法中across的关键，我们恰恰可以从这一个问题出发：为什么这个square不放在对角线上，看形式上，它是放在对角线的上方的，只能说是沿着along对角线，不能说是on对角线。

我们再看一下这种square放置方法的效果是什么：

我们可以仔细看这个和index为3的bin有接触的square，

我们可以把对角线上的坐标映射到x轴以及y轴上

我们可以发现什么，这个square包含的是不是mark1、2的bin和mark4、5的bin之间的互作contact之和！

少了什么？

就是bin3内部的contact！

而1、2和4、5不就是bin3的上下游吗，所以就是来统计评估在bin3上下游的互作情况，而且是指定上游和下游互作强度如何。

我们很容易就可以推论，如果bin3是一个边界的话，那么bin3上游的bin1、2和下游的bin4、5不就正好被绝缘了吗，相当于是正好对应上游的TAD以及下游的TAD，

所以这个时候两个TAD之间的互作绝缘的话，那么bin1、2和bin4、5统计出来的contact应该会很低，所以我们说bin3也就是绝缘边界的话insulation score会很低。

反之，如果bin3不是TAD边界，而是在某一个TAD的内部，那么在一个window_size范围内，其实邻居bin还是在TAD内的，那么依据TAD的定义，内部互作依然是频繁的，所以这个时候bin1、2和bin4、5统计出来的contact应该会很高。

那么两个比较之下，很显然，我们就能够依据insulation score的高低来区分这个bin是在边界上还是在TAD domain内部。

上面都是定性的表述，下面我们来定量描述一下：

我们可以看到这里的bin index为3，也就是我们假设i=3，对于目前要处理的bin；

那么依据我们前面对于bin-square-score一一对应的理论，这里的square应该如何组织呢？

square无非是由4个坐标决定，也就是[Xstart，Xend]x[Ystart，Yend]。

很明显，这里i=3的bin对应的square是[4，5]x[1，2]。

我们可以用一种更加容易理解的方式来重写这个square的坐标：
[3+1，3+2]x[3-2，3-1]，

这里其实就可以联系我们的square_size了，

对标index=i的bin，其实用于收集其对应across的square坐标是[i+1，i+square_size]x[i-square_size，i-1]。

然后因为是对称矩阵，所以square放在对角线右上角还是左下角其实是没有区别的，都可以收集到across这个bin的互作contact。

简单来理解就是bin i右上方邻接的square_size x square_size的正方形。

从这个角度来讲，其实就很好理解我下面这里所写出来的insulation score的伪代码了：

  TAD边界检测算法 (Hi-C数据分析)                         
───────────────────────────────────────────────────────输入:                                                 • hic_matrix: ICE校正矩阵(10kb bins)            • window_size=500kb (绝缘分析窗口)                   输出:                                                 • boundaries: [[start,end],...] (边界区域列表)        
───────────────────────────────────────────────────────1. 初始化参数                                          bin_size ← 10kb                                    window_bins ← window_size/bin_size                 insulation_scores ← array[n]                      2. 计算绝缘分数 (∀ i ∈ [0,n-1]):                       if i∉[window_bins, n-window_bins-1]:                scores[i] ← NaN                                else:                                             ▸ 方块区域square: rows ∈ [i-window_bins, i-1]cols ∈ [i+1, i+window_bins]      ▸ score = (∑方块) / (有效bins数)        3. 归一化:                                           avg ← mean(scores[有效区域])                        scores ← log₂(scores/avg)                         
───────────────────────────────────────────────────────4. 边界检测:                                           a) 计算delta差分(上下游各100kb窗口):                        δ[i] = mean(left_diff) - mean(right_diff)       b) 候选边界:                                        ∀i where δ[i]>0 ∧ δ[i+1]<0                      c) 强度过滤:                                        strength = (left_peak - right_valley)keep if strength ≥ 0.1                          
───────────────────────────────────────────────────────                                                        
5. 输出边界区域:                                       ∀boundary → [boundary±3 bins] (±30kb)

注意：实际坐标或者说bin的index会因为0、1-indexed的差异，或者在实际hic矩阵中bin是作为一个区间还是一个起点——导致上面的具体实现细节有点差异，但是大致算法上的逻辑，基本上应该是和我上面这个理解兼容的。

（2）有了这个曲线之后我们再讨论如何寻找min，也就是谷底

这个时候其实就很简单了，如果我们按照（1）中的方法，已经绘制出了一条insulation score-bin index的曲线，我们其实只需要找这个极小值点就行了（按照我们对TAD边界在insulation score中的表现和定义）。

求极值点其实很简单，学过微积分的都知道，如果是连续函数，就求导，一阶导=0的点，也就是极值点。当然这里我们要注意一下，极值点对应的不一定是valley谷底，也可能是peak山峰，也就是极小值还是极大值。

如果是离散的话，就差分，微分方程对应差分方程，其实就是数值分析那一块的内容。

所以delta法其实就是差分法，同样的，delta=0的地方可能是极大值也可能是极小值，比如说上面图a红线跨0点的地方，对应的黑线部分，可能是insulation score max或min的，其中min的坐标是我们需要的。

然后，注意到我们在计算insulation score的时候，是针对每一个bin进行计算的，所以是bin和score一一对应的；

现在，我们在这里同样计算delta，也是围绕每一个参考bin进行的，所以也是bin和delta一一对应的。

然后这个delta计算其实就很简单了，就是（左侧互作-参考bin）的均值-（右侧互作-参考bin）的均值，至于左侧和右侧到底size多大，还是要定义1个window，相当于还是有一个超参数

当然，还是那个边界的问题，对于delta，我们还是要考虑在参考bin上游的窗口能否取到，理论上delta因为要计算前向值，那么前面没有score的bin，在计算delta向量的时候是否会有边界问题？

当然注意一下，这里的差分delta，其实和我们正常理解的差分不一样，我们正常理解的差分，是用当前bin的score-之前bin的score，但是这里指的是上游bin的score相对参考点的score，所以这里理解的delta的值和我们常规理解的差分其实就差了个方向。

那么我们要从-delta的角度来类比微分，然后理解极小值点和极大值点是什么，

-delta[i]<0,-delta[i+1]>0，valley，对应就是谷底，也就是边界

当然理想情况下是delta=0，

我们还是对照前面官方原始的perl脚本

for(my $y=0;$y<$numHeaders;$y++) {# 获取当前bin的头部信息my $yHead=$inc2header->{ y }->{$y};my $insulation="NA";$insulation=$matrixInsulation->{$yHead} if(exists($matrixInsulation->{$yHead}));# 计算左右边界my $leftBound=($y-$halfInsulationDeltaBinSize);$leftBound=0 if($leftBound < 0);my $rightBound=($y+$halfInsulationDeltaBinSize);$rightBound=$numHeaders-1 if($rightBound >= $numHeaders);# 计算左侧delta值my @leftDelta=();for(my $ly=$leftBound;$ly<$y;$ly++) {my $leftHead=$inc2header->{ y }->{$ly};my $leftInsulation="NA";$leftInsulation=$matrixInsulation->{$leftHead} if(exists($matrixInsulation->{$leftHead}));my $delta="NA";$delta=($leftInsulation-$insulation) if(($insulation ne "NA") and ($leftInsulation ne "NA"));push(@leftDelta,$delta) if($delta ne "NA");}my $leftDeltaAggregrate="NA";if(@leftDelta > 0) {my $leftDeltaStats=listStats(\@leftDelta);$leftDeltaAggregrate=$leftDeltaStats->{ mean };}# 计算右侧delta值my @rightDelta=();for(my $ry=$y+1;$ry<=$rightBound;$ry++) {my $rightHead=$inc2header->{ y }->{$ry};my $rightInsulation="NA";$rightInsulation=$matrixInsulation->{$rightHead} if(exists($matrixInsulation->{$rightHead}));my $delta="NA";$delta=($rightInsulation-$insulation) if(($insulation ne "NA") and ($rightInsulation ne "NA"));push(@rightDelta,$delta) if($delta ne "NA");}my $rightDeltaAggregrate="NA";if(@rightDelta > 0) {my $rightDeltaStats=listStats(\@rightDelta);$rightDeltaAggregrate=$rightDeltaStats->{ mean };}# 计算deltaDelta值my $deltaDelta="NA";$deltaDelta=($leftDeltaAggregrate-$rightDeltaAggregrate) if(($leftDeltaAggregrate ne "NA") and ($rightDeltaAggregrate ne "NA"));$matrixDelta->{$yHead}=$deltaDelta;
}

其实代码逻辑很简单，实现的是同样的原理：

遍历每个bin：
- 使用for(my $y = 0;$ y< $n u m He a d ers;$ y++)循环遍历矩阵中的每个bin。
获取当前bin的绝缘指数：
- 从 $ma t r i x I n s u l a t i o n 哈希中获取当前 bin 的绝缘指数（$ insulation）。
- 如果当前bin的绝缘指数不存在，则设置为"NA"。
计算左右边界：
- 计算当前bin的左边界（ $l e f tB o u n d ）和右边界（$ rightBound）。
- 确保边界不会超出矩阵的范围。
计算左侧delta值：
- 遍历当前bin左侧的bin，计算每个左侧bin的绝缘指数与当前bin的绝缘指数的差值（$delta）。
- 将有效的delta值存储在@leftDelta数组中。
- 如果@leftDelta数组中有数据，计算其平均值（$leftDeltaAggregrate）。
计算右侧delta值：
- 遍历当前bin右侧的bin，计算每个右侧bin的绝缘指数与当前bin的绝缘指数的差值（$delta）。
- 将有效的delta值存储在@rightDelta数组中。
- 如果@rightDelta数组中有数据，计算其平均值（$rightDeltaAggregrate）。
计算deltaDelta值：
- 计算左侧delta平均值与右侧delta平均值的差值（$deltaDelta）。
- 将计算结果存储在$matrixDelta哈希中。

当然使用极小值法找出来边界之后，我们还有一个叫做boundary strength的指标来过滤，

当前bin左侧的局部极大值-右侧的局部极小值

然后计算这个边界的话，其实也很简单，因为我们知道找到的每一个候选的boundary，实际上都是insulation的局部极小值点，按照delta的逻辑和方向的定义，又是delta的右半侧下降曲线交叉零点的地方，

既然是右半侧下降曲线，那么这个下降曲线，向左的边界是局部极大值（delta的局部极大值，不是insulation的），向右的边界是局部极小值。

所以原脚本中的逻辑就是：
通过向左和向右搜索，找到当前点两侧的“delta”值变化趋势的边界，并计算这两个边界之间的差值作为边界强度。如果计算出的边界强度小于噪声阈值，则认为这个边界是无效的，从而跳过它。

对照原始脚本，就是：

$leftSearchIndex = $yIndex - 1;
while( ($leftSearchIndex > 0) and ($matrixDelta->{$inc2header->{ y }->{$leftSearchIndex}} ne "NA") and ($matrixDelta->{$inc2header->{ y }->{$leftSearchIndex-1}} ne "NA") and ($matrixDelta->{$inc2header->{ y }->{$leftSearchIndex-1}} >= $matrixDelta->{$inc2header->{ y }->{$leftSearchIndex}}) 
) {$leftSearchIndex--;  # 向左移动直到Delta不再递增
}
my $leftDeltaBound = $matrixDelta->{$inc2header->{ y }->{$leftSearchIndex}};

向左搜索局部极大值（Peak）：

搜索逻辑：
从候选边界向左移动，找到Delta值的局部最大值（即曲线上升的终点）
终止条件：Delta开始下降（前一个Delta < 当前Delta）或遇到数据缺失
生物学意义：左侧区域的最高绝缘分数变化率

向右搜索局部极小值（Trough）

$rightSearchIndex = $yIndex + 1;
while(($rightSearchIndex < ($numHeaders-2)) and ($matrixDelta->{$inc2header->{ y }->{$rightSearchIndex}} ne "NA") and ($matrixDelta->{$inc2header->{ y }->{$rightSearchIndex+1}} ne "NA") and ($matrixDelta->{$inc2header->{ y }->{$rightSearchIndex+1}} <= $matrixDelta->{$inc2header->{ y }->{$rightSearchIndex}})
) {$rightSearchIndex++;  # 向右移动直到Delta不再递减
}
my $rightDeltaBound = $matrixDelta->{$inc2header->{ y }->{$rightSearchIndex}};

搜索逻辑：
从候选边界向右移动，找到Delta值的局部最小值（即曲线下降的终点）
终止条件：Delta开始上升（后一个Delta > 当前Delta）或遇到数据缺失
生物学意义：右侧区域的最低绝缘分数变化率

然后就是这个边界强度

my $valleyDeltaStrength = "NA";
$valleyDeltaStrength = ($leftDeltaBound - $rightDeltaBound) if(($leftDeltaBound ne "NA") and ($rightDeltaBound ne "NA"));

计算公式：
Strength=LeftPeakDelta−RightTroughDeltaStrength=LeftPeakDelta−RightTroughDelta

综上，再次重申我理解的伪代码：

  TAD边界检测算法 (Hi-C数据分析)                         
───────────────────────────────────────────────────────输入:                                                 • hic_matrix: ICE校正矩阵(10kb bins)            • window_size=500kb (绝缘分析窗口)                   输出:                                                 • boundaries: [[start,end],...] (边界区域列表)        
───────────────────────────────────────────────────────1. 初始化参数                                          bin_size ← 10kb                                    window_bins ← window_size/bin_size                 insulation_scores ← array[n]                      2. 计算绝缘分数 (∀ i ∈ [0,n-1]):                       if i∉[window_bins, n-window_bins-1]:                scores[i] ← NaN                                else:                                             ▸ 方块区域square: rows ∈ [i-window_bins, i-1]cols ∈ [i+1, i+window_bins]      ▸ score = (∑方块) / (有效bins数)        3. 归一化:                                           avg ← mean(scores[有效区域])                        scores ← log₂(scores/avg)                         
───────────────────────────────────────────────────────4. 边界检测:                                           a) 计算delta差分(上下游各100kb窗口):                        δ[i] = mean(left_diff) - mean(right_diff)       b) 候选边界:                                        ∀i where δ[i]>0 ∧ δ[i+1]<0                      c) 强度过滤:                                        strength = (left_peak - right_valley)keep if strength ≥ 0.1                          
───────────────────────────────────────────────────────                                                        
5. 输出边界区域:                                       ∀boundary → [boundary±3 bins] (±30kb)

此处原文获取了对应作为边界的bin之后，再用了上下游30kb，也就是3个bin来延伸对应效果最好的的边界。

如果要严谨对应源代码perl脚本一点，就是

TAD边界检测算法 (Hi-C数据分析)
───────────────────────────────────────────────────────
输入:• hic_matrix: 对称的ICE校正矩阵(等距bin)• window_size=500kb (默认绝缘窗口)• delta_span=window_size/2 (默认差分窗口)
输出:• *.boundaries (边界区域bed)
───────────────────────────────────────────────────────
1. 初始化参数bin_size ← headerSpacing (从矩阵元数据获取)window_bins ← ceil((window_size - bin_size)/bin_size) +1delta_bins ← ceil(delta_span/bin_size)if delta_bins是偶数: delta_bins +=1 (确保奇数窗口)2. 绝缘分数计算:for y in [window_bins, n-window_bins-1]:startY = y - window_binsendY = ystartX = y + 1endX = y + window_bins +1boxData = []for y2 in [startY, endY):for x2 in [startX, endX):if matrix[y2][x2] == -7337: 跳过缺失值if ignoreZero and matrix[y2][x2] == 0: 跳过零值boxData.append(matrix[y2][x2])if 有效数据点 < 总数据点/2: 跳过该binscore = 根据insulationMode计算(boxData) # mean/sum/median3. 归一化处理:valid_scores = 过滤NaN值(scores)avg = mean(valid_scores)for score in scores:if score和avg都不为0:score = log2(score/avg)  # Perl实际使用自然对数换算else:score = score - avg  # 简单中心化4. 边界检测:a) 计算delta差分:for y in [0, numHeaders):left_mean = mean( insulation[y-delta_bins...y] - insulation[y] )right_mean = mean( insulation[y...y+delta_bins] - insulation[y] )delta[y] = left_mean - right_meanb) 候选边界检测:- 遍历delta符号变化点(正→负)- 用matrixDeltaSquare标记(+1表示正delta，-1表示负delta)c) 强度过滤:for 每个候选边界y:left_peak = 向左搜索直到delta不再增加right_valley = 向右搜索直到delta不再减少strength = delta[left_peak] - delta[right_valley]if strength < noiseThreshold(默认0.1): 剔除5. 输出结果:- 原始边界坐标 ± (boundaryMarginOfError * bin_size)

3，Implementation实现

（1）cooltools：

https://cooltools.readthedocs.io/en/latest/notebooks/insulation_and_boundaries.html

contact frequencies between upstream and downstream loci——》其实这个说法和我们对across绝缘的理解是一致的

from cooltools import insulation

其实从API中我们已经很清楚了diamond window这个概念，其实就是我上面提到的square概念，diamond window to calculate the insulation score

（2）GENOVA：

https://rdrr.io/github/robinweide/GENOVA/man/insulation_score.html

TAD boundaries can be identified as local minima in the insulation score

因为GENOVA是用R写的，所以我们可以看一下call_TAD_insulation以及tornado_insulation这两个函数的具体实现细节。

（3）FAN-C

https://vaquerizaslab.github.io/fanc/fanc-executable/fanc-analyse-hic/domains.html

Use fanc insulation to calculate the insulation score from the command line:usage: fanc insulation [-h] [-o OUTPUT_FORMAT][-w WINDOW_SIZES [WINDOW_SIZES ...]] [-r REGION] [-i][--offset OFFSET] [-L] [-N][--normalisation-window NORMALISATION_WINDOW] [-s] [-g][--trim-mean TRIM_MEAN] [-tmp]input [output]