157_Hi-C辅助基因组组装¶

一句话概述¶

Hi-C辅助基因组组装利用染色质三维交互频率信息(近距离片段交互更频繁)将contig/scaffold锚定(scaffold)到染色体水平，常用工具包括3D-DNA、SALSA2和YaHS，结合Juicebox人工校正实现高质量染色体级基因组。

核心知识点总览¶

知识点	说明
Hi-C原理	捕获三维空间中DNA的近距离交互信息
接触频率衰减	同一染色体上距离越近交互频率越高
3D-DNA	Aiden Lab开发的自动化Hi-C scaffolding流水线
SALSA2	基于图论的Hi-C scaffolding工具
YaHS	最新的快速Hi-C scaffolding工具
Juicebox	交互式Hi-C图可视化和手动校正工具
Contact map	基因组位置间Hi-C交互频率的热图
ALLHiC	针对多倍体的Hi-C组装优化工具

各步骤详解¶

第一步：Hi-C辅助组装的原理¶

白话解释： Hi-C就像给染色体拍"合影"——通过化学交联固定空间上靠近的DNA片段，然后测序找出谁和谁靠得近。在同一条染色体上的两段DNA交互频率最高，且距离越近交互越多。利用这个规律，可以把组装好的DNA片段(contigs)按正确的顺序和方向排列到染色体上。

技术原理： 1. 基因组初步组装得到contigs/scaffolds(来自PacBio/ONT + HiFi) 2. Hi-C测序产生配对reads，代表空间近邻的DNA片段对 3. 将Hi-C reads比对到contigs 4. 同一条染色体的contigs之间Hi-C链接多，不同染色体间少 5. 同一染色体内，按距离排列使交互频率单调递减 6. 算法目标：找到最优排列使contact map呈现"块状对角线"模式

流程概述：

原始reads → 初步组装(hifiasm/flye) → contigs
Hi-C reads → BWA比对到contigs → Hi-C链接信息
contigs + Hi-C链接 → scaffolding(YaHS/3D-DNA/SALSA2) → 染色体级scaffold
Juicebox手动校正 → 最终染色体基因组

第二步：Hi-C数据预处理¶

代码示例：

# === Hi-C reads比对到contigs ===
# 使用BWA-MEM2 (注意: Hi-C reads需要特殊处理)

# 1. 建立索引
bwa-mem2 index contigs.fasta

# 2. 分别比对两端reads (不用paired-end模式)
bwa-mem2 mem -t 32 -5SP contigs.fasta hic_R1.fq.gz | \
  samtools view -@ 8 -bhS -F 2316 -o aligned_R1.bam
bwa-mem2 mem -t 32 -5SP contigs.fasta hic_R2.fq.gz | \
  samtools view -@ 8 -bhS -F 2316 -o aligned_R2.bam

# 参数说明:
# -5: 作为split-read比对处理(Hi-C chimeric reads)
# -S: 跳过mate rescue
# -P: 跳过配对
# -F 2316: 过滤secondary/supplementary

# 3. 合并配对信息 (用于YaHS/SALSA2)
# 先排序
samtools sort -@ 8 -n aligned_R1.bam -o sorted_R1.bam
samtools sort -@ 8 -n aligned_R2.bam -o sorted_R2.bam

# 合并为paired BAM
# 使用pairtools或自定义脚本
# 对于YaHS, 可以直接用chromap:
chromap --preset hic \
  -x contigs.index \
  -r contigs.fasta \
  -1 hic_R1.fq.gz \
  -2 hic_R2.fq.gz \
  --SAM -o aligned.sam \
  -t 32

samtools view -bh aligned.sam | samtools sort -n -@ 8 -o hic_sorted.bam

# === 或使用Arima Genomics推荐流程 ===
# 使用Arima mapping pipeline
bwa mem -t 32 -5SP contigs.fasta hic_R1.fq.gz hic_R2.fq.gz | \
  samtools view -@ 8 -bhS - | \
  samtools sort -@ 8 -n -o hic_aligned_sorted.bam

第三步：YaHS — 快速Hi-C Scaffolding¶

白话解释： YaHS(Yet Another Hi-C Scaffolding tool)是目前推荐的scaffolding工具——速度快、内存少、准确性高。它的核心策略是反复迭代：先把小片段连成大块，检查交互模式是否合理，不合理就打断重连。

代码示例：

# === 安装YaHS ===
git clone https://github.com/c-zhou/yahs.git
cd yahs && make

# === 运行YaHS ===
yahs contigs.fasta hic_sorted.bam \
  -o yahs_output \
  --no-contig-ec \
  -q 10

# 参数说明:
# contigs.fasta: 初步组装结果
# hic_sorted.bam: name-sorted的Hi-C比对BAM
# -o: 输出前缀
# --no-contig-ec: 不做contig错误校正(如果contigs质量高)
# -q 10: 最低比对质量

# 输出:
# yahs_output_scaffolds_final.fa  : 最终scaffold序列
# yahs_output_scaffolds_final.agp : AGP格式scaffold组成信息
# yahs_output.bin                 : Hi-C contact信息

# === 生成Juicebox可视化文件 ===
# 将contact信息转换为.hic格式
juicer pre yahs_output.bin \
  yahs_output_scaffolds_final.agp \
  contigs.fasta.fai \
  > alignments_sorted.txt.gz

# 生成.hic文件
java -jar juicer_tools.jar pre \
  alignments_sorted.txt.gz \
  yahs_output.hic \
  <(cat yahs_output_scaffolds_final.fa.fai | awk '{print $1" "$2}')

第四步：3D-DNA流水线¶

代码示例：

# === 安装3D-DNA ===
git clone https://github.com/aidenlab/3d-dna.git
cd 3d-dna

# === 运行3D-DNA ===
# 3D-DNA需要Juicer格式的输入
# 先用Juicer生成merged_nodups.txt

# 方法1: 使用Juicer处理Hi-C数据
# 创建Juicer目录结构
mkdir -p juicer_work/fastq juicer_work/references
cp hic_R1.fq.gz juicer_work/fastq/
cp hic_R2.fq.gz juicer_work/fastq/
cp contigs.fasta juicer_work/references/

# 运行Juicer
juicer.sh -g assembly \
  -d juicer_work/ \
  -D juicer_tools_dir/ \
  -z juicer_work/references/contigs.fasta \
  -p juicer_work/references/contigs.fasta.sizes \
  -t 32

# 方法2: 使用matlock转换BAM到merged_nodups格式
matlock bamfilt -i hic_aligned.bam -o filtered.bam
matlock bam2 juicer filtered.bam merged_nodups.txt

# === 运行3D-DNA scaffolding ===
bash 3d-dna/run-asm-pipeline.sh \
  -r 2 \
  --editor-coarse-resolution 500000 \
  --editor-coarse-region 1000000 \
  --polisher-coarse-resolution 500000 \
  contigs.fasta \
  merged_nodups.txt

# 参数:
# -r 2: 迭代轮数
# 输出: contigs.FINAL.fasta, contigs.FINAL.hic, contigs.FINAL.assembly

第五步：SALSA2 Scaffolding¶

代码示例：

# === 安装SALSA2 ===
conda install -c bioconda salsa2

# === 运行SALSA2 ===
# 输入: contigs + name-sorted BAM + contig长度

# 生成bed文件
bedtools bamtobed -i hic_sorted.bam > hic_links.bed
sort -k4 hic_links.bed > hic_links_sorted.bed

# 运行SALSA2
run_pipeline.py \
  -a contigs.fasta \
  -l contigs.fasta.fai \
  -b hic_links_sorted.bed \
  -o salsa2_output/ \
  -m yes \
  -e GATC,GANTC \
  -p yes

# 参数:
# -a: 输入contigs
# -l: contig长度文件(samtools faidx生成)
# -b: Hi-C链接BED文件
# -e: Hi-C酶切位点(Arima: GATC,GANTC; DpnII: GATC)
# -m yes: 纠正misassembly
# -p yes: 输出contact map

# 输出:
# salsa2_output/scaffolds_FINAL.fasta
# salsa2_output/scaffolds_FINAL.agp

第六步：Juicebox手动校正¶

白话解释： 自动化scaffolding不是完美的——有时会错误地把两条染色体的片段拼在一起，或者把片段方向搞反。Juicebox是一个可视化工具，让你看到contact map上的异常(如不连续的交互块)，然后手动修正。

操作步骤：

# === 生成.hic文件给Juicebox ===
# YaHS输出已经有.hic文件
# 3D-DNA输出也有.hic文件

# 打开Juicebox Assembly Tools (JBAT)
# 下载: https://github.com/aidenlab/Juicebox/releases

# === 手动校正常见操作 ===
# 1. 识别错误拼接: 交互图上两个块之间有"间断"
# 2. 打断(split): 在错误连接处断开
# 3. 移动(move): 将scaffold片段移到正确位置
# 4. 翻转(invert): 翻转方向错误的片段
# 5. 合并(merge): 将应在同一染色体的片段合并

# === 校正后导出 ===
# Juicebox导出.assembly文件
# 使用3D-DNA的post-review脚本生成最终FASTA

bash 3d-dna/run-asm-pipeline-post-review.sh \
  -r contigs.FINAL.review.assembly \
  contigs.fasta \
  merged_nodups.txt

# 或使用juicebox_scripts
python juicebox_scripts/juicebox_assembly_converter.py \
  -a scaffolds.review.assembly \
  -f contigs.fasta \
  -o final_chromosomes.fasta

第七步：质量评估¶

代码示例：

# === 1. BUSCO完整性评估 ===
busco -i final_chromosomes.fasta \
  -l eukaryota_odb10 \
  -o busco_scaffold \
  -m genome \
  -c 16

# === 2. 组装统计 ===
# 使用assembly-stats或QUAST
assembly-stats final_chromosomes.fasta
# N50, L50, scaffold数量, 总长度

# === 3. Hi-C contact map质量 ===
# 好的染色体级组装在contact map上表现为:
# - 清晰的块状对角线(每个块=一条染色体)
# - 块内部交互频率随距离单调递减
# - 块间交互极少

# 生成最终contact map
bwa-mem2 index final_chromosomes.fasta
bwa-mem2 mem -5SP -t 32 final_chromosomes.fasta hic_R1.fq.gz hic_R2.fq.gz | \
  samtools view -bhS - | samtools sort -@ 8 -o final_hic.bam
# 用HiCExplorer生成矩阵
hicBuildMatrix -s final_hic.bam \
  --binSize 500000 \
  -o contact_matrix.cool \
  --threads 16
hicPlotMatrix -m contact_matrix.cool \
  --log1p \
  -o contact_map.png \
  --dpi 300

# === 4. Merqury评估(K-mer based) ===
meryl count k=21 threads=16 reads.fq.gz output reads.meryl
merqury.sh reads.meryl final_chromosomes.fasta merqury_output
# 输出: QV(consensus quality), completeness

实战命令¶

# === 完整流程: contigs → 染色体级基因组 ===

# 1. Hi-C比对
chromap --preset hic -x index -r contigs.fa \
  -1 hic_R1.fq.gz -2 hic_R2.fq.gz \
  --SAM -o hic.sam -t 32
samtools view -bh hic.sam | samtools sort -n -@ 8 -o hic_nsort.bam

# 2. YaHS scaffolding
yahs contigs.fa hic_nsort.bam -o yahs_out -q 10

# 3. 生成Juicebox文件
juicer pre yahs_out.bin yahs_out_scaffolds_final.agp contigs.fa.fai \
  > sorted.txt.gz
java -jar juicer_tools.jar pre sorted.txt.gz output.hic genome.sizes

# 4. Juicebox手动校正 → 导出.assembly

# 5. 最终序列生成
juicer post -o final contigs.fa output.review.assembly

# 6. 质量评估
busco -i final.fa -l eukaryota_odb10 -m genome -c 16
quast final.fa -o quast_output

面试常问点¶

Q1：Hi-C辅助组装的核心原理是什么？¶

A： Hi-C捕获染色质的三维空间邻近性信息。在同一条染色体上的两个位点，其Hi-C交互频率与基因组距离成反比(约1/distance^α)。因此，将contigs排列使得Hi-C contact map呈现平滑衰减的对角线模式，就能确定contigs在染色体上的正确顺序和方向。不同染色体间的交互(trans contacts)远低于染色体内(cis contacts)，这用于区分不同染色体。

Q2：3D-DNA、SALSA2和YaHS各有什么优缺点？¶

A： (1) 3D-DNA: 最成熟，与Juicebox深度集成，自带misassembly校正；但速度慢，流程重(依赖Juicer全流程)。(2) SALSA2: 基于图论方法，内存效率高，错误拼接校正好；但已较少维护，对大基因组不够快。(3) YaHS: 最新最快，内存需求小(适合大基因组)，VGP项目推荐使用；输入简单(只需BAM)。实践推荐：YaHS作为首选，Juicebox校正，必要时与3D-DNA结果对比。

Q3：如何判断Hi-C scaffolding结果的好坏？¶

A：多维度评估：(1) Contact map目视检查：清晰的块状对角线，无明显的off-diagonal异常信号；(2) BUSCO完整性：scaffolding前后应保持不变(scaffolding不该丢基因)；(3) Scaffold数量：应接近物种的已知染色体数(如植物2n信息)；(4) N50/L50：scaffold N50应接近预期染色体大小；(5) QV评分(Merqury)：应高于Q40。

Q4：Juicebox手动校正时常见的错误模式有哪些？¶

A： (1) 误拼(misjoins)：contact map上两个高交互块之间有明显的"断层"(交互突然降为0)——需要在此处打断；(2) 翻转错误(inversions)：一个片段方向反了，表现为沿对角线的"蝴蝶结"图案——需要翻转该片段；(3) 错位(translocations)：属于不同染色体的片段被放在一起——表现为块内有明显的trans信号区域；(4) 遗漏(debris)：很小的scaffold未被整合到主染色体中——尝试手动归位。

Q5：多倍体基因组的Hi-C scaffolding有什么特殊挑战？¶

A：多倍体的同源染色体序列高度相似，导致：(1) Hi-C reads在同源染色体间错误比对，产生假的trans contacts；(2) Scaffolding算法可能错误地将不同亚基因组的contigs拼在一起(嵌合scaffold)；(3) 解决方案：使用ALLHiC(专为多倍体设计)；先用同源基因聚类将contigs分到各亚基因组，再分别scaffolding；提高contig连续性(减少歧义比对)。

易错点¶

1. Hi-C reads比对参数错误¶

错误： 用标准paired-end模式比对Hi-C reads。 正确做法： Hi-C reads的两端通常跨越很远的距离(甚至不同染色体)，不符合paired-end的insert size假设。必须用 -5SP 参数或分别比对两端。

2. 输入contigs含有重复/冗余序列¶

错误： 将包含haplotigs(冗余的等位基因序列)的assembly做scaffolding。 正确做法： 先用purge_dups或purge_haplotigs去除冗余序列，只保留primary assembly。冗余序列会吸收Hi-C reads导致scaffolding质量下降。

3. 不做Juicebox手动校正¶

错误： 直接使用自动化scaffolding的结果。 正确做法： 所有自动化工具都有误判率(通常5-15%的scaffolds需要校正)。Juicebox手动检查是染色体级组装的必要步骤。

4. Hi-C酶切位点指定错误¶

错误： 使用DpnII文库但参数里写了Arima酶切位点。 正确做法： 必须知道Hi-C文库使用的限制酶：DpnII=GATC; MboI=GATC; Arima=GATC,GANTC; Sau3AI=GATC。错误的酶切位点会影响归一化和断点检测。

5. 覆盖度不足¶

错误： 用5x Hi-C数据做scaffolding。 正确做法： Hi-C建议至少30-50x覆盖度(大基因组可能需要更多)。覆盖不足导致远距离contigs间的Hi-C链接太少，无法可靠scaffolding。

补充知识¶

Hi-C数据产生平台¶

平台	酶切位点	特点
Arima Genomics	GATC, GANTC	双酶，覆盖均匀
Dovetail Omni-C	无固定位点(MNase)	不依赖酶切位点
Phase Genomics	多种选择	支持宏基因组Hi-C

典型工作流组合¶

最佳实践 (2024):
  PacBio HiFi → hifiasm → primary contigs
  + Hi-C → YaHS scaffolding → Juicebox校正
  = 染色体级基因组 (contig N50 > 10Mb, scaffold = 染色体数)

数据库和资源¶

VGP (Vertebrate Genomes Project): 参考级基因组组装标准
GenomeArk: VGP数据共享平台
NCBI RefSeq: 高质量参考基因组

计算资源¶

YaHS: 32GB RAM, 16 cores, ~1-2小时(1Gb基因组)
3D-DNA: 64-128GB RAM, 32 cores, ~数小时到1天
Juicebox: 8-16GB RAM的桌面电脑即可

157_Hi-C辅助基因组组装¶

一句话概述¶

核心知识点总览¶

各步骤详解¶

第一步：Hi-C辅助组装的原理¶

第二步：Hi-C数据预处理¶

第三步：YaHS — 快速Hi-C Scaffolding¶

第四步：3D-DNA流水线¶

第五步：SALSA2 Scaffolding¶

第六步：Juicebox手动校正¶

第七步：质量评估¶

实战命令¶

面试常问点¶

Q1：Hi-C辅助组装的核心原理是什么？¶

Q2：3D-DNA、SALSA2和YaHS各有什么优缺点？¶

Q3：如何判断Hi-C scaffolding结果的好坏？¶

Q4：Juicebox手动校正时常见的错误模式有哪些？¶

Q5：多倍体基因组的Hi-C scaffolding有什么特殊挑战？¶

易错点¶

1. Hi-C reads比对参数错误¶

2. 输入contigs含有重复/冗余序列¶

3. 不做Juicebox手动校正¶

4. Hi-C酶切位点指定错误¶

5. 覆盖度不足¶

补充知识¶

Hi-C数据产生平台¶

典型工作流组合¶

数据库和资源¶

计算资源¶

📚 相关文章推荐