原创内容,欢迎转载,转载请注明出处
注意事项
-
序列比对是生物信息学最基础的工作,很多工作都是基于序列比对进行的。序列比对分为局部比对和全局比对,局部比对最经典的就是Blast,也就是本文介绍的内容,全局比对多用于多序列比对,比如说clustal、Muscle、Mafft等工具,这部分内容可以查看网站其他部分的教程。
-
局部比对最常用的工具就是Blast,此外还有很多其他的工具,比如说速度远远快于Blast,用于进行蛋白序列比对的Diamond,中国科学院北京基因组研究所基因组科学与信息重点实验室“百人计划”章张研究员开发的ParaAT,另一个代替blast的工具MMSeq2,用于结构域比对的Hammer,三代测序长reads比对的Minimap2,用于二代测序reads比对的BWA、Bowtie2等。
-
虽然序列比对的工具有很多,但是对于我们这样一个非生物信息学团队而言,大家能够接触到的数据量并不大,Blast完全可以解决大部分的序列比对需求,所以只需要学会并熟悉Blast的基本功能即可。
-
同样,由于比对的需求不大,因此完全可以将Blast安装到本地来使用,由于windows10和windows11的命令行工具越来越好用了,因此,建议大家将Blast安装到本地,在自己的电脑上进行Blast比对。当然也可以到实验室的服务器上面进行比对,比对的结果是完全一致的。
-
在版本上面大家尽量安装最新的版本。
一、Blast在windows上安装使用
1、程序下载
访问blast本地软件包链接 blast_latest ,下载适合自己系统的blast版本。
如果下载速度太慢,可以到实验室的服务器上下载,服务器上面有两个版本分别是:ncbi-blast-2.2.28±win64 和 ncbi-blast-2.14.0±x64-win64,两个版本对于我们这些初级的使用者而言没有区别。
2、安装流程
如果下载以win64.exe结尾的文件,则需要双击安装。如果是以win64.tar.gz结尾则是绿色版本,解压后就可以使用。
安装(或者解压缩)后会生成bin和doc两个子目录,其中bin是程序目录,doc是文档目录,这样就安装完成。
路径中有中文是否会影响运行并不清楚,但是按照习惯还是尽量不要出现中文。
3、添加环境变量
建议直接在百度上搜索方法,这里简单的说一下步骤。
右键点击“我的电脑”-属性,然后点击“高级系统设置”,点击“环境变量”,在用户变量下方点击“Path”这一行,点击下面的“编辑”,点击“新建“,把blast安装的的bin目录的路径拷贝粘贴到里面就行了,最后一路点击”确定。就可以了。
4、测试是否安装配置成功
在windows中,点击“Win+R”在出现的运行窗口中输入cmd或者powershell都可以打开命令行。建议使用powershell,因为powershell的命令跟linux系统的命令一样,如果熟悉了powershell的命令,在后面如果需要使用linux服务器的时候更好上手,如果使用的是cmd命令行,则需要熟悉dos命令,增加了学习成本。
如果安装配置完成,只要打开cmd命令窗口和powershell命令窗口,输入“makeblastdb”,如果出现USAGE,并且出现用法,则说明安装成功了(见下图),如果直接显示Error,则说明安装失败,大概率是环境变量配置的问题。
二、Blast使用:本地数据库的建立(数据库格式化)
Blast的详细用法和详细参数请参考网上其他的Blast使用手册,Blast的高级用法我们通常使用不到,我们最常用的就是用blast来寻找同源序列。因此这里就只介绍最简单的使用方法。
进行Blast比对的时候,不能直接比对fasta格式,必须将想要作为比对对象的一方进行格式化(制作成为数据库),用另一方(fasta格式)与数据库进行比对。
这时我们首先需要使用到makeblastdb程序来格式化。
makeblastdb程序的使用方法如下:
makeblastdb -in dna-database.fasta -out dnadb -dbtype nucl
这里的dna-database.fasta是被比对的数据,里面的序列必须是fasta格式。
dnadb是生成的数据库的名字,这个名字可以随便起,比对的时候会用到这个名字。
-dbtype nucl 代表的是数据库是DNA序列,nucleotide的前面4个字母,如果数据库序列是蛋白序列则后面应该跟着prot(protein前面4个字母)。
下面就是构建蛋白数据库的命令:
makeblastdb -in pro-database.fasta -out prodb -dbtype prot
这里的pro-database.fasta是被比对的数据,里面的序列必须是fasta格式。
prodb是生成的数据库的名字,这个名字可以随便起,比对的时候会用到这个名字。
-dbtype prot 代表的是数据库是DNA序列,nucleotide的前面4个字母,如果数据库序列是蛋白序列则后面应该跟着prot(protein前面4个字母)。
三、Blast使用:DNA序列与DNA序列比对 blastn
实例:
我们通常用来进行ITS测序的引物是ITS1和ITS4,当我们送到公司测序的时候,由于测序质量不同会导致无法完全准确的获取ITS1引物和ITS4引物之间的序列。(具体原因可以查看:食用菌分类与鉴定:ITS测序与系统进化分析)
但是如果我们有单核菌丝的全基因组测序数据,我们就可以用序列比对的方法寻找目标菌株的完整的ITS序列。
本教程以来自于平菇的ITS序列,用来寻找猪肚菇上面的ITS序列。
下面是blastn工具进行序列比对的过程,首先用makeblastdb工具格式化基因组,生成数据库,再用blastn进行比对。-query后面跟着查询序列的文件名称,序列需要是fasta格式的文本文档,文件以什么结尾都可以。-db跟着格式化好的数据库的名字,-out后面跟着要输出的文件名称;-outfmt 6代表利用第6种格式输出,也就是表格格式(tabular),生成的实际上是文本文档,即使这里使用xlsx结尾,生成的还是文本文档,只不过以xlsx结果可以直接双击用表格打开,比较省事。其他输出格式自己百度学习。
-evalue 1e-5是设置一个比对的阈值,只保留可信度高于这个阈值的结果,具体应该用什么值自己设置。
# 创建数据库
makeblastdb -in Pgig.genome.fa -out genome -dbtype nucl
# 利用blastn进行DNA序列与DNA序列比对
blastn -query Post-ITS.fa -db genome -out result.xlsx -outfmt 6 -evalue 1e-5
双击result.xlsx打开比对结果,可以看到表格一共12列,在不同的contig上面有很多个ITS的拷贝。
每一列的含义如下:
query id:查询序列ID标识;
refer id:参考序列ID标识;
identity (%):序列比对的一致性百分比;
alignment length:符合比对的比对区域的长度;
mismatches:比对区域的错配数;
gap openings:比对区域的gap数目;
q.start:比对区域在查询序列(query id)上的起始位点;
q.end:比对区域在查询序列(query id)上的终止位点;
s.start:比对区域在参考序列(refer id)上的起始位点;
s.end:比对区域在参考序列(refer id)上的终止位点;
e-value:比对结果的期望值,将比对序列随机打乱重新组合,和数据库进行比对,如果功能越保守,则该值越低;
bit score:比对结果的bit score值;
四、Blast使用:蛋白序列与蛋白序列比对 blastp
小任务:
萎锈灵carboxin抗性是食用菌遗传转化中常用的一个抗性基因,利用carboxin抗性进行筛选假阳性率低,筛选效果好。carboxin抗性就是succinate dehydrogenase蛋白复合体上面的iron-sulfur亚基的突变体。琥珀酸脱氢酶是三羧酸循环的重要功能部分,由四个亚基构成,分别是黄素蛋白(Fp, SdhA)、铁硫蛋白(Ip, SdhB) 和两种嵌膜蛋白( SdhC和SdhD) 。当琥珀酸脱氢酶B亚基的一个组氨酸(His)突变为亮氨酸(Leu)之后对于萎锈灵敏感性降低,但是仍然能够保持琥珀酸脱氢酶的活性。
本教程以来自于杏鲍菇的CbxR cassette序列作为query序列,用来寻找大球盖菇上面的同源序列。该序列可以利用点突变构建大球盖菇专用的CbxR基因。
由于核酸序列并不保守,蛋白序列相对保守,因此我们应该用杏鲍菇的iron-sulfur基因来对比大球盖菇的iron-sulfur序列,再利用多序列比对来寻找应该突变哪一个位点。
下面是我从网上找的杏鲍菇Pleurotus eryngii ATCC 90797的琥珀酸脱氢酶B亚基的基因序列和蛋白序列,如果要改造成为cbxR基因,需要修改的位点我已经在下面比较出来了。
>KAF9489857
ATGCAGGCGCTCACCTCCAGGTCGTTGGCTCGCTCATCCCGCTCGATTCGTGCTTTCTCCACCTCGTGTGGAAGGTGGCAGGCTGAGCCCCTTCAGAAGCCCGTTCTCCAGAAAGAATTCAAGATCTACCGTTGGAATCCGGATGAACCAGCCAAGAAGCCTCATCTCCAATCGTACACTATTGACTTGAACCAGACGGGCCCCATGATTTTGGATGCTCTTATCAAGATCAAGAACGAAATCGATCCTACGCTCACATTCCGTCGTTCGTGCAGAGAGGGGATTTGCGGCTCGTGCGCGATGAATATTGACGGACAGAACACGCTGGCTTGCCTGTGCCGAATTGACCGTAACGCCAGCAAGGACAGCAAGATCTACCCTTTGCCGCACATGTACATCGTGAAAGACCTCGTACCCGACCTCACGCTTTTCTACAAGCAGTACAAGTCCATCAAGCCTTACCTGCAGAACGACAATGTCCCCGAGAGGGAGCACCTCCAGTCGCCAGAGGACCGCAAGAAGCTGGACGGGATGTATGAATGCATCCTGTGCGCGTGCTGCAGCACATCGTGCCCCAGCTATTGGTGGAACCAAGATGAGTATTTGGGACCAGCTGCTTTGATGGCCGCGTATAGGTGGATTGCGGACTCACGAGATACGTATGGCGCACAACGCAAGGAACATTTCCAGAATGAGATGAGTTTGTTCCGCTGCCAC(CTC)ACCATCTTCAATTGCTCTCGCACTTGTCCAAAGGGCCTCAACCCTGCGAAAGCCATTGCAGAGATCAAGCTCGCGCTTGCTACGGAGTAA
>KAF9489857 iron-sulfur protein subunit
MQALTSRSLARSSRSIRAFSTSCGRWQAEPLQKPVLQKEFKIYRWNPDEPAKKPHLQSYTIDLNQTGPMILDALIKIKNEIDPTLTFRRSCREGICGSCAMNIDGQNTLACLCRIDRNASKDSKIYPLPHMYIVKDLVPDLTLFYKQYKSIKPYLQNDNVPEREHLQSPEDRKKLDGMYECILCACCSTSCPSYWWNQDEYLGPAALMAAYRWIADSRDTYGAQRKEHFQNEMSLFRCH(L)TIFNCSRTCPKGLNPAKAIAEIKLALATE*
这里需要做的就是利用杏鲍菇iron-sulfur蛋白序列,与大球盖菇的预测的所有的蛋白序列进行比较,找到大球盖菇琥珀酸脱氢酶的iron-sulfur亚基的蛋白序列,再找到对应的cds序列和基因组上面的位置。取ATG上游的~1500 bp序列和下游的~300 bp序列作为cbxR cassette,将其中的H突变为L就行了。
# 创建数据库
makeblastdb -in revised_proteome.fa -out proteome -dbtype prot
# 利用blastn进行DNA序列与DNA序列比对
blastp -query PecbxR.fasta -db proteome -out result.xlsx -outfmt 6 -evalue 1e-3
双击result.xlsx打开比对结果,可以看到比对到了一个蛋白,这个蛋白很有可能就是我们的目标蛋白。需要进一步用实验进行验证。
五、Blast使用:blastx、tblastn与tblastx
其他的几个blast程序的使用方法自学。
评论区