Nat Plants 重磅！HapHiC: 对单倍型解析组装进行染色体水平的scaffolding

文摘 2024-08-28 22:34 美国

强烈推荐！有手就行！！使用Hi-C数据在没有参考基因组的情况下对单倍型解析组装进行染色体水平的scaffolding

Zeng, X., Yi, Z., Zhang, X. et al. Chromosome-level scaffolding of haplotype-resolved assemblies using Hi-C data without reference genomes. Nat. Plants (2024). https://doi.org/10.1038/s41477-024-01755-3

Scaffolding对于构建大多数染色体水平的基因组至关重要。高通量染色质构象捕获（Hi-C）技术因其便利性和成本效益已成为主要的scaffolding策略。随着测序技术和组装算法的进步，构建单倍型解析的基因组越来越受到青睐，因为单倍型可以提供额外的等位和非等位变异的遗传信息。ALLHiC是一个广泛使用的等位基因感知的scaffolding工具，专为此目的而设计。然而，它对染色体水平参考基因组的依赖以及较高的染色体错配率仍然阻碍了单倍型解析基因组的解读。在本文中，我们介绍了HapHiC，这是一种不依赖参考基因组的等位基因感知scaffolding工具，在染色体分配以及contig排序和定向方面表现出色。此外，我们通过对各种不利因素进行全面分析，为等位基因感知scaffolding面临的挑战提供了新的见解。最后，在HapHiC的帮助下，我们构建了重要的木质纤维素生物能源作物：巨芒（Miscanthus × giganteus）的单倍型解析三倍体基因组。

构建高质量的参考基因组是物种功能基因组学研究的基础。染色体scaffolding是从头构建真核染色体水平基因组的必要步骤，除非直接组装端粒到端粒的基因组。其目标是确定组装中的contig或scaffold的染色体分配，以及这些序列在染色体上的排序和定向。在早期的基因组研究中，染色体scaffolding通常使用遗传图谱中的连锁群和遗传距离信息来实现。

然而，近年来，高通量染色质构象捕获（Hi-C）技术逐渐取代遗传图谱，成为最广泛使用的染色体scaffolding解决方案，其原因在于其简便、周期短和成本低。Hi-C链接是通过邻近连接和大规模并行测序生成的，用于指示基因组中不同位点之间的染色质相互作用频率。这些信息可以用来推断染色体区域，以及contig或scaffold之间的距离和定向。已经开发了几种基于Hi-C的scaffolding工具，包括LACHESIS、HiRise、3D-DNA、SALSA2和 YaHS，适用于单倍体和单倍体折叠组装。

对于异源二倍体或多倍体，单倍型解析组装由两套或更多套单倍体序列组成。与单倍体折叠组装相比，它提供了额外的遗传信息，如双等位基因或多等位基因，以及非等位变异之间的顺反配置。最近在测序技术和组装算法方面的进展推动了单倍型解析基因组的解读。太平洋生物科学公司（PacBio）的HiFi测序和牛津纳米孔技术公司（ONT）的双链测序都达到了Q30（99.9%）的基础准确度水平，这为更准确的等位基因分相提供了坚实基础。后代数据使用来自父母基因组的短读长序列进行分相，从而实现全基因组水平的分相。最近，hifiasm利用Hi-C测序数据进行染色体水平的分相，无需父母数据。这两种方法在处理二倍体或类二倍体异源多倍体基因组方面表现出高准确性。因此，后续的染色体scaffolding可以独立地在每个相位的单倍型上进行。

同源多倍体在种子植物中普遍存在，特别是在重要经济作物中。在同源多倍体基因组中进行单倍型分相有助于研究作物驯化历史和遗传育种。它也为分析全基因组重复后的等位基因表达优势和基因组进化奠定了基础。然而，组装单倍型解析的同源多倍体基因组比二倍体基因组更具挑战性。后代数据不适用于同源多倍体，而hifiasm的基于Hi-C的算法在同源多倍体基因组组装中产生了不平衡的分相结果。因此，构建单倍型解析同源多倍体基因组的最常见策略是执行等位基因感知的scaffolding，即在染色体scaffolding期间同时使用Hi-C数据将contig分配到不同的单倍型。此外，分别对每个分相的单倍型进行scaffolding可能会导致错误，因为来自多个单倍型的Hi-C数据比对到单个单倍型，而忽略了单倍型之间可能的染色体结构变异。这再次强调了等位基因感知scaffolding的重要性。

ALLHiC是一种广泛使用的Hi-C scaffolding工具，专为等位基因感知scaffolding设计。它有效地识别等位基因序列，并在聚类前移除它们之间的Hi-C链接以减少干扰。ALLHiC在单倍型分相方面表现出色，并已应用于解决多个重要作物的二倍体和自四倍体基因组。然而，这种方法需要来自近缘物种的染色体水平参考基因组，这对于许多物种可能无法获得。尽管可以将单倍型折叠基因组组装和注释为参考，但这显著增加了基因组研究的时间和成本。此外，据观察，使用参考基因组时，ALLHiC 会引入聚类错误。这些限制和缺点在一定程度上阻碍了单倍型解析基因组的构建，特别是在同源多倍体中。

在本文中，我们介绍了HapHiC，这是一种基于Hi-C的scaffolding工具，使得在没有参考基因组的情况下能够进行同源多倍体组装的等位基因感知染色体scaffolding。我们对可能阻碍基因组等位基因感知scaffolding的因素进行了全面分析。与现有方法相比，HapHiC在解决这些挑战时表现出更高的scaffolding连贯性和更低的错误分配率。此外，HapHiC速度快且资源效率高，并已成功验证于不同倍性和分类群的基因组中。利用HapHiC，我们最终构建了重要的木质纤维素生物能源作物巨芒（Miscanthus × giganteus）的单倍型解析基因组。

使用Hi-C数据实现染色体级别单倍型的从头scaffolding

为了克服使用Hi-C数据进行染色体级别单倍型scaffolding的关键挑战，我们开发了HapHiC，这是一种Hi-C scaffolding工具，在处理无参考基因组的单倍型解析组装方面表现出色。我们使用HapHiC构建了复杂的三倍体巨芒基因组。

问题

构建单倍型解析基因组是优先选择的，因为单倍型比折叠基因组提供了更多的等位和非等位变异的信息。这些信息对于诸如等位基因特异性基因表达研究等研究至关重要。目前，高通量染色质构象捕获（Hi-C）技术是最广泛使用的单倍型分相和scaffolding方法。然而，这些任务对于在植物中常见的同源多倍体特别具有挑战性。ALLHiC是一种为执行同源多倍体单倍型分相和scaffolding而设计的算法，但它需要一个注释良好的染色体级别参考基因组。这种资源可能不可用，或可能导致额外的构建成本。此外，依赖参考基因组可能引入scaffolding错误。因此，为了加速基因组的解读，需要一种不依赖参考（即从头）的Hi-C scaffolding工具，它支持不同的倍性水平，并在精度上优于当前工具。

解决方案

单倍型scaffolding中的主要挑战来自等位contig对之间的强Hi-C信号、来自不同单倍型的折叠contig的错误合并以及来自错误连接的嵌合contig。因此，需要一种参考独立的方法来识别等位contig，并增加对各种组装错误的容忍度。除了区分染色体来源并确定contig的排序和定向外，Hi-C数据还提供了宝贵的信息来解决这些问题。我们通过一系列算法创新利用这些信息来识别组装错误和等位contig对的独特Hi-C信号模式，与正确组装的contig相比。染色质信号和错误组装的contig在染色体聚类之前被暂时从数据中移除，以提高对这些问题的容忍度。这些创新一起促成了HapHiC，我们的参考独立的Hi-C scaffolding工具，能够进行染色体级别单倍型的scaffolding（图1）。我们通过对可能阻碍单倍型scaffolding的各种不利因素进行模拟测试，以及在广泛的物种真实基因组组装中的验证，评估了这些因素的实际影响和HapHiC、ALLHiC以及其他软件工具（LACHESIS、3D-DNA、SALSA2和YaHS）的表现。

图1: HapHiC概述。 在染色体分配之前，contig及其连接的Hi-C序列（称为Hi-C链接）经历预处理步骤（黑色交叉表示移除等位contig对之间的Hi-C链接；黑色箭头表示组装校正的断点)。不利的Hi-C信号和错误组装的contig在Markov聚类之前被暂时移除。过滤出的contig然后重新分配到最适合的聚类中。如果聚类数超过预期的染色体数，则进行额外的层次聚类。在此之后，通过我们称之为“快速排序”的加速3D-DNA迭代scaffolding算法对每条染色体的contig进行排序和定向。随后，使用ALLHiC中使用的遗传算法和贪婪方法进一步优化结果。

嵌合contig通常是由于不同染色体或基因组区域之间的错误组装导致的。这意味着错误组装的contig可能同时与实际空间距离很远的不同基因组区域显示强烈的Hi-C信号。这种不一致通常反映在错误组装的contig的邻域密度较低，这使得HapHiC能够将其与正常contig区分开来。此外，等位contig对之间的Hi-C信号与序列比对一致，并沿着基因组位置呈对角线分布。这种模式不同于正常contig之间的Hi-C信号。通过利用这一特征，HapHiC可以识别等位contig对。此外，我们的观察和实际项目及模拟实验的证据表明，等位contig对之间的强信号主要来源于基因组组装过程中碱基水平的切换错误。

影响

本研究中实施的模拟各种不利因素的综合流程为未来的scaffolding工具开发奠定了基础。我们观察到等位contig对之间对角分布的Hi-C信号主要是由组装错误引起的，这为设计新的度量标准以评估组装质量提供了见解。我们的验证实验表明，HapHiC在广泛的分类群和倍性水平中表现出色。此外，HapHiC兼容最近开发的测序技术Pore-C。

尽管HapHiC是一个强大的工具，但它也有局限性。为了获得更好的结果，它需要预先了解生物体的染色体数量。此外，组装质量较差的基因组（尤其是那些包含大面积折叠区域的基因组）即使在HapHiC下也仍然对准确的scaffolding构成挑战，并且需要手动调整。

随着测序技术和组装算法的进步，与单倍型scaffolding相关的挑战预计也会发展。例如，大的折叠区域可能会组装成多个几乎相同的区域，从而复杂化Hi-C和Pore-C序列对齐的准确性。我们未来对HapHiC的开发将继续以问题为导向，并与测序技术及上游组装和对齐方法的进步保持同步。--Xiaofei Zeng & Guoan Chen 南方科技大学

专家意见

“HapHiC在Hi-C scaffolding的几个关键领域进行了改进，包括处理嵌合contig断裂、contig聚类和等位链接移除的工程改进等。这是一篇非常写得很好的方法论文，是对从头组装的研究人员以及特别是多倍体植物基因组工作者的Hi-C工具箱的及时更新。” --唐海宝，福建农林大学

论文背后的故事

在开发HapHiC期间，我们遇到了几个障碍。一个挑战是在组装校正后缓解剩余嵌合contig对scaffolding的负面影响。鉴于诸如Hi-C序列覆盖率等信息已经在组装校正中使用，识别这些contig的额外特征变得困难。当我们将关注点从contig内部的信息转移到contig之间的关系网络时，一个转折点出现了。我们发现可以通过邻域密度识别潜在的嵌合contig，这在检测折叠contig方面也非常有效。令人惊讶的是，这种方法补充了传统使用Hi-C链接密度的方法。这一算法的集成使HapHiC在处理复杂的同源多倍体基因组（如甘蔗基因组）时表现出色。我们克服的每一个挑战都像是解决一个难题，每一次突破都进一步激发了我们的兴奋和动力。X.Z. & G.C.

来自编辑的意见

“现有的单倍型分相基因组scaffolding工具需要来自近缘物种的现有参考基因组。这项研究的突出之处在于它开发了一种不需要参考基因组的等位基因感知染色体scaffolding工具。这对生成单倍型分相基因组组装尤其有用，特别是对于难以分相的同源多倍体植物基因组。” Jun Lyu，资深编辑， 《自然植物》.

Fig. 2: Comprehensive performance analysis of Hi-C-based scaffolding tools in chromosome assignment under various adverse conditions.

null — Fig. 3: Evaluation of Hi-C-based scaffolding tools’ performance in contig ordering and orientation across assemblies with varying contig N50 values.

Fig. 5: Comparative analysis and examples of HapHiC in scaffolding published autotetraploid genomes.

Fig. 6: Comparative genomic analysis of M. × giganteus and other Miscanthus species.

Comparison of inter-allele Hi-C links between hifiasm and HiCanu assemblies.

Installation

# (1) Download HapHiC from GitHub
git clone https://github.com/zengxiaofei/HapHiC.git

# (2) Resolve dependencies
conda env create -f HapHiC/conda_env/environment_py310.yml
# Activate the HapHiC conda environment
conda activate haphic 

# (3) Check whether all dependencies are correctly installed
/path/to/HapHiC/haphic check

# (4) Show all available commands and help message
$ /path/to/HapHiC/haphic -h

Align Hi-C data to the assembly

# (1) Align Hi-C data to the assembly, remove PCR duplicates and filter out secondary and supplementary alignments
bwa index asm.fa
bwa mem -5SP -t 28 asm.fa /path/to/hic_read1_fq.gz /path/to/hic_read2_fq.gz | samblaster | samtools view - -@ 14 -S -h -b -F 3340 -o HiC.bam

# (2) Filter the alignments with MAPQ 1 (mapping quality ≥ 1) and NM 3 (edit distance < 3)
$ /path/to/HapHiC/utils/filter_bam HiC.bam 1 --nm 3 --threads 14 | samtools view - -b -@ 14 -o HiC.filtered.bam

温馨提示

这里的asm.fa可以是单倍体折叠的contigs（例如，hifiasm中的p_ctg）、单倍型分相的unitigs（例如，hifiasm中的p_utg）或一个或多个单倍型的contigs集合（例如，hifiasm中的hap*.p_utg）。此外，asm.fa也可以是其他软件输出的scaffolds。

您可以根据自己的喜好或要求准备BAM文件，但请勿按坐标排序。如果您的BAM文件已经按坐标排序，您需要按序列名称重新排序（samtools sort -n）。

我们不推荐使用Juicer流程进行Hi-C序列比对，尤其是在单倍型分相组装中。

Run HapHiC scaffolding pipeline

单行命令。HapHiC提供了一个单行命令haphic pipeline来执行整个scaffolding流程。所需参数包括：

- asm.fa：您的基因组组装文件，格式为FASTA
- HiC.filtered.bam：在上一步中准备的BAM文件
- nchrs：组装中存在的染色体数量，同时也是预期的输出scaffolds数量。

/path/to/HapHiC/haphic pipeline asm.fa HiC.filtered.bam nchrs

Work with hifiasm

当为分型的hifiasm组装进行scaffolding时，您可以使用hifiasm输出的GFA文件运行HapHiC。这里，“分型的hifiasm组装”指的是通过后代数据或基于Hi-C的算法组装的单倍型解析主链（.hap.p_ctg.gfa）以及分型的unitigs（p_utg.gfa）。

HapHiC使用GFA文件中的序列深度信息在聚类之前过滤掉潜在的折叠contigs/unitigs。如果提供了多个GFA文件，HapHiC假设这些GFA文件是特定于单倍型的（.hap.p_ctg.gfa），并根据此分型信息人工移除或减少单倍型之间的Hi-C连接。请注意，GFA文件中的contigs/unitigs应与FASTA文件中的一致。无论是.gfa还是noseq.gfa文件都是可以接受的。

# (1) For hifiasm primary unitigs, use the GFA file to filter out potential collapsed unitigs before clustering
/path/to/HapHiC/haphic pipeline p_utg.fa HiC.filtered.bam nchrs --gfa p_utg.gfa

# (2) In addition to read depth filtering, HapHiC can also reduce Hi-C links between haplotypes according to phasing information in GFA files for haplotype-resolved primary contigs

# By default, all Hi-C links between haplotypes are completely removed and contigs from different haplotypes will not be clustered into the same group
/path/to/HapHiC/haphic pipeline allhaps.fa HiC.filtered.bam nchrs --gfa "hap1.p_ctg.gfa,hap2.p_ctg.gfa"

# The weight can be set to 0 to ignore the phasing information. You can also set it between 0 and 1 to run HapHiC as a double check. In the latter case, contigs from different haplotypes might be clustered together
/path/to/HapHiC/haphic pipeline allhaps.fa HiC.filtered.bam nchrs --gfa "hap1.p_ctg.gfa,hap2.p_ctg.gfa" ---phasing_weight 0

YaHS-style `scaffolds.raw.agp` (recommended)

为了避免重新比对Hi-C数据，从HapHiC版本1.0.1开始，我们引入了类似YaHS的scaffolds.raw.agp文件。在这个AGP文件中，断开的contigs不会被分配新的ID。相反，其在原始contigs中的起始和结束坐标显示在第七和第八列。通过遵循YaHS提供的方法，您可以生成.assembly和.hic文件，而无需重新比对。

在构建最终的scaffolds后，HapHiC会自动生成一个用于在Juicebox中进行可视化和整理的shell脚本。确保在您的系统上已安装并添加Java和samtools到$PATH中。然后，运行以下命令：

# bash, not sh
bash juicebox.sh

在输出日志文件out_JBAT.log中，您可以找到相应的比例因子，例如，[I::main_pre] scale factor: 2。为了确保Hi-C互作图与Juicebox中的scaffolds和超scaffolds边界正确比对，请在Juicebox的菜单中通过 Assembly > Set Scale 设置您自己的比例因子。

对于大基因组，有必要在juicebox.sh文件中调整Java的内存设置（例如，设置为-Xmx64G或更高），以避免内存不足错误或提高执行速度。

在Juicebox中手动整理后，您将获得修改后的组装文件out_JBAT.review.assembly。要生成最终的scaffolds的FASTA文件：

# Generate the final FASTA file for the scaffolds
/path/to/HapHiC/utils/juicer post -o out_JBAT out_JBAT.review.assembly out_JBAT.liftover.agp asm.fa

Visualization

自HapHiC版本1.0.2以来，我们引入了一个haphic plot命令，用于生成高度可定制的Hi-C互作图。此命令需要两个输入文件：过滤后的BAM文件HiC.filtered.bam和一个包含与BAM文件中contig ID相匹配的scaffold AGP文件。

# For HapHiC scaffolding result
/path/to/HapHiC/haphic plot scaffolds.raw.agp HiC.filtered.bam
# For the AGP file generated after manual curation in Juicebox
/path/to/HapHiC/haphic plot out_JBAT.FINAL.agp HiC.filtered.bam

可视化的Hi-C互作图将输出为contact_map.pdf。如果BAM文件较大，此过程可能会有些缓慢，每10 GiB的BAM文件大约需要几分钟的时间。完成后，程序会生成一个名为contact_matrix.pkl的二进制文件。此文件可以替代HiC.filtered.bam用于更快速的可视化（大约1分钟）。

# Use previously generated `contact_matrix.pkl` for faster visualization
/path/to/HapHiC/haphic plot out_JBAT.FINAL.agp contact_matrix.pkl

此外，您还可以创建一个名为separate_plots.pdf的文件，该文件分别展示了每个scaffold的互作图：

/path/to/HapHiC/haphic plot out_JBAT.FINAL.agp HiC.filtered.bam --separate_plots

# 或者 输入上一步生成的contact_matrix.pkl
/path/to/HapHiC/haphic plot out_JBAT.FINAL.agp contact_matrix.pkl --separate_plots

Order and orient whole scaffolds using a reference genome

HapHiC在版本1.0.4中引入了一个单独的命令haphic refsort，用于根据参考基因组对整个scaffolds进行排序和定向。

首先，您应使用minimap2^[1]将原始contigs（不是scaffolds）与参考基因组比对，生成PAF文件。参考基因组可以来自同一物种或近缘的物种：

# The preset can be `asm5` if the reference genome is well-assembled from the same species
minimap2 -x asm20 ref.fa asm.fa --secondary=no -t 28 -o asm_to_ref.paf
# `haphic refsort` can also be compatible with other aligners, like wfmash
wfmash ref.fa asm.fa -m -n 1 -S 1 -t 28 | cut -f 1-6,8- > asm_to_ref.paf

通过使用haphic refsort命令，您可以基于PAF文件生成一个新的AGP文件：

# By default, scaffolds are output based on the alphabetical order of the chromosome IDs of the reference genome
haphic refsort 04.build/scaffolds.raw.agp asm_to_ref.paf > scaffolds.refsort.agp
# You can manually specify the order by listing them and separating with commas
haphic refsort 04.build/scaffolds.raw.agp asm_to_ref.paf --ref_order chr1,chr2,chr3,chr4,... > scaffolds.refsort.agp

生成的scaffolds.refsort.agp文件可以直接用于Juicebox纠错^[2]和haphic plot的可视化^[3]。请注意，此功能不会更改您的scaffolds，它只通过整体排序和定向来改变展示方式。

以下是一个关于四倍体甘蔗Np-X组装的示例：

Received

13 December 2023

Accepted

01 July 2024

Published

05 August 2024

引用链接

[1] minimap2: https://github.com/lh3/minimap2
[2] Juicebox纠错: https://github.com/zengxiaofei/HapHiC?tab=readme-ov-file#juicebox
[3] 可视化: https://github.com/zengxiaofei/HapHiC?tab=readme-ov-file#visualization

http://mp.weixin.qq.com/s?__biz=MzkzMTM0NTcyNw==&mid=2247487558&idx=2&sn=1c8d0e91fd5f797d9c39102e9e8b2572

X Omics

生物科学与计算机科学的完美碰撞，激发出探索世界的全新视角，让我们一起探索生命科学的新纪元！合作交流：xomics1@gmail.com

Nat Plants 重磅！​HapHiC: 对单倍型解析组装进行染色体水平的scaffolding