Pore-C 多路染色体构象捕获¶

# === Pore-C 数据分析流程 ===

# 1. Basecalling（使用 Dorado）
dorado basecaller sup pod5_dir/ > pore_c.bam  # 超高精度模型识别碱基

# 2. 使用 pore-c-py 进行虚拟消化和比对
# 安装 pore-c-py
pip install pore-c-py                  # ONT官方Pore-C分析工具

# 3. 创建酶切片段参考（virtual digest）
pore_c_py refgenome virtual-digest \   # 虚拟消化参考基因组
  --enzyme NlaIII \                    # 使用的限制性内切酶
  ref.fa \                             # 参考基因组
  ref_NlaIII.bed                       # 输出酶切片段BED

# 4. 比对并提取 contact
pore_c_py reads map \                  # 比对读段
  ref.fa \                             # 参考基因组
  pore_c.bam \                         # 输入读段
  mapped_output/                       # 输出目录

# 5. 提取 pairwise contacts（兼容Hi-C工具）
pore_c_py alignments to-pairs \        # 转换为pairs格式
  mapped_output/aligned.bam \          # 比对结果
  contacts.pairs                       # 输出pairs文件

# 6. 转为 cooler 格式（Hi-C标准分析）
cooler cload pairs \                   # 从pairs文件构建矩阵
  -c1 2 -p1 3 -c2 4 -p2 5 \          # 列定义
  chrom.sizes:10000 \                  # 10kb分辨率
  contacts.pairs \                     # 输入pairs
  pore_c.cool                          # 输出cool文件

# 7. 矩阵归一化
cooler balance pore_c.cool             # ICE归一化

# 8. 生成多分辨率文件
cooler zoomify pore_c.cool             # 自动生成多分辨率mcool

# 分析 Pore-C 多体互作（multi-way contacts）
import pysam                           # BAM文件处理
from collections import Counter        # 计数器

bam = pysam.AlignmentFile("mapped_output/aligned.bam", "rb")

# 按读段ID分组，提取每条读段涉及的染色体
read_contacts = {}                     # 存储每条读段的contact信息
for read in bam:                       # 遍历比对结果
    if read.mapping_quality >= 10:     # 过滤低质量比对
        qname = read.query_name        # 读段名称
        chrom = read.reference_name    # 比对到的染色体
        pos = read.reference_start     # 比对起始位置
        if qname not in read_contacts: # 初始化
            read_contacts[qname] = []
        read_contacts[qname].append((chrom, pos))  # 记录位置

# 统计多体互作
multiway = {k: v for k, v in read_contacts.items() if len(v) >= 3}
print(f"总读段数: {len(read_contacts)}")           # 所有读段
print(f"多体互作读段 (>=3 contacts): {len(multiway)}")  # 多体互作

# 统计互作涉及的染色体组合
chrom_combos = Counter()               # 染色体组合计数
for contacts in multiway.values():     # 遍历多体互作
    chroms = tuple(sorted(set(c[0] for c in contacts)))  # 去重排序
    chrom_combos[chroms] += 1          # 计数

# 打印最常见的多体互作模式
for combo, count in chrom_combos.most_common(10):  # 前10
    print(f"  {' - '.join(combo)}: {count}")        # 打印