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

  1. 云栖社区>
  2. 博客>
  3. 正文

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

徐洲更 2019-05-23 21:43:07 浏览442
展开阅读全文

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

HaploMerger2的分析流程如下

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

图片流程:

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.ctlscoreMatrix.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个任务

  1. 预过滤互为最优的全基因组chainNet联配
  2. 将全基因组联配转成有向图(DGA graph
  3. 使用贪婪算法遍历有向图,寻找错连
  4. 将序列在错连的位置进行打断
  5. 基于有向图输出二倍体组装
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_lastzhm.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的前两个脚本类似,除了两点:

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

运行 hm.batchB3.haplomerger

sh hm.batchB3.haplomerger contigs_wm_A_A_A

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

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

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

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

合并

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

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

这一步会在contigs_wm_A_A_A.contigs_wm_A_A_Ax.result文件下输出三个重要文件,分别是

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

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

sh hm.batchB4.refine_unpaired_sequences contigs_wm_A_A_A

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

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

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

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

最后一步就是合并。 将optiNewScaffolds.fa.gzunpaired_refined.fa.gz合并得到参考单倍型contigs_wm_A_A_A_ref.fa.gz, 将optiNewScaffolds_alt.fa.gzunpaired_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.gzcontigs_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,请引用下面的文章

网友评论

登录后评论
0/500
评论
徐洲更
+ 关注