2022.03.24【基因组组装】|获取比对到参考基因组的contig序列

摘要

很久没有整理工作笔记了，一方面个人有些倦怠，另一方面国内国际发生的事都牵动着许多人，我也不例外。趁着今天项目不多，记录一下最近的解决方案。
上周遇到一个想检测测序样品中是否包含预期的细菌物种。使用nr数据库比对以及metaphlan3进行物种注释都找到了客户的预期物种。然后客户希望通过测序数据组装出一套基因组。要求是组装出来的contig必须是都比对上参考基因组的。这个虽然不难，而且方法还不少，比如可以用之前去宿主的思路保留比对序列。但是这次我用了另外的工具，之前没具体做过，因此进行一个记录。

工具与方法

主要使用blast和seqtk这两个工具，前者进行序列比对，后者进行序列提取。

操作方法

step.1 构建参考基因组数据库

makeblastdb是blast用来构建数据库的命令，-in为输入参考基因组序列，-dbtype为数据库类型，一般就核酸序列和蛋白序列两种。

makeblastdb -in ref_bilvae.fasta -dbtype nucl

step.2 比对序列

blastn是核酸序列比较核酸序列数据库，-query为比对序列，-perc_identity为相似度过滤，可设置0-100，这里设置的99，-outfmt选择输出格式。可输出格式很多，我试了几个，发现选择6号表头格式比较方便获取到query_id序列。

blastn -db ref_bilvae.fasta -query AOLQ_192_scaffolds.fasta -out result.txt -perc_identity 99 -outfmt 6

step.3 获取query_id

这一步比较简单，通过上一步得到的比对结果使用awk可以在第一列获取到query_id，uniq则可以对序列进行去重。

awk -F "\t" '{print $1}' result.txt |uniq > align_queryid.txt

step.4 获取比对序列

最后一步就要用的神奇的seqtk工具（国人李恒大神开发），使用subseq命令直接根据id号提取contig序列，非常方便快捷。

seqtk subseq AOLQ_192_scaffolds.fasta align_queryid.txt >align_contigs.txt

结果展示

最后看一下我们得到的序列

最长的contigs不方便展示，我把最后一些短片段贴了出来。

欢迎大家加vx：bbplayer2021进群交流生信分析。