从测序读段开始组装基因组

先分清输入文件与分析目的#

FASTA 保存序列名称和碱基；FASTQ 还保存逐碱基质量。GFA 是图格式，除了序列片段，还能记录它们的连接。读段的 FASTQ 与组装结果的 FASTA 不能互相替代。

HiFi 和 ONT 都属于长读长测序技术。HiFi 读段的准确性有利于区分相近序列，ONT 超长读段能提供跨越较长重复区域的连续信息；Hi-C 则提供 DNA 片段之间的空间接触信息。联合组装把这些证据用于不同问题，不能把三类文件混在一起当作同一种输入。ONT 测序原理、hifiasm 联合组装方法

本篇命令默认软件已经安装，并且从项目工作目录执行。先为输出建立目录。

1
mkdir -p work/survey work/assembly results logs
2
jellyfish --version
3
genomescope.R --version
4
hifiasm --version
5
seqkit version

软件未找到时，先按官方安装说明配置环境。后面使用的 data/short_R1.fastq、data/short_R2.fastq 是经过质量检查的普通短读段；data/hifi.fastq.gz 是 HiFi 读段。压缩的短读段可先用 gzip -dc 解压到工作目录，原始文件保留。

data/ont_ul.fastq.gz 表示同一个体的 ONT 超长读段，已完成碱基识别和质量检查。不要将不同个体的数据混入一套二倍体组装，也不要仅凭文件名含有 ONT 就认定它符合超长读段输入要求。

联合组装前，分别检查两类长读段#

SeqKit 可以直接读取压缩 FASTQ，汇总各文件的序列条数、总碱基数和读长统计。

1
seqkit stats -a -T \
2
  data/hifi.fastq.gz data/ont_ul.fastq.gz \
3
  > results/long-read-stats.tsv

检查两类数据各自的读长分布与质量报告，尤其要看 ONT 中能跨越目标重复区域的长读段，而不只看最长的一条。SeqKit 的 N50 按碱基总量累计，计算方式与下一篇组装 N50 相同，但这里统计的是原始读段。SeqKit 统计说明

用同一个单倍体基因组大小估计值 G，可以分别报告名义覆盖深度：

B 是对应数据经过质量处理后保留的总碱基数。这个平均值不表示每个区域都被均匀覆盖，也不表示两套同源拷贝各自都获得同样深度。两种技术的读长与误差特点不同，应分别报告，不能只给出相加后的一个覆盖数。

k-mer 频数在数什么#

把一条序列按固定长度切出重叠片段，每次向后移动一个碱基，得到的片段叫作 k-mer。长度至少为 k 的读段，理论上可以取出这些片段：

L 是读段长度，k 是片段长度。含不确定碱基的窗口可能被计数软件跳过；相同片段出现多次，累计到同一个计数中。

频数直方图的横轴是“一个 k-mer 出现了几次”，纵轴是“有多少种 k-mer 出现了这么多次”。低频区域可能混有测序错误，杂合位点与重复序列也会改变曲线形状。只看最高的峰，无法可靠地认定哪个峰对应单拷贝序列。

在理想的单倍体、均匀覆盖且错误已经排除的情况下，可以用一个近似理解基因组大小估计：

分子是可信 k-mer 的总出现次数，分母是单拷贝序列的 k-mer 覆盖深度。二倍体中两个同源拷贝的贡献会改变峰的位置，不能把这个近似公式直接套到任意一个峰上。

用 Jellyfish 计数#

下面把两端普通短读段放在同一次计数中。21 和 100M 展示参数写法，需结合序列长度、复杂度和可用内存调整。

1
jellyfish count \
2
  -m 21 -s 100M -t 8 -C \
3
  -o work/survey/reads.k21.jf \
4
  data/short_R1.fastq data/short_R2.fastq
5

6
jellyfish histo -t 8 work/survey/reads.k21.jf \
7
  > work/survey/reads.k21.histo

-m 设置 k；-t 设置线程数；-C 把一个片段和它的反向互补序列合并计数。-s 100M 设置哈希表容量，表示约一亿个条目，不表示使用 100 MB 内存。实际内存还受 k 和计数存储方式影响。Jellyfish 参数说明

输出 .jf 是计数数据库，.histo 是用于拟合的频数表。检查直方图的高频端是否被输出范围截断；高覆盖或重复较多时，应查看 jellyfish histo --help 中的上限参数，保留分析需要的频数范围。

用 GenomeScope 2 拟合频数曲线#

GenomeScope 2 根据倍性建立 k-mer 频数的混合模型。模型中的不同成分对应不同拷贝数，拟合后估计基因组大小、杂合程度及重复组成。它比直接用总数除峰值考虑了更多情况，但结果仍依赖测序质量、覆盖情况和倍性假设。GenomeScope 2 说明

每个成分使用负二项分布，允许计数的波动大于简单泊松分布所能描述的范围。软件通过非线性最小二乘拟合，让模型曲线与观察到的频数曲线尽量接近。拟合时需要给定倍性，软件不会仅靠这个参数验证物种的真实倍性。GenomeScope 2 方法论文

1
genomescope.R \
2
  -i work/survey/reads.k21.histo \
3
  -o results/genomescope \
4
  -k 21 -p 2

-k 必须与计数时一致，-p 2 表示按二倍体模型拟合。查看输出中的参数报告和拟合图，确认模型曲线是否解释了主要峰形。如果拟合明显偏离观察值，应检查污染、覆盖不均与倍性设定，暂时不要用拟合结果决定后续组装参数。

两种基因组大小估计提供不同证据

流式细胞术通过 DNA 含量估计基因组大小，k-mer 方法从测序序列推断大小。两者都有实验或模型条件。结果不一致时，需要检查各自的误差来源，不能预先认定测序估计一定更准确。

根据数据选择 hifiasm 的组装模式#

hifiasm 在读段重叠关系中保留单倍型差异，再生成连续序列。单倍型可以理解为同一条染色体拷贝上的一组序列。二倍体有两套同源拷贝，组装时要尽量区分它们，并避免把真实重复序列误当成同源拷贝。hifiasm 项目

本篇命令按 hifiasm 0.25.0 核对。先看手上有哪些数据，再选择下列一种模式；使用不同输出前缀，避免混用中间缓存。

1
hifiasm -t 8 -o work/assembly/hifi \
2
  data/hifi.fastq.gz \
3
  2> logs/hifiasm-hifi.log

1
hifiasm -t 8 -o work/assembly/hifi_ont \
2
  --ul data/ont_ul.fastq.gz \
3
  data/hifi.fastq.gz \
4
  2> logs/hifiasm-hifi-ont.log

1
hifiasm -t 8 -o work/assembly/hifi_ont_hic \
2
  --ul data/ont_ul.fastq.gz \
3
  --h1 data/hic_R1.fastq.gz \
4
  --h2 data/hic_R2.fastq.gz \
5
  data/hifi.fastq.gz \
6
  2> logs/hifiasm-hifi-ont-hic.log

1
hifiasm -t 8 -o work/assembly/hic \
2
  --h1 data/hic_R1.fastq.gz \
3
  --h2 data/hic_R2.fastq.gz \
4
  data/hifi.fastq.gz \
5
  2> logs/hifiasm-hic.log

选择正确的输出，再转换为 FASTA#

unitig 是根据明确的读段重叠关系连接出的序列单元，常作为构建 contig 的中间结果。查看文件名时，需要同时确认它属于哪个组装阶段和哪一种运行模式。

命名前面的 hifi_ont_hic 来自自定义输出前缀，后面的 .hic. 来自运行模式。具体产物应与运行版本的输出文档核对。hap1 和 hap2 也不能仅凭名字认定为完整的父源和母源染色体，跨 contig 的来源一致性需要额外证据。

从 HiFi＋ONT 超长读段＋Hi-C 模式的两份 contig 图导出序列，作为下一篇的输入：

1
awk '$1 == "S" {print ">" $2; print $3}' \
2
  work/assembly/hifi_ont_hic.hic.hap1.p_ctg.gfa > results/hap1.contigs.fa
3

4
awk '$1 == "S" {print ">" $2; print $3}' \
5
  work/assembly/hifi_ont_hic.hic.hap2.p_ctg.gfa > results/hap2.contigs.fa

这段命令提取 GFA 的序列行，不能用于不含序列的 noseq.gfa。基础模式应改为自己实际生成的图文件。hifiasm FASTA 导出说明

Verkko 是另一条 HiFi＋ONT 联合路线#

Verkko 同样可以联合高准确度长读段与 ONT 超长读段。它先构建多重 de Bruijn 图，也就是用序列片段的连接关系组织组装，再将 ONT 读段比对到图上，逐步解析重复与分支，最后生成一致序列。Verkko 方法论文

按 Verkko v2.3.2 的接口，HiFi 放在 --hifi 后，ONT 超长读段放在 --nano 后。先查看任务计划，并按计算节点配置资源：

1
verkko -d work/verkko \
2
  --hifi data/hifi.fastq.gz \
3
  --nano data/ont_ul.fastq.gz \
4
  --snakeopts "--dry-run"

确认依赖和资源设置后，去掉 --snakeopts "--dry-run" 才会执行分析。有 Hi-C 时，可按该版本接口补充 --hic1 和 --hic2；这是另一套软件流程，不要求在 hifiasm 后再跑一遍。

Verkko 的主输出包括 assembly.fasta。加入合适分相信息后，还可能产生 assembly.haplotype1.fasta 和 assembly.haplotype2.fasta。应按实际产物接入质检，不套用 hifiasm 的 GFA 文件名。Verkko 输入与输出

完整染色体需要逐条验证

端粒是染色体末端的重复序列；T2T 表示从一端端粒到另一端端粒的完整组装。HiFi＋ONT 可以支持朝这个目标推进，但采用联合组装不等于已经得到 T2T。仍需检查缺口、错误连接与分相质量，也不能把人类端粒重复模式直接套到其他物种。

换组装软件前，先记录问题#

原来的笔记记录过：约 50 倍 HiFi 数据使用 hifiasm 后，Hi-C 互作图中出现了可疑结构，随后尝试 Canu。这次排查的分相状态和 scaffold 方向仍需结合相应接触图与验证记录判断。

Canu 是另一条长读段组装路线。旧文“Canu 在 2021 年停止维护”的说法不成立，官方至少在 2024 年发布过 v2.3。选软件时应核对读段类型、版本和输出目的，再用同一套质量指标比较结果。Canu v2.3 发布记录

下一篇从导出的 FASTA 开始，检查序列连续性与完整性，再说明 Hi-C 挂载和人工校正。