---

桂花双单倍型基因组组装与注释实战全记录

2026年4月12日 组装4-5小时（两个单倍型）

2026年4月12日 挂载10-12小时（两个单倍型）

基因组的组装与注释是一项浩大的系统工程，也是后续群体遗传学、进化分析以及多组学研究的基石。在本项目中，我们针对桂花（Osmanthus fragrans）开展了基因组研究。得益于三代测序技术和深度学习的发展，本次流程采用了目前最前沿的 HiFi + Hi-C 组装策略，以及基于神经网络的基因预测方法。以下是本项目的核心流程全记录。

一、基因组组装

基因组组装是注释流程的基础，其质量直接影响到下游分析的准确性。相比于传统的混合组装，本次我们追求的是极高精度的双单倍型（Diploid/Phased）组装，主要包括数据质控、Contig组装、染色体挂载及手工校正等步骤。

• 数据准备与质量控制
  • 目标：确保原始测序数据的高质量，去除测序接头和低质量嵌合体。
  • 数据类型：基于 PacBio 平台的超长高准确度 HiFi reads，以及提供三维空间互作信息的 Hi-C 数据。
  • 常用工具：利用 FastQC 进行基础质量评估，使用 fastp 对 Hi-C 短读长数据进行质控和过滤。
• 基因组初组装 (Contig Assembly)
  • 目标：在没有参考基因组的情况下，利用长读长数据组装出连续的 Contig 序列，并区分两套单倍体。
  • 方法：利用最新的图算法工具，将 Hi-C 数据与 HiFi 数据联合输入，直接在组装图（Assembly Graph）阶段解开单倍型纠缠。
  • 常用工具：hifiasm。通过其集成的 Hi-C 模式，我们成功输出了两套高纯度的初组装序列（Hap1 与 Hap2），有效避免了传统工具容易产生的杂合序列嵌合问题。
• 染色体挂载 (Scaffolding)
  • 目标：利用 Hi-C 数据的互作频率，将零散的 Contig 锚定、排序并定向到具体的染色体水平。
  • 方法：基于马尔可夫聚类等算法计算 Contig 间的物理距离，构建假染色体（Pseudomolecules）。
  • 常用工具：HapHiC。它对高杂合和多倍体物种非常友好，能自动化生成优化的 .hic 和 .assembly 文件供下游使用。
• 组装后处理与手工校正 (Manual Curation)
  • 目标：纠正自动化挂载中出现的倒位、错拼，剥离冗余碎片，最终提取出桂花标准的 23 条染色体（n=23）。
  • 方法：可视化 Hi-C 互作矩阵，在对角线上寻找完美的“红方块”，通过剪断（Cut）、翻转（Invert）和移动（Move）进行人工“微雕手术”。最后提取纯净的 Top 23 序列并重命名（Chr01-Chr23）。
  • 常用工具：Juicebox（调图神器）、HiCExplorer（绘制最终发表级互作热图）以及 BUSCO（评估基因组完整度）。

二、基因组注释

基因组注释是对组装好的 23 条染色体（物理图谱）进行解读和标记的过程。植物基因组极其庞大且复杂，注释工作必须严格按照：重复序列注释 -> 基因结构注释 -> 基因功能注释的顺序进行。

• 重复序列注释 (Repeat Annotation)
  • 背景：重复序列（如 LTR 等转座子）广泛存在于植物基因组中（通常占比 50% 以上）。如果不加处理，下游软件会把转座子误认为功能基因。
  • 目标：从头构建高质量的转座子文库，并在全基因组范围内进行软/硬屏蔽（Masking）。
  • 常用工具：EDTA（Extensive de-novo TE Annotator）。作为业界金标准，EDTA 内部整合了 LTR_retriever、RepeatModeler、TIR-Learner 等十余款主流软件。通过一键式运行，它可以输出高置信度的 TE 分类文库，并生成屏蔽了重复序列的 .masked 基因组文件。
• 基因结构注释 (Structural Annotation)
  • 背景：在获得屏蔽后的基因组后，需要对真实的基因区（外显子、内含子、UTR等）进行结构预测。传统流程高度依赖转录组（RNA-seq）和近缘物种同源蛋白的反复迭代比对。
  • 目标：快速、精准地在基因组坐标上标注出每一个功能基因的边界。
  • 方法与工具：本次流程摒弃了繁琐的传统 MAKER/BRAKER3 流程，采用了国产最前沿的深度学习注释神器 ANNEVO。
    • 借助高性能 GPU 算力（如 RTX 系列显卡）和 PyTorch 框架，直接调用预训练好的 Embryophyta（有胚植物）模型。
    • 通过分布式的两步法（阶段一：GPU 神经网络序列预测；阶段二：多线程 CPU 基因结构解码），直接从头（ab initio）预测出具备极高生物学意义的基因结构（输出 .gff 文件）。
• 基因功能注释 (Functional Annotation)
  • 目标：为上一步预测出的几万条基因赋予具体的生物学名字和功能描述。
  • 方法：将提取出的蛋白质序列与国际公共数据库进行比对。
  • 常用工具：
    • eggNOG-mapper：极其高效的直系同源簇注释工具，速度快且信息全面。
    • InterProScan：对蛋白结构域和保守基序进行深度扫描。
    • 结合 KEGG 进行代谢通路分析，结合 GO 进行基因本体（分子功能、细胞组分、生物学过程）的注释归类。

结语： 从几十个 GB 的下机数据，到终端里跑出 SFZ_Hap1_Final_Genome.fa 的那一刻；从满屏的代码报错，到最终看到 23 个完美的 Hi-C 互作方块，整个组装与注释的周期充满了挑战与成就感。工欲善其事，必先利其器，合理搭配硬件算力与最前沿的生信算法（如深度学习在注释中的应用），是这趟基因组之旅顺利通关的关键。

---

三、全流程代码

🧬 桂花 (SFZ) 双单倍型基因组组装与注释全流程记录

01. 数据准备 (HiFi 数据提取)

使用 NVIDIA Parabricks 整合包中的 bam2fq 将 BAM 转换为 FASTQ。（相较于samtools bam2fq 或 samtools fastq快10倍左右）。

cd /home/vensin/workspace/hifiasm/SFZ/00.data/02.hifi

# 利用 Parabricks bam2fq 转换 BAM 为 FASTQ
docker run --rm --gpus all \
  --volume $(pwd):/workdir \
  -w /workdir \
  nvcr.io/nvidia/clara/clara-parabricks:4.7.0-1 \
  pbrun bam2fq \
  --in-bam SFZ-yL.bam \
  --out-prefix SFZ-yL.hifi \
  --out-suffix _single.fastq.gz

02. 基因组初组装 (hifiasm)

整合 HiFi 和 Hi-C 数据，直接组装出双单倍型（Hap1 & Hap2）的 Contig。

cd /home/vensin/workspace/hifiasm/SFZ/assembly_run

# 运行 hifiasm (28线程)
nohup hifiasm -o SFZ_yL.asm \
  -t 28 \
  --h1 ../00.data/03.hic/SFZ-y_R1.fq.gz \
  --h2 ../00.data/03.hic/SFZ-y_R2.fq.gz \
  ../00.data/02.hifi/SFZ-yL.hifi_single.fastq.gz > run_hifiasm.log 2>&1 &

# 将 GFA 组装图格式转换为 FASTA 序列
gfatools gfa2fa SFZ_yL.asm.hic.hap1.p_ctg.gfa > SFZ.hic.hap1.p_ctg.fasta
gfatools gfa2fa SFZ_yL.asm.hic.hap2.p_ctg.gfa > SFZ.hic.hap2.p_ctg.fasta
gfatools gfa2fa SFZ_yL.asm.hic.p_ctg.gfa > SFZ.hic.p_ctg.fasta

03. 染色体挂载前期准备 (Hi-C 比对)

为 HapHiC 准备过滤后的 BAM 文件。

# 激活环境
conda activate haphic
cd /home/vensin/workspace/hifiasm/SFZ/assembly_run

# 构建 BWA 索引
bwa index SFZ.hic.hap1.p_ctg.fasta
samtools faidx SFZ.hic.hap1.p_ctg.fasta
bwa index SFZ.hic.hap2.p_ctg.fasta
samtools faidx SFZ.hic.hap2.p_ctg.fasta

# Hi-C 数据比对 (Hap1 & Hap2)
bwa mem -5SP -t 28 SFZ.hic.hap1.p_ctg.fasta \
    ../00.data/03.hic/SFZ-y_R1.fq.gz \
    ../00.data/03.hic/SFZ-y_R2.fq.gz | \
    samtools view - -@ 28 -S -h -b -F 3340 -o hap1_HiC.bam

bwa mem -5SP -t 28 SFZ.hic.hap2.p_ctg.fasta \
    ../00.data/03.hic/SFZ-y_R1.fq.gz \
    ../00.data/03.hic/SFZ-y_R2.fq.gz | \
    samtools view - -@ 28 -S -h -b -F 3340 -o hap2_HiC.bam

# 过滤低质量比对信号
/home/vensin/software/HapHiC/utils/filter_bam hap1_HiC.bam 1 --nm 3 --threads 28 | \
    samtools view - -b -@ 28 -o hap1_HiC.filtered.bam
/home/vensin/software/HapHiC/utils/filter_bam hap2_HiC.bam 1 --nm 3 --threads 28 | \
    samtools view - -b -@ 28 -o hap2_HiC.filtered.bam

04. 自动化染色体挂载 (HapHiC)

# 运行 Hap1 挂载
conda activate haphic
cd /home/vensin/workspace/hifiasm/SFZ/scaffolding_run/HapHiC_hap1

haphic pipeline /home/vensin/workspace/hifiasm/SFZ/assembly_run/SFZ.hic.hap1.p_ctg.fasta \
    /media/vensin/软件1/hifiasm/SFZ/scaffolding_run/hap1_HiC.filtered.bam 23 --threads 28

# 运行 Hap2 挂载
conda activate haphic
cd /home/vensin/workspace/hifiasm/SFZ/scaffolding_run/HapHiC_hap2

haphic pipeline /home/vensin/workspace/hifiasm/SFZ/assembly_run/SFZ.hic.hap2.p_ctg.fasta \
    /media/vensin/软件1/hifiasm/SFZ/scaffolding_run/hap2_HiC.filtered.bam 23 --threads 28

05. 序列格式化 (修改 Fasta Header)

将 HapHiC 输出的 >groupX 批量重命名为正式的 >ChrXX 格式。

mkdir -p /home/vensin/workspace/hifiasm/SFZ/final_genomes

# 处理 Hap1
cat /home/vensin/workspace/hifiasm/SFZ/scaffolding_run/HapHiC_hap1/04.build/scaffolds.fa | \
awk '{
    if($0 ~ /^>group[0-9]+$/){
        num=substr($1,7); 
        printf(">Chr%02d\n", num)
    }else{
        print $0
    }
}' > /home/vensin/workspace/hifiasm/SFZ/final_genomes/SFZ_hap1_scaffolds.fa

# 处理 Hap2
cat /home/vensin/workspace/hifiasm/SFZ/scaffolding_run/HapHiC_hap2/04.build/scaffolds.fa | \
awk '{
    if($0 ~ /^>group[0-9]+$/){
        num=substr($1,7); 
        printf(">Chr%02d\n", num)
    }else{
        print $0
    }
}' > /home/vensin/workspace/hifiasm/SFZ/final_genomes/SFZ_hap2_scaffolds.fa

06. 组装质量综合评估

包含了二代数据 QV 评估（Yak）和基因组完整度评估（BUSCO）。

A. Yak QV 评估

mkdir -p /home/vensin/workspace/hifiasm/SFZ/evaluation/qv
cd /home/vensin/workspace/hifiasm/SFZ/evaluation/qv

# 构建 Illumina K-mer 数据库
yak count -k 31 -b 32 -t 28 -o SFZ_illumina.yak \
    /home/vensin/workspace/hifiasm/SFZ/00.data/01.2nd/raw_R1.fq.gz \
    /home/vensin/workspace/hifiasm/SFZ/00.data/01.2nd/raw_R2.fq.gz

# 评估双单倍体 QV 值
yak qv -t 28 SFZ_illumina.yak /home/vensin/workspace/hifiasm/SFZ/final_genomes/SFZ_hap1_scaffolds.fa > SFZ_hap1_yak.qv.txt
yak qv -t 28 SFZ_illumina.yak /home/vensin/workspace/hifiasm/SFZ/final_genomes/SFZ_hap2_scaffolds.fa > SFZ_hap2_yak.qv.txt

B. BUSCO 完整度评估 (核心执行命令示例)

conda activate busco

# 示例：针对 Hap2 最终序列测试有胚植物数据库
busco -i /home/vensin/workspace/hifiasm/SFZ/scaffolding_run/HapHiC_hap2/04.build/SFZ_Hap2_Final_Genome.fa \
      -l embryophyta_odb12 \
      -o BUSCO_SFZ_Hap2_embryophyta \
      -m genome \
      -c 28

07. 基因组重复序列屏蔽 (EDTA)

利用 EDTA 从头注释转座子 (TE) 并生成 Masked 基因组。

conda activate EDTA

# 强行将引起报错的代码 x[0] 替换为正确的 x.iloc[0]
sed -i 's/family = x\[0\]/family = x.iloc\[0\]/g' /home/vensin/anaconda3/envs/EDTA/share/TIR-Learner3/app/check_TIR_TSD.py

mkdir -p /home/vensin/workspace/hifiasm/SFZ/annotation/repeat_hap2
cd /home/vensin/workspace/hifiasm/SFZ/annotation/repeat_hap2

# 软链接最终基因组
ln -s /home/vensin/workspace/hifiasm/SFZ/scaffolding_run/HapHiC_hap2/04.build/SFZ_Hap2_Final_Genome.fa ./

# 后台运行 EDTA
nohup EDTA.pl \
    --genome SFZ_Hap2_Final_Genome.fa \
    --species others \
    --step all \
    --sensitive 1 \
    --anno 1 \
    --threads 28 \
    > edta_hap2_run.log 2>&1 &

---