Hi-C 数据分析与 3D 基因组¶

# === Hi-C 数据分析标准流程（使用 HiC-Pro） ===

# 1. 创建配置文件（关键参数）
cat > config.txt << 'EOF'
BOWTIE2_IDX_PATH = /path/to/bowtie2_index   # 基因组索引路径
REFERENCE_GENOME = hg38                       # 参考基因组
GENOME_SIZE = /path/to/chrom_sizes.txt       # 染色体大小文件
LIGATION_SITE = GATCGATC                      # DpnII酶切+连接后的序列
MIN_MAPQ = 10                                 # 最小比对质量
EOF

# 2. 运行 HiC-Pro（从FASTQ到接触矩阵）
HiC-Pro -i rawdata/ \                  # 输入FASTQ目录
  -o hicpro_out \                      # 输出目录
  -c config.txt                        # 配置文件

# === 使用 cooler 生成多分辨率矩阵 ===

# 3. 将 pairs 文件转为 cooler 格式
cooler cload pairs \                   # 从pairs文件加载
  -c1 2 -p1 3 -c2 4 -p2 5 \          # 指定列号（chr1,pos1,chr2,pos2）
  chrom_sizes.txt:1000 \              # 分辨率1kb
  valid_pairs.txt \                    # 有效配对文件
  output.cool                          # 输出cool文件

# 4. 生成多分辨率 mcool 文件
cooler zoomify output.cool \           # 输入cool文件
  -r 1000,5000,10000,25000,50000,100000,500000,1000000  # 多种分辨率

# 5. 矩阵平衡（ICE 归一化）
cooler balance output.cool             # ICE归一化消除偏差

# === A/B Compartment 分析 ===
# 6. 使用 cooltools 计算 eigenvector
cooltools eigs-cis \                   # 计算顺式特征向量
  output.mcool::resolutions/100000 \   # 100kb分辨率
  --phasing-track gene_density.bw \    # 用基因密度校正方向
  -o compartments                       # 输出前缀

# === TAD 识别 ===
# 7. 使用 cooltools 计算绝缘分数
cooltools insulation \                 # 计算绝缘分数
  output.mcool::resolutions/25000 \    # 25kb分辨率
  25000 50000 100000 \                 # 窗口大小
  -o insulation.tsv                    # 输出文件

# 用 Python 可视化 Hi-C 接触矩阵
import cooler                          # Hi-C数据格式库
import matplotlib.pyplot as plt        # 绑图库
import numpy as np                     # 数值计算

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

# 提取某个区域的接触矩阵
region = "chr1:50000000-55000000"      # 感兴趣的区域
matrix = clr.matrix(balance=True).fetch(region)  # 提取归一化矩阵

# 绘制热图
plt.figure(figsize=(8, 8))            # 设置图形大小
plt.imshow(
    np.log1p(matrix),                  # 对数变换使信号更明显
    cmap="YlOrRd",                     # 黄-橙-红配色
    interpolation="none"               # 不插值
)
plt.colorbar(label="log(contact)")     # 添加颜色条
plt.title(f"Hi-C Contact Map: {region}")  # 标题
plt.savefig("hic_heatmap.png", dpi=150)   # 保存图片