使用ALLMAPS进行辅助组装

简介

在从头组装过程中，确定基因组的scaffolds/contig的顺序和朝向是重建染色体非常关键的一步。这一步可以由多种辅助组装策略完成：例如遗传图谱, Hi-C, BioNano光学图谱，10X Chicago 。

一个物种可能会有多个遗传图谱，可以是不同项目中的不同定位群体结果，可以是不同软件如R/QTL, MSTMAP和JOINMAP的分析结果。遗传图谱会因重组率，偏分离(segregation distortion) , PAV(presence-absence variation)和染色体比对多态位点不同而发生变化。每一种图谱都能够提供不同的证据，举个例子，一个scaffold可能在一个图谱中无法被锚定，但是在另一个图谱中可以进行锚定，将这些图谱进行整合就能提最后染色体组装的精确度。

如果只用一个图谱，对scaffold进行排序只是计算量大一点而已，你需要根据图谱中分子标记在每个scaffold的平均距离进行排序就行。ALLMAPS, 正如名字说的那样，就是能够使用所有的图谱证据的工具，它能够计算scaffold的朝向，使其和已有的图谱的共线性关系最大化。他有如下亮点：

• 可重复性：清晰的可计算目标使得让多种输入图谱的共线性关系能够最大化
• 灵活性：允许为输入的图谱设置权重，更好的处理冲突
• 强大性：能使用多种遗传图谱，只需要做最小的转换
• 易用性：让基因组构建(FASTA和AGP输出)和基因组升级(CHAIN输出)流程化

通常，解决scaffold顺序和朝向(OO)是一种NP问题。因为基因组组装和遗传图谱中都可能存在一些错误，我们的目标只能是找到一个近似解。ALLMAPS将该问题转换成旅行商问题，然后用遗传算法优化scaffold OO. 使用遗传算法优化是为了避免在局部最优上出现瓶颈。遗传图谱最常见的错误是倒置和异位(inversion and translocation)

安装可以用CONDA。

conda install jcvi

实际演示

先现在一个测试数据集，用于练习了解软件的总体流程。比较尴尬的是，这个数据集放在了Dropbox上，我们都知道这是一个国内不存在的网站，所以我把他转存到了我的腾讯微云上,链接：https://share.weiyun.com/5nwjljN .

第一步：准备输入文件。 首先你得要提供物理图谱和遗传图谱的对应关系，格式为

Scaffold ID, scaffold position, LG, genetic position

你可以在excel表格中进行操作，然后另存为csv为文件，如下所示。你可以发现有一些scaffold里就只有一个遗传标记进行对应，你可以思考下该scaffold最后会如何处理

当然，我们解压缩后的zip文件里就有这些内容，你可以用一些命令工具(less, head, tail)看下数据是如何存放的。

第二步: 将两个图谱进行合并，最后会得到一个权重文件(weights.txt)和输入的bed文件

python -m jcvi.assembly.allmaps merge JMMale.csv JMFemale.csv -o JM-2-test.bed

第三步：对权重文件"weights.txt" 进行调整。weights.txt默认每个输入的图谱的权重都是1。 作者有一个建议就是，你通过检查最后的报告和诊断图，有监督地来重新对每个遗传图谱进行权重赋值。

$ cat weights.txt 
JMFemale 1
JMMale 1

第四步： 对scaffold进行排序，搭建成准染色体水平

python -m jcvi.assembly.allmaps path -w weights.txt JM-2-test.bed scaffolds.fasta

结果解读

上面第四步会在当前文件夹下生成许多结果文件，

JM-2-test.agp
JM-2-test.bed
JM-2-test.chain
JM-2-test.chr.agp
JM-2-test.chr.fasta
JM-2-test.fasta
JM-2-test.fasta.sizes
JM-2-test.lifted.bed
JM-2-test.summary.txt
JM-2-test.tour
JM-2-test.unplaced.agp
JM-2-test.unplaced.fasta

可以分为如下几类

首先是组装后的基因组序列

• "JM-2-test.fasta"
• "JM-2-test.agp": 每个scaffold的顺序和朝向，用于上传Genbank
• "JM-2-test.chain": 用于新旧坐标间的转换，比如说你在之前坐标注释得到GFF文件就可以用lifOver转换到新坐标系下

然后可视化报告。每个染色体都会得到一个对应的pdf文件，可视化展示如下：左图主要关注交叉的线，表示某些marker存在矛盾。右图关注是否有单上升/下降趋势，图中的斜率反应的是物理距离(x轴)相对遗传距离(y轴)的变化，可以认为是重组率。低重组率不容易确定在同一个重组区间内的scaffold的位置和朝向

白色部分为可信区间，灰色部分是存疑区间。

除了默认的输出外，你还可以通过movie得到软件运行过程每次迭代后组装情况，不过要额外安装ffmpeg和parallel

python -m jcvi.assembly.allmaps movie -w weights.txt JM-2-test.bed scaffolds.fasta chr23

最后是总结性报告。尽管可视化会给你比较直观的结果，但是具体的参数还是要看JM-2-test.summary.txt. 例如遗传标记的密度，以及有多少序列被锚定到基因组上。

当然ALLMAPS还有一些比较高级的操作：

• 拆分嵌合contig：嵌合contig指的是一段区域能够比对到多个连锁群中或者染色体中，一个常见的来源就是不同染色体中的重复区域由于过于相似在组装的时候坍缩成了一个。ALLMAPS也提供split进行拆分。
• 估计gap长度： 默认会用100个N在填充两个scaffold连接的区域，方便Genbank识别其为未知区域。你可以根据遗传距离在不同物理位置上的对应关系，预测不同gap的近似大小然后进行填充，子命令是estimategaps。
• 在ALLMAPS中使用多种遗传图谱.

光学图谱在着丝粒区间表型不好，主要是里面的串联重复过多，使得限制性内切酶数目减少，遗传图谱中在这个区域依旧有比较多的标记

局限性：

• 标记数目和密度。过多计算量大，太小准确性低
• 基因组组装和建图软件的潜在错误。基因组组装的原始质量对ALLMAPS的影响非常大，尤其是高度片段化的结果可能会导致错误结果。
• 基因组结构特征也会有影响，例如染色体倒置、异位和片段重复。

参考资料：

• ALLMAPS: robust scaffold ordering based on multiple maps
• https://github.com/tanghaibao/jcvi/wiki/ALLMAPS