2023.10.22

输入数据来自第三步：

3.HiFi + HiC 数据通过 hifiasm 组装嵌合体基因组

HiFi reads 经过 hifiasm 组装至 contigs 水平后，接下来将借助 HiC 数据，将一个个的 contig 定位到染色体上，组装成染色体水平的基因组。

HiC 数据之前在 hifiasm 组装 contigs 的时候已经用过一次了，当时的作用是为了让 hifiasm 能够更好的 phase 两套单倍型，期望得到更好的组装效果，那么现在就是HiC数据在组装基因组中发挥最重要作用的时刻：将 contigs 锚定到染色体上。用到的工具就是 **chromap + yahs。**这两个软件都是近两年发表的，相比与之前使用的 Juicer + 3dDNA 流程，组装效果可以说又快又好，一个普通的哺乳动物基因组，两到三天即可完成组装。

chromap 论文链接： https://doi.org/10.1038/s41467-021-26865-w

yahs 论文链接：https://doi.org/10.1093/bioinformatics/btac808

  1. chromap 进行高效 Hi-C 数据比对

    需要对 pri.hic.asm.hic.p_ctg.fa 文件建立索引：

    # 普通 fasta 文件索引
samtools faidx pri.hic.asm.hic.p_ctg.fa

# chromap 专属索引
chromap -i -r pri.hic.asm.hic.p_ctg.fa -o pri.hic.asm.hic.p_ctg.index

    HiC 数据比对：

    chromap --preset hic \\
	-r ./pri.hic.asm.hic.p_ctg.fa \\
	-x ./pri.hic.asm.hic.p_ctg.index \\
	--remove-pcr-duplicates \\
	-1 ./hic.clean_R1.fastq.gz \\
	-2 ./hic.clean_R2.fastq.gz \\
	--SAM \\
	-o aligned.sam \\
	-t 28

# -1, -2 两参数后为质控后的 HiC FASTQ 文件。
# --remove-pcr-duplicates 去除测序过程中的 PCR 重复冗余。

    对输出的 SAM 文件排序并转 BAM 输出：

    samtools view -bh aligned.sam | samtools sort -@ 28 -n > aligned.bam

# 按照read的名字进行排序和按照位置排序或未排序，对 scaffolding 的结果会有一些不同。
# 成功运行完成后可删除 aligned.sam 文件，减少存储占用

  2. yahs 进行 scaffolding

    初步进行 sacffolding：

    yahs -e GATC \\
	pri.hic.asm.hic.p_ctg.fa \\
	aligned.bam

# -e GATC 指的是 HiC 测序时使用的 Mbol 酶切位点:GATC
# 如果对 contigs 组装结果自信，也可加参数 --no-contig-ec，加参数后将不对 contigs 进行纠错，即不打断 contigs 组装。
# yahs.out_scaffolds_final.fa 即为初步组装的结果。

使用 yahs 软件安装包自带的工具 juicer + juicer_tools.1.9.9_jcuda.0.8.jar 对组装结果进行可视化：（juicer_tools 下载链接：https://github.com/aidenlab/juicer/wiki/Download#version-199）

juicer pre -a \\
	-o out_JBAT \\
	yahs.out.bin \\
	yahs.out_scaffolds_final.agp \\
	pri.hic.asm.hic.p_ctg.fa.fai

# yahs.out.bin, yahs.out_scaffolds_final.agp 即为上一步的输出文件
# -o out_JBAT 为这一步的输出前缀，其中输出的 out_JBAT.txt 为接下来的输入文件
# 输出的 out_JBAT.assembly 用于后续的 juicebox 查看
# 输出的 out_JBAT.liftover.agp 用于后续的二次矫正
# 此步骤运行输出的 log 信息中包含类似 <(echo "assembly 1416099967") 信息，
# 其中 1416099967 为基因组大小或者大小的二分之一。
# 因为现阶段 yahs v1.2a.1 不支持 2G 以上的基因组，故超过 2G 的基因组可缩放 1/2, 
# 后续在 JuiceBox 中查看时设置缩放比例为 2 即可。

java -Xmx100G -jar ~/software/juicer_tools.1.9.9_jcuda.0.8.jar pre \\
	out_JBAT.txt out_JBAT.hic <(echo "assembly 1416099967")
# 生成的 out_JBAT.hic 文件与out_JBAT.assembly 可一起导入 JuiceBox 中查看组装结果

在 JuiceBox 中查看，若组装没错误，则初步组装结果 yahs.out_scaffolds_final.fa 即为 Scaffolding 的最终版本，若发现有组装不合适的地方，需要手动调整后，导出 out_JBAT.review.assembly 文件进行二次矫正

**juiceBox 手动调整，教程如下：**导出纠正后的 out_JBAT.review.assembly 文件

Juicebox Assembly Tools · aidenlab/Juicebox Wiki (github.com)

juicer post 进行二次矫正：

juicer post \\
	-o out_JBAT \\
	out_JBAT.review.assembly \\
	out_JBAT.liftover.agp \\
	pri.hic.asm.hic.p_ctg.fa

输出的 out_JBAT.FINAL.fa 即为最终的 scaffolding 结果