文章目录
摘要工具与方法操作方法step.1 构建参考基因组数据库step.2 比对序列step.3 获取query_idstep.4 获取比对序列结果展示摘要
很久没有整理工作笔记了，一方面个人有些倦怠，另一方面国内国际发生的事都牵动着许多人，我也不例外。趁着今天项目不多，记录一下最近的解决方案。
上周遇到一个想检测测序样品中是否包含预期的细菌物种。使用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：bbplayer进群交流生信分析。
如果觉得《.03.24【基因组组装】|获取比对到参考基因组的contig序列》对你有帮助，请点赞、收藏，并留下你的观点哦！
