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

Hi-C = 给基因组的3D结构"拍照"

原理（像做三明治）：
1. 交联（Cross-link）: 用甲醛把空间上靠近的DNA片段"粘在一起"
2. 酶切（Digest）: 用限制性内切酶把DNA切成片段
3. 连接（Ligate）: 让被"粘在一起"的片段连接起来
4. 反交联+测序: 解开粘连，测序这些连接产物

结果：
  如果基因A和基因B在空间上靠近 → 它们的DNA会被连在一起 → 测序时会一起出现
  接触频率 ∝ 空间距离的反比

Hi-C数据的层级结构：
  ┌─────────────────────────────────┐
  │ Compartment（Mb级）             │  A=活跃，B=抑制
  │   ┌───────────────────────┐     │
  │   │ TAD（100kb-1Mb级）    │     │  拓扑关联域
  │   │   ┌─────────────┐    │     │
  │   │   │ Loop（kb级）│    │     │  染色质环
  │   │   └─────────────┘    │     │
  │   └───────────────────────┘     │
  └─────────────────────────────────┘

# ========== 安装Juicer ==========
# 克隆Juicer仓库
git clone https://github.com/aidenlab/juicer.git  # 下载Juicer
cd juicer

# 安装Juicer Tools
wget https://github.com/aidenlab/Juicer/releases/download/v2.20.00/juicer_tools_2.20.00.jar  # 下载jar包
mv juicer_tools_2.20.00.jar scripts/common/juicer_tools.jar  # 放到脚本目录

# 安装依赖
conda install -c bioconda bwa samtools  # 安装BWA和SAMtools

# 准备基因组
mkdir -p references/hg38
# 放入参考基因组和BWA索引
bwa index references/hg38/hg38.fa  # 建BWA索引

# 生成限制性酶切位点文件
python scripts/common/generate_site_positions.py \
    MboI \                      # 限制性酶名称（MboI或HindIII）
    hg38 \                      # 基因组名称
    references/hg38/hg38.fa     # 基因组文件
# 生成 hg38_MboI.txt

# ========== 准备目录结构 ==========
mkdir -p project/fastq  # FASTQ文件目录
# 把Hi-C的FASTQ文件放入fastq目录（R1和R2）

# ========== 运行Juicer ==========
bash scripts/juicer.sh \
    -d project/ \               # 项目目录
    -g hg38 \                   # 基因组名称
    -s MboI \                   # 限制性酶
    -z references/hg38/hg38.fa \  # 基因组FASTA
    -p references/hg38/hg38.chrom.sizes \  # 染色体大小
    -y references/hg38/hg38_MboI.txt \  # 酶切位点文件
    -t 16                       # 16线程

# Juicer流程包括：
# 1. BWA比对
# 2. 过滤无效连接产物
# 3. 去重
# 4. 生成统计报告
# 5. 生成.hic文件

# 输出文件：
# project/aligned/inter.hic        ← 主要结果文件
# project/aligned/inter_30.hic     ← MAPQ≥30过滤版
# project/aligned/inter.txt        ← 文本格式接触对
# project/aligned/stats_dups.txt   ← 统计信息

# ========== 用Juicer Tools提取矩阵 ==========
# 提取特定分辨率的接触矩阵
java -jar juicer_tools.jar dump observed \
    KR \                        # 标准化方法：KR=Knight-Ruiz（推荐）
    project/aligned/inter_30.hic \  # .hic文件
    chr1 chr1 \                 # 区域：chr1 vs chr1（cis接触）
    BP \                        # 单位：BP=碱基对
    25000 \                     # 分辨率：25kb
    chr1_contact_matrix.txt     # 输出文件

# 提取跨染色体接触
java -jar juicer_tools.jar dump observed \
    KR \
    project/aligned/inter_30.hic \
    chr1 chr2 \                 # chr1 vs chr2（trans接触）
    BP 100000 \                 # 100kb分辨率
    chr1_chr2_contact.txt

# ========== 调用TAD ==========
java -jar juicer_tools.jar arrowhead \
    -k KR \                     # 标准化方法
    -r 10000 \                  # 分辨率10kb
    project/aligned/inter_30.hic \  # .hic文件
    tad_output/                 # 输出目录

# ========== 调用Loop ==========
java -jar juicer_tools.jar hiccups \
    -k KR \                     # 标准化方法
    -r 5000,10000,25000 \       # 多分辨率
    project/aligned/inter_30.hic \
    loop_output/                # 输出目录

# ========== A/B Compartment分析 ==========
java -jar juicer_tools.jar eigenvector \
    KR \                        # 标准化
    project/aligned/inter_30.hic \
    chr1 \                      # 染色体
    BP 100000 \                 # 100kb分辨率
    eigenvector_chr1.txt        # 输出特征向量

#!/usr/bin/env python3
"""Python分析Hi-C接触矩阵"""

import numpy as np  # 数值计算
import matplotlib.pyplot as plt  # 画图
from matplotlib.colors import LogNorm  # 对数颜色标准化

# ========== 读取接触矩阵 ==========
def read_contact_matrix(filename, resolution, chrom_size):
    """读取Juicer dump输出的接触矩阵"""
    n_bins = chrom_size // resolution + 1  # 计算bin数量
    matrix = np.zeros((n_bins, n_bins))  # 初始化空矩阵

    with open(filename, 'r') as f:
        for line in f:
            parts = line.strip().split('\t')
            i = int(parts[0]) // resolution  # 行索引
            j = int(parts[1]) // resolution  # 列索引
            value = float(parts[2])  # 接触频率
            if i < n_bins and j < n_bins:
                matrix[i][j] = value
                matrix[j][i] = value  # 对称矩阵

    return matrix

# ========== 可视化接触矩阵 ==========
def plot_hic_matrix(matrix, title="Hi-C Contact Map", output="hic_map.png",
                    region=None, vmax=None):
    """画Hi-C接触热图"""
    if region:
        matrix = matrix[region[0]:region[1], region[0]:region[1]]  # 截取区域

    plt.figure(figsize=(10, 10))
    plt.imshow(
        matrix,
        cmap="Reds",  # 红色色板
        norm=LogNorm(vmin=1, vmax=vmax or np.percentile(matrix[matrix>0], 98)),  # 对数标准化
        interpolation="none"
    )
    plt.colorbar(label="Contact Frequency", shrink=0.8)
    plt.title(title, fontsize=16)
    plt.xlabel("Genomic Position (bins)", fontsize=12)
    plt.ylabel("Genomic Position (bins)", fontsize=12)
    plt.tight_layout()
    plt.savefig(output, dpi=300)
    plt.close()
    print(f"Hi-C热图已保存到 {output}")

# 使用示例
# matrix = read_contact_matrix("chr1_contact_matrix.txt", 25000, 248956422)
# plot_hic_matrix(matrix, "Chr1 Hi-C Contact Map (25kb)", region=(0, 400))

#!/usr/bin/env python3
"""使用cooler分析Hi-C数据（推荐的现代方法）"""

# pip install cooler cooltools
import cooler  # Hi-C数据格式库
import cooltools  # Hi-C分析工具包

# ========== 加载.cool文件 ==========
# .hic转.cool
# hic2cool convert input.hic output.mcool -r 0  # 转换格式

clr = cooler.Cooler("output.mcool::resolutions/10000")  # 加载10kb分辨率
print(f"分辨率: {clr.binsize}")
print(f"染色体: {clr.chromnames}")

# 获取接触矩阵
matrix = clr.matrix(balance=True).fetch("chr1:0-10000000")  # 获取chr1前10Mb的矩阵
print(f"矩阵大小: {matrix.shape}")

# ========== A/B Compartment分析 ==========
# 用cooltools计算
# eigvals, eigvecs = cooltools.eigs_cis(clr, n_eigs=3)  # 计算特征向量

# ========== TAD边界检测 ==========
# insulation_table = cooltools.insulation(clr, [10000, 25000, 50000])  # 绝缘分数

========================================
3D基因组Hi-C分析Juicer 速查表
========================================

【Hi-C层级结构】
Compartment(Mb)      → A=活跃区，B=抑制区（PCA分析）
TAD(100kb-1Mb)       → 拓扑关联域（Arrowhead算法）
Loop(kb)             → 染色质环（HiCCUPS算法）

【Juicer流程】
1. 比对               → BWA MEM
2. 过滤               → 去除无效配对
3. 去重               → 去除PCR重复
4. 生成.hic           → 多分辨率接触矩阵

【Juicer Tools关键命令】
dump                 → 提取接触矩阵
arrowhead            → TAD调用
hiccups              → Loop调用
eigenvector          → A/B Compartment

【标准化方法】
KR                   → Knight-Ruiz（推荐）
VC                   → Vanilla Coverage
SCALE                → 平衡标准化

【常用分辨率】
1Mb                  → Compartment分析
100kb                → TAD概览
25kb                 → TAD精细
10kb                 → Loop分析
5kb                  → Loop精细

【质量指标】
有效配对比例          → >40%
Cis/Trans比值        → >1（cis多于trans）
测序深度             → >200M reads（低分辨率）

【面试考点】
Q: TAD是什么？
A: 拓扑关联域，基因组中自我互作强于区间互作的区域

Q: A/B Compartment怎么区分？
A: PCA第一特征向量正值=A（活跃），负值=B（抑制）

Q: Hi-C和4C/5C的区别？
A: Hi-C检测全基因组所有互作，4C只看一个位点
========================================