几种单倍型基因组组装方法的比较

发表于 2022-07-12 更新于 2022-07-13 分类于 Learn ， Bioinfo 阅读次数： 1480 本文字数： 31k 阅读时长 ≈ 28 分钟

摘要

二倍体个体具有两套遗传信息，一套来自父本，一套来自母本。在大多数二倍体基因组组装中，来自同源染色体的两个同源拷贝被折叠在一起，最终得到一个马赛克 (mosaic) 序列，也被称为伪单倍型 (pseudo-haplotypes)，即一套基因组表示两个单倍型信息，缺失了近 50% 的等位变异信息。单倍型基因组信息对于研究基因组如何影响表型差异至关重要。常规二倍体基因组组装缺失的单倍型信息将会影响后续基因注释的准确性，忽略了两个同源染色体之间的差异。本研究介绍了 Gamete binning、FALCON-Phase、Trio binning 和 ALLHiC 4 种用于单倍型基因组组装的方法，旨在获得二倍体两个完整的单倍型基因组，并从连续性、完整性、交换错误率等方面评估其单倍型组装质量。分析整理 4 种组装方法获得不同水平的单倍型基因组，根据不同的测序数据选择合适的单倍型组装方法，比较得出更适合单倍型基因组组装的组装方法，发现 Gamete binning 的整体效果最佳，可以获得染色体水平的单倍型解析基因组，N50 达到了 25Mbp。 ALLHiC 虽然可以得到染色体水平的组装，但对于二倍体基因组只获得了单倍型混合基因组，并没有将两个单倍型完全解析出来。而 FALCON-Phase 和 Trio binning 虽然获得了两个单倍型，但两者组装序列连续性较低，N50 均没有达到兆 (Mb) 级。整理 4 种组装方法的优缺点，尝试寻找一种适合单倍型基因组组装的组装方法，为后续的基因注释等下游分析提供保障。得到两个完整的单倍型基因组序列，有助于了解单倍型之间的等位基因差异影响个体表型差异的机制。

关键字：二倍体；单倍型基因组组装；方法评估；等位变异；染色体水平组装；伪单倍型

1 前言

1.1 项目背景

单倍型 (Haplotype) 是指共存于单条染色体上的一系列遗传变异位点的组合（李双双等 2018），是一组来自相同亲本的等位序列组合，描述的是一段单条染色体上的序列差异。单倍型基因组组装旨在获得单倍型解析 (haplotype-resolved) 的组装序列，即从二倍体物种获得两个完整的单倍型基因组序列；单倍型基因组之间存在的序列差异可能导致某些表型或性状发生改变，如等位基因表达差异、等位基因甲基化差异等。

大多数二倍体基因组组装工具只是简单的将读段 (reads) 组装到重叠群 (contigs) 水平，忽略了同源染色体之间的差异，得到的组装序列被称为伪单倍型 (pseudo-haplotypes) 序列，其中混合了两个单倍型的序列。这种单倍型混合的组装可能影响下游分析 (Al Bkhetan et al 2021)，从而导致错误的生物学解释。二倍体基因组的两组同源染色体，一组来自父本，另一组来自母本；彼此之间存在复杂的等位基因变异。得到两个完整的单倍型基因组序列，有助于了解单倍型之间的等位基因差异影响个体表型差异的机制，如最近在荔枝 (Litchi chinensis Sonn) 单倍型基因组中发现每个单倍型之间存在显著的覆盖差异，具有特定于极早熟和晚熟的不同模式 (Edger 2022, Hu et al 2022)，这一研究证实了单倍型基因组之间的等位差异会影响荔枝的果实成熟期。完全相同的一组变体的不同配置可能会导致表型和疾病易感性的不同结果，这往往是通过将组装序列对齐人类参考基因组获得，所以精确的人类单倍型解析基因组对于了解基因组与个体表型之间的关系至关重要 (Cao et al 2015)。

为了得到高精度的单倍型基因组序列，一方面可以从测序样本入手，另一方面则是开发新的组装算法，从杂合二倍体中尽可能的恢复其两个单倍体基因组序列。从测序样本入手，可以考虑近交、染色体分选 (chromosome sorting)、Strand-seq。对目标物种进行近交，减少单倍型变异，获得纯合个体的基因组组装，但这限制于可育；况且，人为干预得到的纯合个体可能导致基因组不能代表自然种群中发现的变异 (Koren et al 2018)。染色体分选 (Yang et al 2011) 是通过 Phase-Seq (Chen et al 2017) 技术分离同源染色体的两条姐妹染色单体，并分别进行测序组装，从而获得染色体层次的单倍型基因组，但这依赖于染色体分选时的荧光强度，并不能保证染色体的两条姐妹染色单体被正确分离。Strand-seq 是一种单细胞测序技术，使用二代测序技术独立分析每条染色体的姐妹染色单体遗传模式 (Falconer and Lansdorp 2013)；仅对 DNA 模板链进行测序，保留了单个同源物的结构连续性 (Kronenberg et al 2021)，但是其数据的生成难度限制了其广泛应用。

FALCON-Unzip 无法完整的输出单倍型序列，只能进行部分定相，从主要重叠群 (primary contigs) 中识别杂合变异，并根据杂合位置之间的相位信息输出一个单倍型参考基因组序列和单倍型特异性重叠群 (haplotigs)(Chin et al 2016)。FALCON-Unzip 输出的单倍型参考基因组序列混合了两个单倍型的序列，可能会影响下游的基因注释准确性；输出的单倍型特异性重叠群也只是杂合区域的局部单倍型信息。FALCON-Phase 作为 FALCON-Unzip 的扩展，利用 Hi-C 测序数据的超远程连接信息将相位块（phase block）连接，从而生成更长的单倍型序列，但并不能定相整个染色体。

DipAsm 可以利用 PacBio HiFi 测序数据和 Hi-C 数据重建二倍体的两个单倍型基因组 (Garg et al 2020)。在人类基因组中得到了验证，但由于测试物种杏 (Prunus armeniaca) 没有公开的 PacBio HiFi 测序数据，故本文没有将其一同比较。此外，能够利用 PacBio HiFi 数据和 Hi-C 数据的单倍型组装算法还有 hifiasm，不过 hifiasm 需要两个亲本的 Hi-C 或 Illumina 短读段双末端测序数据和子代的 PacBio HiFi 测序数据 (Cheng et al 2021)，这也限制了 hifiasm 的发展。相比之下，DipAsm 在不需要亲本测序数据的情况下，尽可能的恢复二倍体生物的单倍型解析基因组。PacBio HiFi 测序数据单分子读段的准确度高达 99.9%(Wenger et al 2019)，在保证长读段的同时兼顾了测序准确性，可以获得连续性更好、适合结构变异检测的基因组组装。

本文比较了 4 种单倍型基因组组装方法（除 ALLHiC）产生的二倍体两组单倍型基因组序列。准确的单倍型基因组信息为研究同源染色体等位差异提供了可能性。

1.2 Trio binning 方法概述

Trio binning (Koren et al 2018) 在组装之前，利用亲本基因组之间的基因组差异将子代全基因组测序读段分离成单倍型特异性读段集，并分别进行基因组组装，从而获得目标个体的两个单倍型基因组。二倍体个体分别从两个亲本各继承一个重组单倍型。二倍体基因组存在的复杂等位基因变异阻碍了其单倍型基因组序列的组装，而 Trio binning 在组装前解决等位基因变异来简化单倍型组装，即在组装前识别来自父母 Illumina 短读段测序数据的特异性 k-mers，将子代 PacBio CLR 长读段划分为特定于单倍型的集合，而后对每个单倍型进行独立组装，从而形成完整的单倍型解析基因组，避免了单倍型之间等位变异的干扰 (Koren et al 2018)。但这取决于目标个体的杂合性，后代杂合度越高，单倍体分型效果越佳。Trio binning 方法集成于 Canu 程序，可在 Canu 程序中直接调用。

1.3 ALLHiC 方法概述

ALLHiC (Zhang et al 2019) 是专为多倍体构建单倍型、染色体级别组装设计的一种新算法；该算法能够使用 Hi-C 配对末端读段和创新的修剪 (prune) 和优化 (optimize) 步骤构建多倍体基因组的等位基因感知、染色体规模的单倍型基因组组装 (Zhang et al 2019)。ALLHiC 能够实现同源多倍体基因组的染色体水平的从头组装，分离每个等位基因，同时也适用于异源多倍体和杂合二倍体基因组的单倍型染色体水平组装。ALLHiC 算法包括五个步骤：修剪 (pruning)，分区 (partition)，修复 (rescue)，优化 (optimization) 和构建 (building)(Zhang et al 2019)。修剪步骤提供等位基因感知的组装技术；分区步骤使用聚类算法将重叠群 (contigs) 划分为与染色体条数相对应数量（用户定义的 k 值）的组；修复步骤是一个检查的过程，搜索最佳 Hi-C 信号，在初始修剪步骤中由于缺乏有效信号而未处理的重叠群被进一步划分到特定组中；优化步骤使用遗传算法对每个组内的重叠群进行排序和定向；最终从重叠群中构建染色体水平或支架 (scaffold) 水平组装。

1.4 Gamete binning 方法概述

Gamete binning (Campoy et al 2020) 是一种基于单倍体配子的单细胞测序方法，能够将全基因组测序读段分离成单倍型特异性读段集；在组装每个单倍型的读段后，使用源自配子的遗传图谱将 contigs 构建到染色体水平。Gamete binning 首先从目标个体中分离配子核，对数百个单倍体配子基因组进行高通量单细胞测序；配子基因组中序列变异的分离能够直接将所有变体定相为两个单倍型，并给予配子基因组中的重组模式构建基因图谱。基于遗传连锁群利用此外定相好的等位变异分离子代 PacBio CLR 长读段为两个不同的读段组（每组代表一个单倍型），而后对每个单倍型进行独立组装，并使用配子衍生的遗传图谱将每个单倍型支架 (scaffolding) 到染色体水平 (Campoy et al 2020)。

1.5 FALCON-Phase 方法概述

FALCON-Phase (Kronenberg et al 2021) 是一个定相 (phasing) 工具，使用超远程 Hi-C 染色质作用数据将部分定相 (partially phased) 二倍体组装的相位块 (phase block) 扩展到染色体或支架 (scaffold) 水平。PacBio CLR 长读段测序组装一般得到的是主要重叠群（同时代表两个单倍型信息）和较短的交替单倍体 (haplotigs)。FALCON-Phase 输入一个部分定相的长读段重叠群水平的组装 (primary contigs)，并使用来自同一样本的 Hi-C 测序数据扩展重叠群的定相。这里引入一个新的概念–相位块 (phase block)，通过将交替的单倍体与其相关的主要重叠群对齐来定义。通过相位块的对齐区分主要重叠群的定相区域与非定相区域，然后使用 Hi-C 测序数据对沿着每个重叠群处于相同相位（相同亲本同源物）的单倍型块进行分类，最后通过将折叠的序列整合到两个单倍型中来扩展组装单倍型序列，以获得两组重叠群集。FALCON-Phase 的分型效果取决于目标个体的杂合度和提供的主要重叠群准确度，杂合度越高，分型效果越佳。

1.6 研究目的

分别利用 Trio binning、ALLHiC、Gamete binning、FALCON-Phase 4 种组装方法获得杏 (Prunus armeniaca) 的单倍型基因组序列，并从连续性、完整性、准确性等方面评估 4 种组装方法的单倍型组装质量。根据不同的测序数据选择合适的单倍型基因组组装方法，获得两个完整的单倍型基因组序列。整理 4 种组装方法的优缺点，尝试寻找一种适合单倍型基因组组装的组装方法，为后续的基因注释等下游分析提供保障。

2 材料与方法

2.1 原始测序数据统计

根据 4 种组装方法所需测序数据的要求，我们采用 Prunus armeniaca（杏，栽培品种：“Rojo Pasión”）作为评估组装效果的物种，所有的原始测序数据均可从 ENA 数据库中获得，登录号为 “PRJEB37669”，其中 Currot 和 OrangeRed 为目标个体的两个亲本。对原始测序数据进行了初步统计（表 1），包含两个亲本高于 60x 的 Illumina 双端测序，这将用于 Trio binning 分型以及后期验证各种方法分型质量；目标个体 F1 接近 80x 的 PacBio CLR 长读段测序 (19.9 Gbp)、Hi-C (57.6 Gbp) 以及花粉 10x-Genomics 单细胞测序数据 (61.7 Gbp)。

表 1 原始测序数据统计
Table 1 The statistics of raw sequencing data

数据来源	数据类型	读段 / 读对 (条 / 对)	总碱基数 /bp	测序深度
F1 果实组织	PacBio CLR	2,587,118	19,929,106,572	~ 80x
F1 叶片组织	Hi-C	191,194,281	57,579,855,538	~ 230x
F1 花粉	10x Genomics	220,312,318	61,687,449,040	~ 247x
Currot 叶片组织	Illumina	53,802,680	16,226,087,299	~ 65x
OrangeRed 叶片组织	Illumina	52,115,882	15,720,264,273	~ 63x

2.2 基于 Kmer 的基因组特征分析

在进行基因组组装之前，了解基因组的基本特征有助于我们后续组装参数的调整。GenomeScope v2.0 可以在没有参考基因组的情况下，分析原始测序读段的 k-mer 频率以估计目标基因组的主要特征，如基因组大小、杂合度和重复率 (Ranallo-Benavidez et al 2020)。首先，利用 Jellyfish v2.2.10 (Marçais and Kingsford 2011) 工具以 kmer=21 统计来自 10x-Genomics Illumina 短读段测序数据 k-mers 的出现频率，并将其输入到 GenomeScope v2.0 程序，参数设置为 “–max_kmercov=1000000”（过滤高频出现的 k-mers）；这将会过滤在测序数据中出现次数超过 1000000 次的 k-mers，从而避免人为引入的碱基错误。从 k-mers 频率分布图，我们得知预测的单倍型基因组大小约为 250 Mbp，杂合率为 1.06%，重复率为 41%（图 1 a）。在 coverage=1 的位置出现大量唯一 k-mer（红色线）可能是 PCR 错误导致。

图 1 k-mers 频率分布图与 PacBio 长读段长度分布
Fig. 1 Frequency distribution of k-mers and PacBio long read length distribution
a) 为 K-mers 频率分布图，x 轴表示 k-mer 出现次数，y 轴则指示出现 x 次的 k-mer 数量。从左至右依次为杂合峰、纯合峰、重复峰。其中 len 为预测的单倍型基因组大小；uniq 表示基因组中非重复序列比例；aa 为纯合率；ab 为杂合率。b) PacBio CLR 测序读段 (19.9 Gbp) 的长度分布。

2.3 二倍体基因组初步组装

对 F1 的 19.9 Gbp 原始 PacBio CLR 测序数据 (图 1 b)，使用 Canu (Koren et al 2017) v2.2 构建初步的二倍体组装。参数设置为 “‘batOptions=-dg 3 -db 3 -dr 1 -ca 500 -cp 50’ genomeSize=250m corOutCoverage=200”，其中参数‘batOptions=-dg 3 -db 3 -dr 1 -ca 500 -cp 50’是针对多倍体基因组设计的优化参数，尽可能的组装杂合基因组序列，但这也适用于杂合度高的杂合二倍体组装。Genomescope2 预测 F1 基因组杂合率为 1.06%，这说明其基因组组装有一定困难。最终，得到约为 389 Mbp 的基因组（二倍体基因组的单倍型混合表示），是预测单倍型基因组大小 (250 Mbp) 的 1.56 倍，其中含有将近 52.4% 的重复序列。

基因组的某些部分可能具有非常高的杂合性，导致该部分两种单倍型的重叠群被组装为单独的主要重叠群 (primary contigs)，而不是作为重叠群或相关的单倍体。Canu 获得的基因组是二倍体的单倍型融合表示 (Roach et al 2018)(haplotype-fused representation)，即代表两组等位基因序列。其算法会将高度杂合的区域组装成单独的重叠群（对于二倍体来讲，可以得到两条染色单体），而不是预期的单个单倍型融合的重叠群，从而导致基因组大小显著大于单倍型基因组的实际大小。Purge haplotigs 程序可以为单倍型组装生成去冗余的单倍型融合重叠群 (Roach et al 2018)。最终 389 Mbp 的基因组去除 166 Mbp 单倍型特异性重叠群 (haplotigs)，得到约 219 Mbp 的单倍型融合组装（其中丢失了约 3.4 Mbp 可能存在的错误），与 Genomescope2 预测的单倍型基因组大小相差不大。

2.4 基于 Trio binning 方法获得单倍型组装序列

首先使用 Meryl 对亲本近 60x 的未组装短读段测序数据进行 k-mers (kmer=19) 统计以获得 k-mers 频率分布，并去除可能存在的错误 k-mers（低拷贝）和重复序列（高拷贝），只留下来自特异性纯合或杂合基因组序列的 k-mers (Koren et al 2018)。接着，对 F1 代 PacBio CLR 测序数据中的每一条长读段进行 k-mers 统计，并与亲本特异性 k-mers 数据库进行比对，将长读段分给具有最匹配的单倍型特异性 k-mer 的单倍型读段集。最后分别利用 Canu v2.2 组装工具对两个单倍型进行组装，从而得到两个重叠群 (contigs) 水平的单倍型分辨的单倍型基因组。

2.5 基于 ALLHiC 方法获得染色体水平的单倍型基因组序列

将 Canu 得到的单倍型融合基因组和 Hi-C 数据输入 ALLHiC 程序，参数设置为 “-k 16 -e GATC”，最终获得一个含有 16 条染色体以及未识别的 scaffolds 的 fasta 文件 (groups.asm.fasta)。对于简单二倍体，Zhang 等人 (Zhang et al 2019) 提供了一个流程脚本 (https://github.com/tangerzhang/ALLHiC/blob/master/bin/ALLHiC_pip.sh) 可直接对基因组进行支架 (scaffolding)。为了验证 ALLHiC 管道是否能够分离同源染色体，以参数 “-k 8 -e GATC” 运行 ALLHiC 管道，对 Purge haplotigs 处理的主要重叠群进行支架 (scaffolding)，获得 8 条染色体的基因组序列。

2.6 基于 Gamete binning 方法获得染色体水平的单倍型组装序列

两个完整的单倍型基因组可从 ENA 数据库中获得，组装 ID (Assembly ID) 分别是 GCA_903112645 (Currot haplotype) 和 GCA_903114435 (OrangeRed haplotype)。

2.7 使用 FALCON-Phase 工具获得两个全长伪单倍型序列

FALCON-Phase 程序需要输入主要重叠群 (primary contigs) 和单倍型特异性重叠群 (haplotigs)。FALCON 从 PacBio CLR 长读段组装得到主要重叠群 (p_contigs) 和次级重叠群 (a_contigs)；接着 FALCON-Unzip 算法利用主要重叠群中的杂合信息进行单倍型定相，最终生成一组更新的主要重叠群和单倍型特异性重叠群 (Chin et al 2016)（可利用长读段测序数据进行自身校正）。单倍型特异性重叠群比 FALCON 生成的次级重叠群体现更连续的单倍型特异性序列信息。最终得到 237,925,913 bp (N50: 364,944 bp) 的主要重叠群和 117,123,257 bp (N50: 74,370 bp) 的单倍型特异性重叠群，并将两者输入 FALCON-Phase 程序，FALCON-Phase 识别 Hi-C 测序数据中的超远程链接信息对基因组进行进一步分型，最终得到两个重叠群 (contigs) 水平的全长伪单倍型 (pseudo-haplotypes)。

2.8 组装错误纠正

Pilon 是一款集成基因组组装改进和变异检测的软件，通过将基因组组装与原始 Illumina 测序数据进行比对，可以检测并校正单碱基差异 (single base differences)、插入缺失 (indels)、局部组装错误 (local misassemblies)、填充间隙 (gap filling) 等 (Walker et al 2014)。利用 bowtie2 v2.4.5 以默认参数构建每个单倍型组装与子代 10x Genomics 测序数据的对齐 Bam 文件，并通过 samtools v1.13 进行排序，以及构建索引，最后以参数 “–fix bases –mindepth 0.85” 运行 pilon v1.24。

2.9 组装质量评估 - 组装连续性

N50 是评估基因组连续性的重要指标。首先需要将组装中所有序列按照长度从长到短依次排列，并按照该顺序将序列长度依次相加，当相加的长度达到序列总长度的 50% 时，最后一条序列的长度便是 N50。N90 是当相加的长度达到序列总长度的 90% 时，最后一条序列的长度。N50，N90 值越大，则表明组装连续性越高。

2.10 组装质量评估 - 组装完整性

基于 OrthoDB v10 的有胚植物单拷贝同源物数据库 (embryophyta_odb10)，对每个单倍型基因组进行 BUSCO v3.0.2 评估，可以显示基因组的完整率 (C)、重复率 (D)、缺失率 (M) 等指标。

2.11 单倍型组装质量评估

Merqury 可以在没有参考基因组的情况下，基于 k-mers 频率分布对组装进行完整性和质量评估；更重要的是其可以对单倍型分辨的两个单倍型基因组进行单倍型组装质量评估。首先，利用 Meryl 统计两个亲本（Currot 和 OrangeRed）Illumina 测序数据和子代 10x Genomics 数据的 k-mers 频率分布。

在 Merqury 中，hap-mers 被定义为亲本遗传的单倍型特异性 k-mers，代表着子代从亲本中继承的 k-mers，仅在单个单倍型中出现一次或多次，用集合的思想表示为亲本各自特异的 k-mers 与子代 k-mers 的交集。Hap-mers 用于确定 Merqury 中定义的相位块 (phase block)，相位块被定义为源自相同单倍型的一组一致性的标记 (Rhie et al 2020)。Merqury 最终可以输出拷贝数谱图 (spectra-cn plot)、组装谱图 (spectra-asm plot)、评估单倍型组装质量的圆圈图 (hap-mers blob plot)(Rhie et al 2020)。

杂合二倍体的拷贝数谱主要由两个峰组装，一个代表在二倍体基因组为 1 个拷贝 (1-copy) 的 k-mers（杂合，在单个单倍型上出现）；另一个在基因组中代表 2 个拷贝 (2-copy) 的 k-mers（纯合，在两个单倍型上或在一个单倍型中出现两次）(Rhie et al 2020)。最终生成一组多重性直方图，并以拷贝数着色。

组装谱与拷贝数谱类似，通过 k-mers 在不同的单倍型来着色，主要显示两个单倍型中的 k-mers 分布。可以检测到只存在于一个单倍型组装中的 k-mers、在两个单倍型组装中同时存在的 k-mers (shared)、仅在原始测序数据中发现的 k-mers (Rhie et al 2020)。其次，组装谱可以指示 k-mers 的完整性，从而评估两个单倍型组装的完整性（如果其中出现了黑色峰，则说明单倍型组装不完整，只有部分单倍型被识别，剩下的仅在原始测序数据中被发现）。圆圈图 (blob plot) 可以直观地显示两个单倍型的分型效果，通过单倍型类别着色图中的圆圈，圆圈的位置将指示单倍型基因组组装的好坏（圆圈越靠近坐标轴，定相效果越好）。

3 结果与分析

Gamete binning、Trio binning、FALCON-Phase 都得到了不同程度的单倍型分辨的单倍型基因组序列。Gamete binning 的单倍型组装序列是从 ENA 数据库获得，已被校正；而 Trio binning、FALCON-Phase 得到的两个单倍型基因组均通过 Pilon v1.24 进行校正，纠正错误组装。Pilon 可以从原始测序数据中恢复组装中缺失的基因或修正组装错误，尤其是在 Trio binning 中表现突出。根据单拷贝基因完整性评估 (BUSCO)，Trio binning 得到的两个单倍型在 pilon 校正之前均存在约 6.0% 的缺失，检测到的单拷贝基因占比仅有 87.0%；经 pilon 校正之后，缺失率降低至 2.6% 以下，单拷贝基因占比也达到了 91.7% 以上（表 3），说明 Pilon 从原始测序数据中有效的恢复了由组装错误导致的缺失。ALLHiC 方法并没有得到两个单倍型分辨的基因组，但是我们通过不同的 k 值获得两个单倍型混合基因组（k=8 和 k=16），ALLHiC 提供的二倍体模型 (ALLHiC_pip.sh) 并不足以支撑单倍型解析的单倍型基因组组装；正如我们预期的一样，ALLHiC（k=16，二倍体，染色体条数为 8）得到的基因组并没有两两同源的染色体，更像是 ALLHiC (k=8) 组装基因组的碎片（图 S4）。输入 k=8 管道的主要重叠群是经过 Purge haplotigs 以参数 “-l 7 -m 60 -h 120” 处理（去除冗余的单倍型特异性重叠群，保留主要重叠群）。

本文主要从基因组连续性 (N50)、基因组组装完整性 (BSUCO)、单倍型组装质量 (Merqury) 等方面评估 Gamete binning、Trio binning、FALCON-Phase、ALLHiC 4 种单倍型组装方法对二倍体物种杏 (Prunus armeniaca) 的单倍型基因组组装质量。

表 2 每个单倍型基因组组装统计
Table 2 Genome assembly statistics per haplotype

方法	单倍型	序列 / 条	N50/bp	N90/bp	总碱基数 /bp
Gamete binning	Currot	93	25,835,044	20,833,752	215,952,222
Gamete binning	OrangeRed	104	25,508,227	20,530,804	215,239,966
FALCON-Phase	phased.0	1,321	362,109	75,530	238,169,392
FALCON-Phase	phased.1	1,321	362,109	75,530	238,042,337
Trio binning	F1-haplotypeCurrot	1,648	484,325	49,148	237,732,501
Trio binning	F1-haplotypeOrangeRed	1,664	489,043	46,145	236,509,352
ALLHiC	k=8	8	24,510,357	22,706,253	218,080,681
ALLHiC	k=16	16	25,351,798	15,217,824	367,847,792

3.1 评估基因组组装连续性

N50 是评估基因组组装连续性的重要指标。N50 越大，表示组装连续性越高。此外，对组装的基因组大小、序列条数进行统计（表 2）。每种方法（除 ALLHiC）获得的两个单倍型之间差距不大，说明获得的两个单倍型基因组组装在连续性上有一定的相似性。获得的单倍型基因组大小在 210 Mbp-240 Mbp 之间，符合 GenomeScope2 预测的 250 Mbp（图 1 a）单倍型基因组大小的预期（误差在可接受范围之内）。染色体水平的单倍型基因组 N50 都在 20 Mbp 以上，而 FALCON-Phase 与 Trio binning 获得的单倍型组装 N50 均在 500 kbp 以下。从连续性上来看，效果最佳的是 Gamete binning 方法。如果对 FALCON-Phase 和 Trio binning 的单倍型基因组组装进行 Scaffolding，可以一定程度上改善组装的连续性，但也有可能引入未知的错误。ALLHiC (k=8) 是在去除单倍型特异性重叠群 (haplotigs) 之后运行的结果，整体的基因组大小也在正常范围之内；理想状态下，ALLHiC (k=16) 基因组中含有目标个体的两套遗传信息，即基因组大小应该是单倍型基因组的两倍，但是 k=16 得到的基因组大小只有单倍型基因组的 1.69 倍，说明二倍体基因组约 30% 的等位变异没有被分离。在完全分离等位变异的情况下，ALLHiC (k=16) 可以通过与参考基因组进行点阵图分析分离同源染色体（图 S4），但是本文得到的 k=16 基因组并没有这个特征。

3.2 评估基因组组装完整性

使用 OrthoDB (Zdobnov et al 2021) 的单拷贝直系同源物数据库，BUSCO 可以基于预期基因含量提供具有生物学意义的指标来估计基因组的完整性和冗余性 (Manni et al 2021)。杏属于有胚植物，选择有胚植物直系同源数据库 (embryophyta_odb10)，其中共有 1614 个保守的单拷贝基因。4 种单倍型基因组组装方法得到的基因组完整性 (C) 均在 93.2% 以上，且缺失率 (M) 在 5.9% 以下（表 3），说明各种方法基本恢复了单倍型基因组中的单拷贝基因。FALCON-Phase 与 ALLHiC (k=16) 的单倍型基因组组装存在 9.0% 以上的重复 (D)（表 3）。尽管 ALLHiC 并没有分离二倍体的单倍型基因组序列，但其组装的基因组序列完整性较高；将 k 值设置为二倍体染色体条数 (k=8)，则可以获得 Scaffolding 之后的基因组 Scaffold 水平的组装。可以进一步提升组装序列的连续性。

表 3 单倍型基因组的 BUSCO 评估（经 pilon 校正）
Table 3 BUSCO assessment of haplotype genomes (pilon-corrected)

方法	单倍型	完整率 (C)/%	单拷贝率 (S)/%	重复率 (D)/%	缺失率 (M)/%
Gamete binning	Currot	97.9	96.1	1.8	1.4
Gamete binning	OrangeRed	97.6	95.3	2.3	1.3
FALCON-Phase	phased.0	96.7	87.6	9.1	2.2
FALCON-Phase	phased.1	96.7	87.7	9.0	2.2
Trio binning	F1-haplotypeCurrot	97.1	92.6	4.5	1.8
Trio binning	F1-haplotypeOrangeRed	96.3	91.7	4.6	2.6
ALLHiC	k=8	97.9	94.7	3.2	1.1
ALLHiC	k=16	93.2	82.2	11.0	5.9

3.3 评估单倍型基因组组装的准确性

利用两个亲本 Illumina 全基因组测序数据来验证单倍型基因组的准确性。Merqury 提出基于 hap-mers 评估单倍型的准确性，相较于只考虑亲本特异性 k-mers 更加严格。亲本遗传的单倍型特异性 k-mers (hap-mers) 是子代从亲本中继承而来的 k-mers，仅在一个单倍型中出现一次或多次，可作为评估单倍型准确性的指标。经统计，来自亲本 Currot 和 OrangeRed 的亲本遗传的单倍型特异性 k-mer 数量分别为 15,033,008 条和 19,212,772 条。

Gamete binning 方法得到的两个单倍型基因组均在相应亲本基因组中找到 92.70% 以上的 hap-mers。Trio binning 方法生成的两个单倍型更是在 95.85% 以上（表 4），表明大部分的单倍型序列都被正确分配给了单倍型。FALCON-Phase 得到的两个单倍型的 hap-mers 在两个亲本中的占比相似，说明其单倍型为单倍型融合基因组，并没有将基因组中的单倍型序列分离开。ALLHiC (k=8) 与 FALCON-Phase 类似，但其中只包含 49.84% 左右的 hap-mers，说明 ALLHiC (k=8) 的单倍型混合基因组并不完整；而 ALLHiC (k=16) 虽然恢复了更多的单倍型信息，但仍是混合的单倍型基因组。

表 4 各个单倍型中来自不同亲本 hap-mers 的统计
Table 4 Statistics of hap-mers from different parents in each haplotype

方法	单倍型	Currot.hap-mer	占比 /%	Orange.hap-mer	占比 /%
Gamete binning	Currot	14,006,321	93.17	679,810	3.54
Gamete binning	OrangeRed	531,039	3.53	17,809,062	92.70
FALCON-Phase	phased.0	8,575,279	57.04	10,805,973	56.24
FALCON-Phase	phased.1	8,552,823	56.89	10,812,787	56.28
Trio binning	F1-haplotypeCurrot	14,408,428	95.85	81,132	0.42
Trio binning	F1-haplotypeOrangeRed	59,666	0.40	18,495,382	96.27
ALLHiC	k=8	7,492,212	49.84	9,535,258	49.63
ALLHiC	k=16	13,590,882	90.41	17,478,066	90.97

Hap-mers 圆圈图 (hap-mers blob plot) 可以直观的显示各条组装序列的单倍型来源，从而评估单倍型分型质量。从两个亲本中统计各自特异的 hap-mers（表 4），作为圆圈图的 x 轴和 y 轴，显示理想状态下的单倍型分布；圆圈颜色对应着不同的单倍型，大小为序列长度。理想状态下的圆圈图是两种颜色的圆圈（两个单倍型）无限接近于两轴（如 Trio binning 的圆圈图）。Gamete binning 与 Trio binning 的分型效果最佳（图 2），几乎将所有的单倍型序列正确分配给了每个单倍型。FALCON-Phase 则是两个单倍型序列的混合（紫色圆圈，红色与蓝色的重合）。

图 2 hap-mers 圆圈图
Fig.2 Hap-mer blob plots
两轴分别是来自两个亲本 Currot 和 OrangeRed 原始测序数据的 hap-mers 计数；圆圈颜色表示不同的单倍型，大小表示序列长度。理想状态下，两种不同颜色的圆圈分别无限接近于两轴，不重叠。a) Gamete binning 方法两个单倍型的 hap-mers 圆圈图。b) FALCON-Phase 方法两个单倍型的 hap-mers 圆圈图。c) Trio binning 方法两个单倍型的 hap-mers 圆圈图。

3.4 评估单倍型基因组冗余性和杂合性

BUSCO 通过检测组装中是否存在给定数据库中的单拷贝基因来评估组装的完整性。Merqury 通过组装基因组和原始测序数据两者之间的 k-mers 计数进行统计比较，绘制 k-mer 拷贝数谱 (k-mer copy number spectra)，从而指示拷贝数分布，可以在一定程度上显示组装基因组的冗余性和杂合性，以及缺失率 (Mapleson et al 2017)。k-mer 拷贝数谱，即 k-mer 拷贝数直方图。一般情况下，二倍体基因组的 k-mer 拷贝数谱会出现两个峰，根据 k-mer 在每个单倍型中出现的次数来决定拷贝数。第一个峰为单拷贝峰 (1-copy)，表示只在其中一个单倍型中发现的 k-mers，也被称为杂合峰；第二个峰为双拷贝峰 (2-copy)，表示同时在两个单倍型中找到，也就是纯合峰。

由于 ALLHiC 方法并没有得到两个单倍型解析的基因组，故本节只对 Gamete binning、FALCON-Phase 和 Trio binning 三种方法得到的单倍型解析基因组绘制 k-mer 拷贝数谱（图 3），图中不同的颜色表示不同的拷贝数。FALCON-Phase 中出现了 4 拷贝 (4-copy) 峰（蓝色峰下的紫色峰），说明 FALCON-Phase 方法得到的单倍型中存在一定的错误重复，正如 BUSCO 统计其两个单倍型的重复率均在 9.0% 以上（表 3）。Gamete binning 与 Trio binning 方法运用的是相同的策略，都是在组装前将子代长读段测序数据分成两组不同的读段，再进行分别组装；从图 3 可看出两者的拷贝数差异不大。反观 FALCON-Phase 的两个单倍型存在大量缺失，图 3 显示组装中的单拷贝 k-mers 在原始测序数据中大量出现（黑色峰），说明 FALCON-Phase 只恢复了测序数据中的部分单倍型信息。此外，FALCON-Phase 的两个单倍型中几乎没有单拷贝峰，说明组装得到的两个单倍型序列几乎都是纯合的（k-mer 拷贝数谱中只有 2 拷贝峰），且两个单倍型的序列高度相似，所以我们猜测 FALCON-Phase 组装缺失的可能是单倍型的杂合序列（黑色峰）。K-mer 拷贝数谱还显示各个单倍型组装中均存在不同程度的碱基错误，因为在单倍型组装中发现了在原始测序数据中不存在的 k-mers（在 0 覆盖度下均存在红色竖线）。

图 3 K-mer 拷贝数谱图
Fig. 3 K-mer copy number spectrum plots
x 轴表示 k-mer 出现的次数，y 轴表示出现 x 次数的 k-mers 数量；不同的颜色表示不同的拷贝数，黑色表示在测序数据中统计的 k-mer 数量，在一定程度上显示单倍型组装的缺失。第一个峰（杂合）表示基因组中的 1 拷贝 (1-copy) k-mers，仅在一个单倍型中发现的 k-mers；第二个峰（纯合）代表来自纯合序列或单倍型特异性重复的 2 拷贝 (2-copy) k-mers，即同时在两个单倍型中存在的 k-mers。a) Gamete binning 方法两个单倍型的 k-mer 拷贝数谱。b) FALCON-Phase 方法两个单倍型的 k-mer 拷贝数谱。c) Trio binning 方法两个单倍型的 k-mer 拷贝数谱。

3.5 评估单倍型基因组组装的完整性

k-mer 拷贝数谱图是从所有的拷贝数进行评估两个单倍型基因组的质量，k-mer 组装谱图 (spectra-asm) 则是分别统计两个单倍型组装的拷贝数，并以单倍型类别进行着色（图 4）。一个完美组装的二倍体基因组预计具有特异于代表杂合等位序列的每个单倍型 k-mer 平衡 (Rhie et al 2020)。Gamete binning 与 Trio binning 除了存在可能的碱基错误之外，几乎是完美的；其中 1 拷贝峰 k-mers 对每个单倍型组装都是特异的，2 拷贝峰（绿色峰）则是两个单倍型组装共有的 k-mers。这也验证了前面提到的高于 92.70% 的 k-mer 完整性（表 4）。FALCON-Phase 单倍型组装中并没有发现特异于单个单倍型的 k-mers，所有统计的 k-mers 在两个单倍型中共享，且在 1 拷贝和 2 拷贝中都存在不同程度的缺失（黑色峰）。

图 4 K-mer 组装谱图
Fig. 4 K-mer assembly spectrum plots
从左至右，第一个峰是纯合峰，以单倍型着色；紫色（红色与蓝色重叠）峰是两个单倍型中特异性的 k-mers，绿色则代表在两个单倍型中同时存在的 k-mers。a) Gamete binning 方法两个单倍型的 k-mer 组装谱。b) FALCON-Phase 方法两个单倍型的 k-mer 组装谱。c) Trio binning 方法两个单倍型的 k-mer 组装谱。

3.6 统计单倍型基因组的交换错误

交换错误 (switch error) 是对等位变异定相到错误单倍型的统计，是评估基因分型质量的重要指标。在单倍型中检测到的来自其他单倍型的 k-mers，被定义为交换错误。Merqury 将存在交换错误的重叠群打断成相位块 (phase block)，相位块被定义为至少具有两个来自相同单倍型特异性 k-mers (hap-mers) 的连续序列 (Rhie et al 2020)。这里设置的阈值是每 20 kbp 允许存在 100 次交换，即允许 0.5% 的交换错误发生。此外，还对相位块进行统计 N50、序列总长度等信息，方便定位交换错误在基因组中发生的位置。单倍型解析基因组组装的整个单倍型可以算做一个单独的块 (block)(Rhie et al 2020)；理想状态下，在不存在交换错误时，统计的相位块 N50 与组装序列 N50 相同，并且相位块的总碱基数将与组装序列总碱基数持平。

4 种方法得到的单倍型组装中只有 Gamete binning、Trio binning 和 FALCON-Phase 得到的是单倍型解析的基因组组装，三者相位块的总长度与其组装序列总长度相近，但是相位块的 N50 相差较大。经统计，Gamete binning 的单倍型解析基因组交换错误率最低，在 0.3515% 以下。Trio binning 也控制在 0.5000% 以下，但是 FALCON-Phase 的交换错误率超出了预先设置的阈值 (0.5000%)，高达 0.6161%，这也验证了 FALCON-Phase 的单倍型之间都混合了另一个单倍型的等位变异序列。而 ALLHiC 无论是 k=8 还是 k=16，交换错误率均在 1.6491% 以上，这也说明了 ALLHiC 获得的杂合二倍体基因组是单倍型混合的组装。

表 5 交换错误率统计
Table 5 The statistics of switch error rate

方法	单倍型	块数	总长度 /bp	最长 /bp	N50/bp	交换错误率 /%
Gamete binning	Currot	658	206,632,055	4,213,615	886,066	0.3515
Gamete binning	OrangeRed	599	207,177,951	6,596,059	1,188,055	0.2800
FALCON-Phase	phased.0	2,854	215,706,318	1,117,862	150,569	0.6173
FALCON-Phase	phased.1	2,853	215,499,112	1,117,860	150,410	0.6161
Trio binning	F1-haplotypeCurrot	2,292	215,912,840	2,630,238	407,484	0.4995
Trio binning	F1-haplotypeOrangeRed	2,174	216,187,017	2,645,633	409,347	0.3155
ALLHiC	k=8	2,323	210,205,802	2,290,945	194,128	2.1485
ALLHiC	k=16	4,106	356,383,539	2,277,856	158,390	1.6491

3.7 总结

结合以上几种单倍型基因组评估指标（表 6），发现 Gamete binning 与 Trio binning 这两种方法得到的单倍型基因组较为完整，其中 Gamete binning 方法整体表现最好。FALCON-Phase 方法分离得到的二倍体两组单倍型序列准确度不高，其中交换错误率超出了阈值（表 5），说明单倍型序列中存在另一单倍型序列。ALLHiC 方法则没有从基因组初步组装中分离出完整的单倍型基因组序列，但是可以利用 ALLHiC 对基因组组装进行 scaffolding，从而提高基因组组装连续性。

表 6 4 种单倍型基因组组装方法的优缺点
Table 6 Advantages and disadvantages of 4 haplotype genome assembly methods

评估指标	最好	较好	较差	最差
组装连续性	Gamete binning	ALLHiC	Trio binning	FALCON-Phase
组装完整性	Gamete binning	Trio binning	ALLHiC	FALCON-Phase
单倍型准确性	Trio binning	Gamete binning	FALCON-Phase	ALLHiC
交换错误	Gamete binning	Trio binning	FALCON-Phase	ALLHiC

4 讨论

本文最终得到了 Gamete binning、Trio binning 和 FALCON-Phase 3 种单倍型基因组组装方法的单倍型解析基因组，分别从组装序列连续性、完整性、分型质量等方面整体评估单倍型组装质量。利用 Trio binning 方法在组装前将子代长读段测序数据以单倍型进行分组，再分别进行独立组装，在很大程度上避免了等位变异对单倍型组装的影响。Gamete binning 借鉴了 Trio binning 方法的分型策略，只是前期单倍型信息来源不同。Trio binning 利用亲本各自特异的 k-mers 将子代长读段测序数据进行分组（每组作为单独的单倍型），依赖于目标个体的两个亲本 Illumina 短读段测序数据，但这对于部分难以获得亲本的样本是行不通的；Gamete binning 则是从上百个配子 10x Genomics 高通量单细胞测序数据中获取用于分型的单倍型信息，在没有亲本测序数据的情况下完成对子代 PacBio 长读段测序数据的每个连锁群进行分组（每个连锁群代表一条染色体；每组表示一个染色单体信息），再对每个连锁群进行独立的单倍型组装，并基于遗传图谱将单倍型基因组挂载到染色体水平（Trio binning 只能得到重叠群水平的组装，染色体水平需要通过其他 scaffolding 工具获得）。FALCON-Phase 更多的是算法上的突破，利用 Hi-C 测序数据中特有的定相信息对基因组进行定相，在一定程度上可以降低定相的计算复杂度 (Kronenberg et al 2021)。ALLHiC 同样从 Hi-C 测序数据中获取用于长读段组装序列分型的信息，但这严重依赖于基因组的杂合性与提供的初步重叠群组装的准确性；杏 (Prunus armeniaca) 的杂合性为 1.06%，还不足以 ALLHiC 识别长读段组装序列中所有的等位变异，根据 ALLHiC 对二倍体物种茶 (Camellia sinensis) 的应用，杂合度在 2.00% 以上可能会有更好的分型效果 (Zhang et al 2021)。ALLHiC 用于将组装序列挂载到染色体水平会是一个不错的选择。

Gamete binning 和 Trio binning 的分型质量最佳（表 6），但是 Trio binning 的组装序列连续性较差，不过这可以通过 scaffolding 工具来提升 N50 的长度。FALCON-Phase 和 ALLHiC 方法得到的单倍型组装都混合了另一单倍型序列，我们预测 FALCON-Phase 的分型质量与前期提供的主要重叠群和交替的单倍体重叠群质量有关；经统计，其中只有约 49% 的基因组被解压缩 (unzip)，不足以支撑单倍型的分离，这也与基因组的杂合性有关（基因组的解压缩率与杂合性成正比），且 FALCON-Phase 组装序列中出现大于 9.0% 的重复序列，这可以考虑利用 Purge haplotigs (Roach et al 2018) 程序去除冗余。ALLHiC 目前只从长读段初级组装中获得二倍体单倍型混合的基因组，并没有得到预期的单倍型解析基因组，故没有将其同另 3 种组装方法一同比较，而是对不同的 k 值进行比较（详见附录 A）。ALLHiC 的组装质量严重依赖于输入的长读段初级组装，以及合理的 k 值设置。

本文几种单倍型组装方法比较并不全面，用于单倍型解析基因组组装的测序技术还有 PacBio HiFi 高保真测序，DipAsm 与 hifiasm 均可以通过 HiFi 测序数据获得单倍型解析基因组；由于没有获得公开的杏 (Prunus armeniaca) PacBio HiFi 测序数据，故没有将其加入一同比较，这也是本文的遗憾与不足。PacBio HiFi 测序数据比 PacBio CLR 测序数据具有更准确的优势，可以生成连续性更好、准确度更高的基因组组装。同时，PacBio HiFi 测序数据由于高准确度，可以在较低测序深度（如 30x-50x）即可获得连续性较好的基因组；反观 PacBio CLR 测序需要至少 80x 才能满足基因组组装的需求。

最后，就适用性而言，Trio binning 对新手最友好，集成于 Canu 程序，但耗资源、运行慢，依赖于集群作业。Gamete binning 是诸多相关程序的整合，涉及的依赖也相应增加，并且需要上百个配子单细胞测序用于分型，但 Gamete binning 可确定单倍型序列来自于哪个亲本。FALCON-Phase 组装序列的连续性，与输入的主要重叠群和单倍型特异性重叠群密切相关，可利用 FALCON + FALCON-Unzip 或 Canu + Purge-haplotigs 对 PacBio CLR 原始测序数据组装获得。ALLHiC 对于简单二倍体而言，可获得染色体水平单倍型混合的基因组，或许可以通过 ALLHiC 的多倍体模型获得二倍体的单倍型解析基因组，但这也依赖于高质量染色体水平参考基因组生成的等位基因重叠群表 (allelic contig table)。

[1] Circular Consensus Sequencing (CCS) 测序模式产生的测序序列，表示 DNA 分子的循环连续测序序列。

[2] Continuous Long Read 的简称，表示 DNA 分子的连续长测序序列。

[3] 根据 GenomeScope2 预测的单倍型基因组大小 (250 Mbp) 计算，即总碱基数除以基因组大小。

[4] https://github.com/tangerzhang/ALLHiC/blob/master/bin/ALLHiC_pip.sh

[5] 各个单倍型统计的 hap-mers 数量与亲本原始测序数据的 hap-mers 数量比值。

参考文献

[1] 李双双，邓传良，张迎新，胡赞民，范成明，陈宇红。单倍型分析技术研究进展. 生物工程学报 , 2018, 34: 852–861

[2] Al Bkhetan Z, Chana G, Ramamohanarao K, Verspoor K, Goudey B. Evaluation of consensus strategies for haplotype phasing. Brief Bioinform, 2021, 22: 1–12

[3] Cheng H, Concepcion GT, Feng X, Zhang H, Li H. Haplotype-resolved de novo assembly using phased assembly graphs with hifiasm. Nat Methods, 2021, 18: 170

[4] Chin CS, Peluso P, Sedlazeck FJ, Nattestad M, Concepcion GT, Clum A, Dunn C, O’Malley R, Figueroa-Balderas R, Morales-Cruz A, Cramer GR, Delledonne M, Luo C, Ecker JR, Cantu D, Rank DR, Schatz MC. Phased diploid genome assembly with single-molecule real-time sequencing. Nat Methods, 2016, 13: 1050–1054

[5] Campoy JA, Sun H, Goel M, Jiao WB, Folz-Donahue K, Wang N, Rubio M, Liu C, Kukat C, Ruiz D, Huettel B, Schneeberger K. Gamete binning: chromosome-level and haplotype-resolved genome assembly enabled by high-throughput single-cell sequencing of gamete genomes. Genome Biol, 2020, 21: 1–20

[6] Cao H., Wu H., Luo R, Huang S, Sun Y, Tong X, Xie Y, Liu B, Yang H, Zheng H, Li J, Li B, Wang Y, Yang F, Sun P, Liu S, Gao P, Huang H, Sun J, Chen D. et al. De novo assembly of a haplotype-resolved human genome. Nat Biotechnol, 2015, 33: 617–622

[7] Chen X, Yang H, Wong WH. Phased Genome Sequencing Through Chromosome Sorting. Methods Mol Biol, 2017, 1551: 171–188

[8] Edger PP. The power of chromosome-scale, haplotype-resolved genomes. Mol Plant, 2022, 15: 393–395

[9] Falconer E, Lansdorp PM. Strand-seq: a unifying tool for studies of chromosome segregation. Semin Cell Dev Biol, 2013, 24: 643–652

[10] Garg S, Fungtammasan A, Carroll A, Chou M, Schmitt A, Zhou X, Mac S, Peluso P, Hatas E, Ghurye J, Maguire J, Mahmoud M, Cheng H, Heller D, Zook JM, Moemke T, Marschall T, Sedlazeck FJ, Aach J, et al. Chromosome-scale, haplotype-resolved assembly of human genomes. Nat Biotechnol, 2020, 39: 309–312

[11] Hu G, Feng J, Xiang X, Wang J, Salojärvi J, Liu C, Wu Z, Zhang J, Liang X, Jiang Z, Liu W, Ou L, Li J, Fan G, Mai Y, Chen C, Zhang X, Zheng J, Zhang Y, Peng H, et al. Two divergent haplotypes from a highly heterozygous lychee genome suggest independent domestication events for early and late-maturing cultivars. Nat Genet, 2022, 54: 73–83

[12] Kronenberg ZN, Rhie A, Koren S, Concepcion GT, Peluso P, Munson KM, Porubsky D, Kuhn K, Mueller KA, Low WY, Hiendleder S, Fedrigo O, Liachko I, Hall RJ, Phillippy AM, Eichler EE, Williams JL, Smith TPL, Jarvis ED, Sullivan ST, et al. Extended haplotype-phasing of long-read de novo genome assemblies using Hi-C. Nat Commun, 2021, 12: 1935

[13] Koren S, Rhie A, Walenz BP, Dilthey AT, Bickhart DM, Kingan SB, Hiendleder S, Williams JL, Smith TPL, Phillippy AM. De novo assembly of haplotype-resolved genomes with trio binning. Nat Biotechnol, 2018, 36: 1174–1182

[14] Koren S, Walenz BP, Berlin K, Miller JR, Bergman NH, Phillippy AM. Canu: Scalable and accurate long-read assembly via adaptive κ-mer weighting and repeat separation. Genome Res, 2017, 27: 722–736

[15] Mapleson D, Accinelli GG, Kettleborough G, Wright J, Clavijo BJ. KAT: a K-mer analysis toolkit to quality control NGS datasets and genome assemblies. Bioinformatics, 2017, 33: 574

[16] Manni M, Berkeley MR, Seppey M, Zdobnov EM. BUSCO: Assessing Genomic Data Quality and Beyond. Curr Protoc, 2021, 1: e323

[17] Marçais G, Kingsford C. A fast, lock-free approach for efficient parallel counting of occurrences of k-mers. Bioinformatics, 2011, 27: 764

[18] Ranallo-Benavidez TR, Jaron KS, Schatz MC. GenomeScope 2.0 and Smudgeplot for reference-free profiling of polyploid genomes. Nat Commun, 2020, 11: 1432

[19] Roach MJ, Schmidt SA, Borneman AR. a: Purge Haplotigs: allelic contig reassignment for third-gen diploid genome assemblies. BMC Bioinformatics, 2018, 19: 1–10

[20] Rhie A, Walenz BP, Koren S, Phillippy AM. Merqury: Reference-free quality, completeness, and phasing assessment for genome assemblies. Genome Biol, 2020, 21: 1–27

[21] Walker BJ, Abeel T, Shea T, Priest M, Abouelliel A, Sakthikumar S, Cuomo CA, Zeng Q, Wortman J, Young SK, Earl AM. Pilon: An Integrated Tool for Comprehensive Microbial Variant Detection and Genome Assembly Improvement. PLoS One, 2014, 9: e112963

[22] Wenger AM, Peluso P, Rowell WJ, Chang PC, Hall RJ, Concepcion GT, Ebler J, Fungtammasan A, Kolesnikov A, Olson ND, Töpfer A, Alonge M, Mahmoud M, Qian Y, Chin CS, Phillippy AM, Schatz MC, Myers G, DePristo MA, Ruan J, et al. Accurate circular consensus long-read sequencing improves variant detection and assembly of a human genome. Nat Biotechnol, 2019, 37: 1155–1162

[23] Yang H, Chen X, Wong WH. Completely phased genome sequencing through chromosome sorting. Proc Natl Acad Sci USA, 2011, 108: 12–17

[24] Zhang X, Chen S, Shi L, Gong D, Zhang S, Zhao Q, Zhan D, Vasseur L, Wang Y, Yu J, Liao Z, Xu X, Qi R, Wang W, Ma Y, Wang P, Ye N, Ma D, Shi Y, Wang, H, et al. Haplotype-resolved genome assembly provides insights into evolutionary history of the tea plant Camellia sinensis. Nat Genet, 2021, 53: 1250

[25] Zdobnov EM, Kuznetsov D, Tegenfeldt F, Manni M, Berkeley M, Kriventseva E V. OrthoDB in 2020: evolutionary and functional annotations of orthologs. Nucleic Acids Res, 2021, 49: D389

[26] Zhang X, Zhang S, Zhao Q, Ming R, Tang H. Assembly of allele-aware, chromosomal-scale autopolyploid genomes based on Hi-C data. Nat Plants, 2019, 5: 833-845

附录 A

A.1 本文涉及软件及其版本号

表 S1 软件及其版本号
Table S1 Softwares and versions

序号	软件	版本号	下载链接
1	GenomeScope	2.0	https://github.com/tbenavi1/genomescope2.0
2	Jellyfish	2.2.10	https://github.com/gmarcais/Jellyfish/releases/tag/v2.2.10
3	Canu	2.2	https://github.com/marbl/canu/releases/tag/v2.2
4	Purge haplotigs	1.1.2	https://bitbucket.org/mroachawri/purge_haplotigs/src/master/
5	Trio binning	2.2	https://github.com/marbl/canu/releases/tag/v2.2
6	ALLHiC	0.9.13	https://github.com/tangerzhang/ALLHiC
7	FALCON	1.4.2	https://github.com/PacificBiosciences/pb-assembly
8	FALCON-Unzip	1.3.1	https://github.com/PacificBiosciences/pb-assembly
9	FALCON-Phase	1.1.0	https://github.com/PacificBiosciences/pb-assembly
10	Pilon	1.24	https://github.com/broadinstitute/pilon/releases/tag/v1.24
11	BUSCO	3.0.2	https://anaconda.org/thiesgehrmann/busco
12	Meryl	1.3	https://github.com/marbl/meryl/releases/tag/v1.3
13	Merqury	1.3	https://github.com/marbl/merqury/releases/tag/v1.3
14	Minimap2	2.24	https://github.com/lh3/minimap2/releases/tag/v2.24
15	Samtools	1.15.1	https://github.com/samtools/samtools/releases/tag/1.15.1
16	Bowtie2	2.4.5	https://anaconda.org/bioconda/bowtie2

A.2 本文涉及的代码与脚本

A.2.1 基于 Kmer 的基因组特征分析

GenomeScope2 v2.0 运行 Shell 脚本：

k=21
h=1000000
name=rojpas
fq1=trimmed_bamtofastq_R1.fastq
fq2=trimmed_bamtofastq_R2.fastq
## Unzip
<path>/pigz-2.6/unpigz -p 8 ${fq1}.gz
<path>/pigz-2.6/unpigz -p 8 ${fq2}.gz
## count and histo
<path>/jellyfish-2.2.10/bin/jellyfish count ./ -C -o ${name}_${k}mer -t 2 -s 30G \
-m ${k} ${fq1} ${fq2}
<path>/jellyfish-2.2.10/bin/jellyfish histo -t 2 \
-o ${name}_${k}mer.histo ${name}_${k}mer -h $h
## Genomescope2 
<path>/genomescope2 -i ${name}_${k}mer.histo -o Result_m1e6 -k $k -p 2 -m $h

A.2.2 二倍体基因组初步组装

Canu v2.2 组装 PacBio CLR 长读段获取主要重叠群（Linux 命令终端运行）：

1
2
3

$ canu -p rojpas -d out_canu_result genomeSize=250m corOutCoverage=200 \
		 "batOptions=-dg 3 -db 3 -dr 1 -ca 500 -cp 50" -pacbio $fasta

Purge haplotigs v1.1.2 去除主要重叠群中的单倍型特异性重叠群：

## step0 – alignment
assembly=  ## Canu获得的主要重叠群
pacbio=   ## 目标个体的PacBio CLR长读段测序fasta文件
minimap2 -ax map-pb ${assembly} ${pacbio} \
    	| samtools view -hF 256 - \
    	| samtools sort -@ 8 -m 1G -o aligned.bam -T tmp.ali
## step1 - hist
purge_haplotigs hist -b aligned.bam  -g ${assembly} -t 8
## step2 - cov step1-Low Cutoff-Mid Cutoff-High Cutoff 
##（查看step1生成的直方图进行判断-l、-m、-h的值）
purge_haplotigs  cov  -i aligned.bam.gencov  \
		-l 7 -m 60 -h 120 -o coverage_stats.csv
## step3 - purge
purge_haplotigs purge -g ${assembly} -c coverage_stats.csv -b aligned.bam

A.2.3 基于 Trio binning 方法获得单倍型组装序列

Canu v2.2 运行 trio binning 组件命令（Linux 命令终端运行）：

$ canu -p F1 -d Result genomeSize=250m \
		-haplotypeCurrot 0_RawData/CurrotParental/*.fastq.gz \
		-haplotypeOrangeRed 0_RawData/OrangeRedParental/*.fastq.gz \
		-pacbio 0_RawData/pacbioF1/*.fasta \
		hapUnknownFraction=0.01 \
		corMhapSensitivity=high corMinCoverage=0 corOutCoverage=200 \
		correctedErrorRate=0.105 useGrid=true

A.2.4 基于 ALLHiC 方法获得染色体水平的单倍型基因组序列

以 “k=16” 运行 ALLHiC 管道（Linux 命令终端运行）：

1
2
3

$ sh ALLHiC_pip.sh -r assembly.fasta \
		-1 hic_R1.fastq.gz -2 hic_R2.fastq.gz -k 16 -e GATC -t 16 -b 500k

以 “k=8” 运行 ALLHiC 管道（Linux 命令终端运行）：

1
2
3

$ sh ALLHiC_pip.sh -r assembly_purge.fasta \
		-1 hic_R1.fastq.gz -2 hic_R2.fastq.gz -k 8 -e GATC -t 16 -b 500k

A.2.5 基于 Gamete binning 方法获得染色体水平的单倍型组装序列

由 Gamete binning 方法获得杏 (Prunus armeniaca) 的两个单倍型基因组序列下载链接分别是 http://www.ebi.ac.uk/ena/data/view/GCA_903112645 (Currot) 和 http://www.ebi.ac.uk/ena/data/view/GCA_903114435 (OrangeRed)

A.2.6 使用 FALCON-Phase 工具获得两个全长伪单倍型序列

FALCON 参数 (fc_run.cfg)：

input_fofn=input.fofn
input_type=raw
pa_fasta_filter_option=pass
target=assembly
skip_checks=False
LA4Falcon_preload=false
pa_DBsplit_option=-x500 -s200
ovlp_DBsplit_option=-x500 -s200
pa_HPCTANmask_option=
pa_REPmask_code=0,300;0,300;0,300
genome_size = 244847837
seed_coverage = 30
length_cutoff = -1  
pa_HPCdaligner_option=-v -B128 -M24
pa_daligner_option= -k19 -e0.80 -l1000 -h256 -w8 -s100
falcon_sense_option=--output-multi --min-idt 0.70 --min-cov 4 --max-n-read 200
falcon_sense_greedy=False
ovlp_HPCdaligner_option=-v -B128 -M24 
ovlp_daligner_option= -k24 -e.92 -l1800 -h1024 -s100
length_cutoff_pr=1000
overlap_filtering_setting=--max-diff 100 --max-cov 100 --min-cov 2

FALCON-Unzip 参数 (fc_unzip.cfg)（利用 PacBio 长读段进行校正）：

max_n_open_files = 100
input_fofn=input.fofn
input_bam_fofn=input_bam.fofn
polish_include_zmw_all_subreads = true

FALCON-Phase 参数 (fc_phase.cfg)：

min_aln_len=100000
iterations=10000000
enzyme=GATC
output_format=pseudohap

A.2.7 组装错误纠正

Pilon 纠正错误组装（Shell 脚本）：

## index
bowtie2-build -f $asm $asm --threads 4
## Aligned
bowtie2 -x $asm -1 $R1 -2 $R2 -p 4 \
		| samtools view -@ 20 -bS – 
		| samtools sort -@ 20 -o ${name}.aligned.sorted.bam -
## samtools index
samtools index -@ 4 ${name}.aligned.sorted.bam
## polish
java -Xmx200G -jar <path>/pilon-1.24-0/pilon.jar \
		--genome $asm --bam ${name}.aligned.sorted.bam \
		--fix bases --mindepth 0.85 --output ${name}.polished –vcf

A.2.8 组装质量评估 - 组装完整性

BUSCO v3.0.2 评估每个单倍型的组装完整性（Linux 命令终端运行）：

1
2
3

$ run_busco -i ${name}.fasta -l <path>/database_BUSCO/embryophyta_odb10 \
		-o ${name}-busco -m genome --cpu 8

A.2.9 单倍型组装质量评估

运行 Meryl 和 Merqury（Shell 脚本）：

## 统计illumina数据中的k-mers频率分布，对亲本与子代都运行以下程序以生成各自的meryl数据库
## meryl 数据库均以k=21构建（即k-mers集合）
read1=R1
read2=R2
name=ChildF1 ## Currot or OrangeRed
meryl k=21 count output ${read1}.meryl ${read1}.fastq.gz
meryl k=21 count output ${read2}.meryl ${read2}.fastq.gz
meryl union-sum output ${name}.meryl ${read1}.meryl ${read2}.meryl
## 获取hap-mers
## 输入亲本与子代的meryl数据库，自动生成单倍型特异性k-mers（hap-mers）
## Currot.meryl和OrangeRed.meryl分别是两个亲本k-mers，ChildF1.meryl是子代原始测序数据统计得来的k-mers
sh hapmers.sh Currot.meryl OrangeRed.meryl ChildF1.meryl
## 运行Merqury
sh merqury.sh ChildF1.meryl Currot.hapmer.meryl OrangeRed.hapmer.meryl \ 
		hap1.fasta hap2.fasta $asm_name

A.3 QV

较高的 QV 表示更准确的一致性 (consensus)，Q30 对应于 99.9% 的准确度，Q40 则表示准确度高达 99.99%。

表 S2 一致性质量
Table S2 QV

方法	单倍型	QV	错误率
Gamete binning	Currot	31.9748	0.000635
Gamete binning	OrangeRed	31.7811	0.000664
FALCON-Phase	phased.0	34.7795	0.000333
FALCON-Phase	phased.1	34.7788	0.000333
Trio binning	F1-haplotypeCurrot	32.2219	0.000600
Trio binning	F1-haplotypeOrangeRed	32.0279	0.000627
ALLHiC	k=8	35.0411	0.000313
ALLHiC	k=16	32.5396	0.000557

A.4 Spectra-cn plots(ALLHiC)

图 S1 ALLHiC 拷贝数谱图
Fig. S1 ALLHiC spectra-cn plots

A.5 Spectra-asm plots(ALLHiC)

图 S2 ALLHiC 组装谱图
Fig. S2 ALLHiC spectra-asm plots
通过上面的拷贝数谱和组装谱图，可以看出 k=8 的组装缺失了一半的单拷贝 k-mers，并在 reads 中找到（在 1-copy 峰处红色与黑色几乎重合）；而 k=16 的光谱图显示其为两个单倍型的混合物，并存在不同程度的缺失，其中只组装了部分单倍型，在 2-copy 中出现了 1-copy 的峰，显示混合单倍型组装部分折叠 (partially collapsed) 了单倍型。

A.6 Hap-mers blob plot(ALLHiC)

图 S3 ALLHiC hap-mers 圆圈图
Fig. S3 ALLHiC hap-mers blob plots
从 hap-mers 圆圈图可以很明显的看出无论是 k=8 还是 k=16，得到的基因组并没有很好的定相，下方两个图正是将其组装序列根据交换错误 (switch error) 打断成相位块 (phase block) 之后重新定相的结果，这也说明输入的基因组中含有两个单倍型信息，只是并没有按照染色体进行区分，而是糅合在一起。

A.7 ALLHiC (k=16) 与 ALLHiC (k=8) 的点阵图

图 S4 点阵图
Fig. S4 Dot plot
横轴为 ALLHiC (k=16) 结果，纵轴为 ALLHiC (k=8) 结果