本文只是按照自己的需求翻译了HaploMerger2提供的手册部分内容。HaploMerger2的帮助文档写的非常好，一定要花点时间去读啊！
HaploMerger2的分析流程如下
- 重建单倍体组装中的等位基因关系
- 检测并纠正二倍体组装中的错连(mis-join)
- 重建2个单倍型组装
- 进一步对单倍型组装进行桥接
- 检测并移除串联错配
- 补洞
图片流程:
对于如何使用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
hm.batchA2是利用ChainNet算法创建互为最优的单覆盖(reciprocally-best single-coverage)全基因组联配
sh ./hm.batchA2.chainNet_and_netToMaf contigs_wm
这两步都可以通过修改其中的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
这一步可以进行调整的参数如下
- 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.
