2020.10.13丨Blast本地比对及可视化

• 这次接到一个客户的需求，他自己组装了一个大肠杆菌的序列（之前客户用它比对，所以简称A基因组），并且想用这个基因组对几个样品进行比对和注释，但注释结果不太好，所以问我们能不能比对用这个客户组装的基因组，注释用NCBI上的标准的参考基因组（B基因组）的注释文件。
• 我的第一反应就是不行，即使物种类似，不同基因组肯定在长度、位置上有差异，导致注释信息不匹配。本来都打算直接用B基因组给客户从比对开始做。但领导提醒我，开始之前还是需要先回复客户，得到客户的意见反馈后再开始。如何给客户回答这个问题，就是我本次的工作。
• 首先，在NCBI上检索，在核苷酸一栏检索得到客户组装的A基因组\B基因组的序列信息
• 查看Genbank基本信息，发现组装序列为4546849 bp，而标准基因组有4641652 bp 。10万bp的差异，基本不能用了。接下来就是要对两个序列进行比对，看看有没有可能缺失片段是连续的，如果正好连续且都在尾部，那前面注释信息似乎也可以用，不过这个概率太小了。
• 本来打算线上比对，结果软件提示序列过长，不支持，因此，我们需要下载序列到本地进行比对。
• 首先，下载A、B基因组序列，选择的是fasta格式
• 接着，安装Blast，这里使用miniconda3进行快速简易安装

  • $ source ~/miniconda3/bin/activate
    $ conda install blast

     

• 安装好之后，通过帮助页面确认可以使用。
• blast的比对原理是query序列去和数据库进行比对，由于我这里是两条相似的基因组，因此不需要单独下载其他数据库。这里，我们将标准参考基因组（B基因组）进行建库,建库的参数是 makeblastdb

  • makeblastdb -in B_genome.fasta -dbtype nucl -parse_seqids -out B_genome

    建库成功后会得到6个不同格式的文件

• 然后，我们就可以以B基因组为数据库，对A基因组进行比对了。这里要说明的是，blast生成的结果文件有多种格式，由于后面可视化工具有格式要求，因此我生成的是XML格式第五个版本。

  • blastn -db B_genome -query A_genome.fasta -out result.xml -outfmt 5

     

• 得到结果
• 
• 虽然得到了结果，不过看起来太费劲，到底哪一段比对上了，哪一段没有比对上，都需要更直观地展示。这里就要推荐大家使用一个在线可视化工具，Kablamm
• 他的功能主要就是展示blast的比对结果，这里，我们直接选择load results，可以得到上传及参数选择解面
• 这里要求上传文件格式为XML，之前我们比对时已经设置好，所以不用再重新运行。至于参数选择该如何选择，这要看你的比对结果，我给大家分别展示一下默认参数，以及修改参数后的图表结果。
• 这是参数较低时的结果，部分序列会匹配到其他位置上。有颜色的地方是比对上的，白色就是没比对上的，颜色越深，bit score分数越高。看到的白线实际是两条。两条白线中间代表比对的序列，不过由于片段较小，大图看起来就只有一条白线，可以明显看到筛选参数较低会出现错配情况
• 
• 设置合理的筛选参数后（我这里提高了bit score值），可以看到大部分序列都是相对应的。但中间仍有很多没有匹配上的序列片段。
• 
• 到这里，我们基本可以解释两条基因组即使物种类似，也不能用A基因组比对后再用B基因的注释文件去进行注释。