Bionano数据既可以自己从头de-novo出图谱,也可以依据参考序列混合组装出图谱。
===从头组装====
-
AutoNoise + SplitBNX: 这一步会将bnx和参考的cmap文件进行比对,估算出噪声系数,然后把bnx进行拆分便与后续比对
-
Pairwse: 这一步进行molecules之间的两两比较,寻找overlap, 结果存放在"align"文件夹下
-
Assembly: 根据两两比对结果,通过OLC算法进行组装,结果在"contigs/exp_unrefined"下,合并后的文件为"EXP_UNREFINED.cmap",同时还会将"EXP_UNREFINED.cmap"和参考基因组的cmap进行比对,结果放在"contigs/exp_unrefined/exp_unrefined/alignref",此外还将拆分后Bnx文件和参考基因组的cmap文件进行比对,结果放在"contigs/alignmolvref"下
-
refineA: 第3步得到的图谱先会进行第一次优化输出结果在"contigs/exp_refineA", 这一步不会使用所有的原始数据,而是pairwise阶段用的质量比较好的分子,所以速度会快一些
-
refineB: 在第4步的基础上进行第二次的优化, 输出结果在"contigs/exp_refineB0"和"contigs/exp_refineB1".这一步会使用所有的输入原始数据,速度稍微慢一些。第一轮和第二轮的结果会将地覆盖度的区域进行打断,然后更新标记和标记的位置。
-
Extension and merge: 将上一步的contig回贴到参考基因组的map,进行延伸和合并,这一步可以迭代3-5次。中间结果在"contigs/exp_extensionX_X"和"contigs/exp_mrgX"
-
Final refinement: 最后一步的优化
-
SVDecetion: 在有参考基因组的前提下,最后还会寻找一些大规模的结构变异。
由Solve工具中Pipeline文件夹下的脚本pipelineCL.py进行控制。
python /gpfs03/home/jingjing/software/Solve_3.3.1/Pipeline/10252018/pipelineCL.py -T 40 -N 4 -f 1 -i 5 -b ../raw/output.bnx -l Assembly -t /gpfs03/home/jingjing/software/Solve_3.3.1/RefAligner/7915.7989rel -a /gpfs03/home/jingjing/software/Solve_3.3.1/RefAligner/7915.7989rel/optArguments_nonhaplotype_noES_noCut_saphyr.xml
注:Solve的好多脚本python脚本不能在python3下面运行,要切换到python2下。
或者加上参数-r
perl /gpfs03/home/jingjing/software/Solve_3.3.1/Pipeline/10252018/fa2cmap_multi_color.pl -i /gpfs03/bioinfo/20210607_K326/assembly_k326.polish_contig.fa -e cttaag 1 -o .
python /gpfs03/home/jingjing/software/Solve_3.3.1/Pipeline/10252018/pipelineCL.py -T 40 -N 4 -f 1 -i 5 -r assembly_k326.polish_contig_CTTAAG_0kb_0labels.cmap -y -b ../raw/output.bnx -l Assembly -t /gpfs03/home/jingjing/software/Solve_3.3.1/RefAligner/7915.7989rel -a /gpfs03/home/jingjing/software/Solve_3.3.1/RefAligner/7915.7989rel/optArguments_nonhaplotype_noES_noCut_saphyr.xml
注: XML文件命名解释和选择
irys/saphyr: 数据来源仪器
DLE1: DEL1标记系统
human: 物种是人类
BG: big genome. 大于5G,主要是优化内存使用
noES: no extend and split: 即便等位基因里有超过30kbp的结构变异,也不要将他们分开
haplotype/nonhaplotype: 单倍型优化指的是将那些含有超过500bp或者更大的SV差异的等位基因进行分开,不推荐用于非人类基因组组装是用haplotype,这会导致组装结果过度碎片化,组装的基因组会变大.
nocut: 不对CMPR(complex multipath regions)进行拆分,所谓的CMPR指的是长度超过130kbp的高度重复序列,因为相似度过高,组装的时候不知道如何处理。
看网上大咖们的推荐是:对于非人类的物种,
除非是小鼠的SV缺失,大部分情况都用nonhaplotype,用于后续的Hybridscaffold Pipeline(HS);要noES;是否cut看情况而定。大部分人喜欢不cut。
最后输出结果是"contigs/exp_refineFinal1/EXP_REFINEFINAL1.cmap"。
而结果好坏则要看"exp_informaticsReportSimple.txt",两个核心标准。
一般来说:Bionano图谱应该占原来的物理图谱的90%以上。
Bionano图谱的N50 会由于物种不同有很大差异,动物一般在1M以上,我们是植物的结果,N50是33M。DLE系统差异更明显。
=======混合组装图谱========================
一般来说,包含下面几步:
为序列数据产生 insilico 图谱
将序列和Bionano基因组图谱进行比较,找到两者之间的冲突并尝试解决
将不冲突的图谱合并成hybrid scafold
在序列图谱和hybridscaffold之间形成联配
得到scaffold的AGP和FASTA文件。
perl /gpfs03/home/jingjing/software/Solve_3.3.1/HybridScaffold/10252018/hybridScaffold.pl -n /gpfs03/bioinfo/20210607_K326/assembly_k326.polish_contig.fa -b EXP_REFINEFINAL1.cmap -c hybridScaffold_config.xml -r /gpfs03/home/jingjing/software/Solve_3.3.1/RefAligner/7915.7989rel/RefAligner -o Hybrid_assembly -B 2 -N 2 -f
注意:
这个命令对perl的版本有要求,5.10.X and 5.14.X and 5.16.X and 5.18.X。(哎,版本还太高了,重新建了个环境)
另外:需要修改XML其中fasta2cmap的enzyme部分。
可以看出,混合组装后的N50为64M,scaffold数量为3996。
本文使用 文章同步助手 同步
网友评论