TAD结构识别与可视化¶

TAD = 基因组的"社区划分"

想象一个城市：
  ├── 社区A（TAD-A）：区内居民频繁互动
  │     └── 增强子A → 调控基因A（在同一社区内）
  ├── 围墙（TAD边界）：有CTCF守门人
  └── 社区B（TAD-B）：区内居民频繁互动
        └── 增强子B → 调控基因B（在同一社区内）

TAD的生物学意义：
  1. 确保增强子只调控同一TAD内的基因
  2. 防止不同区域的调控元件"串线"
  3. TAD边界破坏 → 增强子可能错误激活邻近TAD的基因 → 疾病

TAD特征：
  大小：100kb - 2Mb（平均约800kb）
  数量：人类基因组约2000-3000个TAD
  保守性：不同细胞类型间TAD边界高度保守
  边界特征：CTCF结合 + 管家基因 + CpG岛

#!/usr/bin/env python3
"""TAD边界识别 - 绝缘分数法"""

import numpy as np  # 数值计算
import matplotlib.pyplot as plt  # 画图
from scipy.signal import argrelmin  # 找局部最小值

def calculate_insulation_score(matrix, window_size=10):
    """
    计算绝缘分数（Insulation Score）
    原理：在每个bin处用一个方形窗口计算区域内的平均互作强度
    绝缘分数低的位置 = TAD边界（两侧互作弱）
    """
    n = matrix.shape[0]  # 矩阵大小
    scores = np.zeros(n)  # 初始化分数数组

    for i in range(window_size, n - window_size):  # 遍历每个bin（跳过边缘）
        # 提取以i为中心的窗口区域
        window = matrix[i - window_size:i, i:i + window_size]  # 左上-右下的方形窗口
        scores[i] = np.nanmean(window)  # 计算窗口内平均互作

    # 归一化为log2比值
    mean_score = np.nanmean(scores[scores > 0])  # 全局平均
    log2_scores = np.log2(scores / mean_score)  # log2归一化
    log2_scores[np.isinf(log2_scores)] = 0  # 处理无穷值

    return log2_scores

def find_tad_boundaries(insulation_scores, delta_window=5, threshold=-0.1):
    """
    从绝缘分数中找TAD边界
    边界 = 绝缘分数的局部最小值点
    """
    # 找局部最小值
    minima = argrelmin(insulation_scores, order=delta_window)[0]  # 找局部极小值

    # 过滤：只保留低于阈值的最小值
    boundaries = [m for m in minima if insulation_scores[m] < threshold]

    return boundaries

# ========== 使用示例 ==========
# 假设已有Hi-C接触矩阵
# matrix = np.loadtxt("chr1_contact_matrix.txt")

# 模拟一个Hi-C矩阵（实际使用时替换为真实数据）
np.random.seed(42)
n = 200  # 200个bin
matrix = np.random.rand(n, n)
matrix = (matrix + matrix.T) / 2  # 对称化

# 模拟TAD结构：让某些区域内的互作更强
for start, end in [(0, 50), (50, 100), (100, 150), (150, 200)]:
    matrix[start:end, start:end] += 5  # TAD内互作增强

# 计算绝缘分数
scores = calculate_insulation_score(matrix, window_size=10)

# 找TAD边界
boundaries = find_tad_boundaries(scores, delta_window=5, threshold=-0.5)
print(f"检测到 {len(boundaries)} 个TAD边界")
print(f"边界位置: {boundaries}")

# ========== 可视化 ==========
fig, axes = plt.subplots(2, 1, figsize=(14, 10), gridspec_kw={'height_ratios': [3, 1]})

# 上图：Hi-C接触热图
from matplotlib.colors import LogNorm
im = axes[0].imshow(matrix, cmap="Reds", aspect="auto",
                     norm=LogNorm(vmin=0.1, vmax=np.percentile(matrix, 95)))
for b in boundaries:  # 在边界处画线
    axes[0].axvline(b, color="blue", linewidth=0.5, alpha=0.7)
    axes[0].axhline(b, color="blue", linewidth=0.5, alpha=0.7)
axes[0].set_title("Hi-C Contact Matrix with TAD Boundaries", fontsize=14)
plt.colorbar(im, ax=axes[0], shrink=0.8)

# 下图：绝缘分数
axes[1].plot(scores, color="steelblue", linewidth=1)
axes[1].axhline(-0.5, color="red", linestyle="--", alpha=0.5, label="Threshold")
for b in boundaries:
    axes[1].axvline(b, color="blue", linewidth=0.5, alpha=0.7)
axes[1].set_xlabel("Genomic Position (bins)", fontsize=12)
axes[1].set_ylabel("Insulation Score", fontsize=12)
axes[1].legend()

plt.tight_layout()
plt.savefig("tad_analysis.png", dpi=300)
plt.close()
print("TAD分析图已保存")

#!/usr/bin/env python3
"""使用cooltools进行TAD分析"""

# pip install cooler cooltools
import cooler
import cooltools
import pandas as pd

# 加载cooler文件
clr = cooler.Cooler("output.mcool::resolutions/25000")  # 25kb分辨率

# ========== 计算绝缘分数 ==========
insulation_df = cooltools.insulation(
    clr,
    window_bp_sizes=[100000, 250000, 500000]  # 多尺度窗口
)

# 查看结果
print(insulation_df.head())
# 列包含：chrom, start, end, log2_insulation_score_100000, ...
# 以及 boundary_strength_100000, is_boundary_100000 等

# ========== 提取边界 ==========
boundaries = insulation_df[insulation_df["is_boundary_250000"]].copy()
print(f"250kb窗口检测到 {len(boundaries)} 个TAD边界")

#!/usr/bin/env python3
"""TAD结构可视化 - 三角形热图"""

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Polygon
from matplotlib.collections import PatchCollection

def plot_tad_triangle(matrix, tad_starts, resolution=25000, output="tad_triangle.png"):
    """
    画TAD三角形热图（发表常用格式）
    将正方形接触矩阵旋转45度显示为三角形
    """
    n = matrix.shape[0]

    fig, ax = plt.subplots(figsize=(16, 6))

    # 旋转矩阵45度
    # 转换坐标：(i,j) → (x, y) = ((i+j)/2, (j-i))
    patches = []
    colors = []

    for i in range(n):
        for j in range(i, min(i + 100, n)):  # 限制显示距离
            x = (i + j) / 2
            y = j - i

            diamond = Polygon([
                (x, y), (x + 0.5, y + 1),
                (x + 1, y), (x + 0.5, y - 1)
            ])
            patches.append(diamond)
            colors.append(matrix[i, j])

    pc = PatchCollection(patches, cmap="Reds")
    pc.set_array(np.array(colors))
    pc.set_clim(0, np.percentile(colors, 95))
    ax.add_collection(pc)

    # 画TAD边界
    for start in tad_starts:
        ax.axvline(start, color="blue", linewidth=1, alpha=0.7)

    ax.set_xlim(0, n)
    ax.set_ylim(0, 100)
    ax.set_xlabel(f"Genomic Position ({resolution/1000:.0f}kb bins)")
    ax.set_title("TAD Structure (Triangle View)")

    plt.colorbar(pc, ax=ax, label="Contact Frequency", shrink=0.5)
    plt.tight_layout()
    plt.savefig(output, dpi=300)
    plt.close()
    print(f"TAD三角形热图已保存到 {output}")

========================================
TAD结构识别与可视化 速查表
========================================

【TAD识别方法】
绝缘分数法           → cooltools.insulation()（推荐）
Arrowhead            → Juicer Tools arrowhead
TopDom               → R包，快速TAD检测
HiCseg               → 贝叶斯分割法
DI方向性指数         → Dixon et al. 2012原始方法

【推荐分辨率】
TAD概览              → 50-100kb
TAD精细              → 25kb
Sub-TAD              → 10kb
Loop                 → 5-10kb

【绝缘分数参数】
窗口大小             → 250kb-500kb（常用）
边界阈值             → 绝缘分数局部最小值
多尺度分析           → 同时用多个窗口大小

【TAD特征】
大小                 → 100kb-2Mb（平均~800kb）
数量                 → 人类~2000-3000个
边界特征             → CTCF + Cohesin + 管家基因
保守性               → 不同细胞类型高度保守

【可视化工具】
JuiceBox             → .hic文件交互式浏览
HiGlass              → 在线Hi-C可视化
pyGenomeTracks       → 编程式多轨道可视化
cooltools             → Python分析+可视化

【面试考点】
Q: TAD边界有什么特征？
A: CTCF结合位点、管家基因、CpG岛、转录活跃

Q: TAD边界破坏会怎样？
A: 增强子可能错误激活邻近TAD的基因，导致疾病

Q: TAD识别用什么分辨率？
A: 一般25-50kb，太高分辨率需要更多测序量
========================================