HaploMerger2: 从高杂合二倍体基因组组装中重建单倍型

本文只是按照自己的需求翻译了HaploMerger2提供的手册部分内容。HaploMerger2的帮助文档写的非常好，一定要花点时间去读啊！

HaploMerger2的分析流程如下

重建单倍体组装中的等位基因关系
检测并纠正二倍体组装中的错连(mis-join)
重建2个单倍型组装
进一步对单倍型组装进行桥接
检测并移除串联错配
补洞

图片流程:

HaploMerger2的流程图

对于如何使用HaploMerger2获取高质量分型单倍型组装，作者从Canu和Falcon/Falcon-unzip的结果分别给了建议

如果是Canu，尽量避免基因组坍缩，尽可能获取完整的分离的二倍型组装；然后用HaploMerger2获取高质量的参考基因组和可选单倍型组装；然后用HapCUT2或者其他分型工具基于参考单倍型组装去获取高质量单倍型组装，
对于Falcon结果，将Falcon的输出结果p-tigs和a-tigs进行合并，然后将结果传给HaploMerger2

软件安装

在https://github.com/mapleforest/HaploMerger2/releases下载HaploMerger2，解压缩后会得到很多的文件夹

bin: HaploMerger2的核心程序
chainNet
gapCloser
Gmaj
lastz
SSPACE-STANDARD
winMasker
project_template: 配置文件

project_template里的运行脚本用的都是相对路径，也就是../bin/软件名, 这其实很不方便进行脚本调用，因此，推荐用下面的方法将../bin替换成你的实际的bin路径。（如果担心操作失误，可以备份project_template)

# 备份
cp -r project_template project_template_bck
cd project_template
# 确认要被更改的信息
grep '\.\.' hm*
# 用sed进行更改
sed -i 's/\.\./\/opt\/biosoft\/HaploMerger2/'  hm*

注: 我的路径是/opt/biosoft/HaploMerger2，你需要按照自己的情况进行调整。

最后还需要保证lastz, chainNet, winMasker, SSPACE, GapCloser这些软件应该在环境变量中。

项目实战

最低要求: 基因组压缩后的文件。

因为我没有处理过10k, 20k, 50k文库用于scaffold的项目，而且D流程会删除基因组序列，所以本次实战只运行A, B这2个流程.

在HaploMerger2目录下新建一个项目文件，命名为example,之后下载分析所用的数据。然后将你的基因组序列移动到example里，例如 contig.fa

mkdir project_your_name
cd project_your_name
cp /path/to/your/contig.fa .

从project_template拷贝流程对应的脚本和配置文件

cp ../project_template/hm.batchA* .
cp ../project_template/hm.batchB* .
cp ../project_template/hm.batchD* .
cp ../project_template/all_lastz.ctl .
cp ../project_template/scoreMatrix.q .

第一步：对重复序列进行软屏蔽

这是必须要做的一步，能提高后续分析的准确性。推荐用windowmasker。

构建能重复使用的屏蔽文库

windowmasker \-checkdup true \-mk_counts \-in contigs.fa \-out masking-library.ustat \-mem 6500

然后用构建的文库对基因组进行重复序列软屏蔽，也就是将重复序列从大写改成小写。

windowmasker \-ustat masking-library.ustat \-in contigs.fa \-out contigs_wm.fa \-outfmt fasta \-dust true

然后将输出文件进行压缩

gzip contigs_wm.fa

第二步: 移除组装结果中主要的错连

如果你发现运行脚本的时候只用了几秒钟就结束了，那么基本上就是你的环境没有配置好。

hm.batchA1是对文件进行拆分，然后用LASTZ进行all-vs-all的全基因组比对。比对结果会存放在contigs_wm.contigs_wmx.result/raw.axt(contigs_wm是前缀), 根据你设置的特异性不同，占用的空间大小也不同，在运行完第二个脚本(A2),可以删除。

sh ./hm.batchA1.initiation_and_all_lastz contigs_wm

运行A1

hm.batchA2是利用ChainNet算法创建互为最优的单覆盖(reciprocally-best single-coverage)全基因组联配

sh ./hm.batchA2.chainNet_and_netToMaf contigs_wm

运行A2

这两步都可以通过修改其中的threads参数来提高运行速度，其中LASTZ每个线程要求约1G内存，而chainNet需要至少4-8G内存。

如果为了降低非同源联配，可以修改hm.batchA1中的identity参数，默认是80.如果你的基因组杂合度也就是4~5%左右，那么设置88%~92%能够获得更加严格的过滤。

这两步有两个非常重要的配置文件，all_lastz.ctl和scoreMatrix.q。其中all_lastz.ctl控制LASTZ，例如你可以将--step=1改成--step=20来提高LASTZ的运行速度。默认的参数其实已经挺不错了，如果想要修改的话，建议先去看看LASTZ的手册。

scoreMatrix.q里面记录的是核酸替换的得分矩阵，用于LASTZ和chainNet中。默认的得分矩阵来自于Chinese amphioxus的二倍体组装(GC 41%, 杂合度4~5%)。为了更好的区分等位基因和非等位基因的联配，建议参考后面的如何自定义LASTZ得分矩阵部分进行构建。

hm.batchA3会调用5个Perl程序，完成5个任务

预过滤互为最优的全基因组chainNet联配
将全基因组联配转成有向图(DGA graph
使用贪婪算法遍历有向图，寻找错连
将序列在错连的位置进行打断
基于有向图输出二倍体组装

sh hm.batchA3.misjoin_processing contigs_wm

运行A3

这一步可以进行调整的参数如下

scoreSchme: 计分规则，目前使用默认的score即可，分数越高说明联配质量越高
filterScore_1, escapeFilter, filterScore_2和NsLCsFilter: 大部分低得分的联配可能都是假的，应该被归为噪音。有三个参数用于进行控制过滤规则。

对于第二步，作者建议迭代2~3次，保证尽可能的找到错连的部分。同时还建议尝试下调aliFilter和overhangFilter（从50kb到40kb即可，低于30kb要十分小心）

假如你迭代了三次，那么每次的结果应该是contigs_wm > contigs_wm_A > contigs_wm_A_A > contigs_wm_A_A_A

第三步: 创建二倍型组装

先运行hm.batchB1.initiation_and_all_lastz和hm.batchB2.chainNet_and_netToMaf

sh hm.batchB1.initiation_and_all_lastz contigs_wm_A_A_A
sh hm.batchB2.chainNet_and_netToMaf contigs_wm_A_A_A

前两个脚本和之前的batchA的前两个脚本类似，除了两点:

二倍体输入文件应该是第一步的输出结果
chaiNet联配会输出MAF文件，可以用Gmaj查看

运行hm.batchB3.haplomerger

sh hm.batchB3.haplomerger contigs_wm_A_A_A

第三个脚本做了以下四件事情

预过滤互为最优的全基因组chainNet联配
将全基因组联配转成有向图(DGA graph)
使用贪婪算法遍历和线性化DGA
基于线性化的DGA输出参考单倍型和可选单倍型

hm.batchB3和hm.batchA3的参数类似。作者发现，提高filterScore_2 和escapeFilter能够避免将原本不是等位基因的两条序列进行合并，但是却会导致单倍型组装结果丢失正确的联配信息，从而导致包含了更多的冗余序列。部分的联配中主要是N和重复序列，使用NsLCsFilter能够过滤这些联配。

同时HaploMerger还能够将2个存在重叠区域(overlapping region)的序列进行合并，

合并

参数minOverlap控制合并操作的最小联配长度（不含N, gaps和InDels）。默认设置是0，也就是将任何可能合并的序列都进行合并。由于低分联配已经被filterScore_2过滤，所以最小联配长度实际上是大于filterScore_2。提高minOverlap会让结果更加可信，不过肯定也会损失可能的连接。

因为任何合并都会将两个单倍型合并，或者你不想混合单倍型或在单倍型间产生交换(generate switches between haplotypes)，那么你可以设置minOverlap= 99999999

这一步会输出三个重要文件，分别是

optiNewScaffolds.fa.gz: 有联配支持的新参考序列
optiNewScaffolds_alt.fa.gz: 有联配支持的可选参考序列
unpaired.fa.gz: 没有联配支持的序列

第四个脚本是优化没有联配支持的序列

sh hm.batchB4.refine_unpaired_sequences contigs_wm_A_A_A

unpaired.fa.gz里存放的是没有联配支持的序列，有以下几个原因会导致该情况

一些序列在二倍型组装中没有对应的等位序列
等位序列/单倍型坍缩成单个序列
一些真实的联配由于信号太弱被过滤
一些联配包含N和重复序列，导致它被过滤
一些联配和串联重复，倒置和易位，所以在线性化的DGA中被忽略
在图遍历是，部分序列会被过滤（例如序列中的无法联配自由末端）

考虑到这些情况，我们会认为这些序列中大部分是冗余序列。所以bactchB4脚本就是找到这些冗余序列，然后进行过滤。该脚本会将unpaired.fa.gz和optiNewScaffolds.fa.gz进行比对，然后移除其中潜在冗余序列。这会得到新的fasta文件，即unpaired_refined.fa.gz。几个可供调整的参数

threads: 控制线程数
identity: 控制特异性和敏感性
maskFilter: 过滤只有N和重复序列的scaffold
redudantFilter: 过滤在optiNewScaffolds.fa.gz有对应等位序列的scaffold

最后一步就是合并。将optiNewScaffolds.fa.gz 与 unpaired_refined.fa.gz合并得到参考单倍型contigs_wm_A_A_A_ref.fa.gz, 将optiNewScaffolds_alt.fa.gz 与 unpaired_refined.fa.gz合并得到contigs_wm_A_A_A_alt.fa.gz

sh hm.batchB5.merge_paired_and_unpaired_sequences contigs_wm_A_A_A

最后的结果是contigs_wm_A_A_A_ref.fa.gz和contigs_wm_A_A_A_alt.fa.gz

如何自定义LASTZ得分矩阵

作者为了方便我们推断LASTZ_D的得分矩阵，封装了一个脚本 lastz_D_Wrapper.pl

我们需要根据基因组大小将基因组序列分成2个部分，一部分只有5~15%的序列，而另一部分则是剩下的序列。一般而言，你应该选择最长的序列用于第一部分。

lastz_D_Wrapper.pl --target=part1.fa.gz --query=part2.fa.gz --identity=90

这里的参数主要是--identity=90, 表示lastz_D只会使用大于90%相似度的联配，该值越高表示严格度越高。作者推荐当杂合率高在4%~5%的时候选择90%，在1%~3%的时候设置为94%~96%。

运行完之后，将part1.part2.raw.time.xx_xx.q的内容复制到scoreMatrix.q中。

注1: lastz_D速度很慢而且不支持并行，因此只要用10%以内的序列作为target就行。序列多了反而找到等位基因的概率低了。

注2: 有些时候lastz_D还会诡异的停住，然后输出奇怪的结果。此外，不同的序列还会有不同的推断结果。因此，part1.fa.gz可能要选择不同的序列，然后运行程序。之后，你可以从不同的得分矩阵中筛选结果。

注3: 这一步和基因组复杂度密切相关，杂合度越高，运行时间越长。

参考文献

如果在分析中用到了HaploMerger2，请引用下面的文章

Shengfeng Huang, et al. HaploMerger2: rebuilding both haploid sub-assemblies from high-heterozygosity diploid genome assembly. Bioinformatics. 2017.
Shengfeng Huang, et al. HaploMerger: reconstructing allelic relationships for polymorphic diploid genome assemblies. Genome Res. 2012, 22(8):1581-1588.