Archives

BLASTparser : 转换BLAST -m5/ -m0 格式成Table格式

一、 准备工作:

1.1 系统平台:这里使用的是Ubuntu Linux 12.04 LTS  64bit 版本;

1.2 Windows版本: 如果你说只有Windows,安装个虚拟机嘛

推荐VirtualBox  https://www.virtualbox.org/

网上有相当多的图文教程

比如:http://blog.csdn.net/tangyajun_168/article/details/7063448

1.3 BLAST 本地下载安装

wget ftp://ftp.ncbi.nih.gov/blast/executables/LATEST/ncbi-blast-2.2.28+-x64-linux.tar.gz

注意 blast 和 blast+ 在参数方面有很大的不同。而且blast 好像2.2.26 以后就没有更新,而blast+一直在更新,而且blast+可以自己定制(- outfmt 6 ,7, 10)输出格式。

这里测试都是使用blast+。

二、测试数据

2.1 参考序列:

使用The Consensus CDS (CCDS)数据库中的人类的氨基酸序列和核酸序列(20130430 版本),关于CCDS以后会在kbase 标签下介绍。

下载地址:ftp://ftp.ncbi.nlm.nih.gov/pub/CCDS/current_human/

2.2查询序列:

找一个比较小的测试集,选了人类的编码激酶的蛋白和核酸数据(kinome)。

http://kinase.com/kinbase/FastaFiles/Human_kinase_rna.fasta

http://kinase.com/kinbase/FastaFiles/Human_kinase_protein.fasta

三、执行BLAST

这里我把测试了BLASTp,BLASTn,BLASTx

因为自定义输出格式(仅限于table格式)大家用的不多,所以这里就给了几个示例。

3.1. Formatdb

../bin/ncbi-blast-2.2.28+/makeblastdb -in  CCDS_protein.20130430.faa  -title CCDS -dbtype prot  -taxid 9606  -out  CCDS_protein.20130430  -logfile CCDS_protein.20130430.log

../bin/ncbi-blast-2.2.28+/makeblastdb -in  CCDS_nucleotide.20130430.fna -title CCDS -dbtype nucl -taxid 9606  -out  CCDS_nucleotide.20130430  -logfile CCDS_nucleotide.20130430.log

3.2  BLASTP

../bin/ncbi-blast-2.2.28+/blastp  -query Human_kinase_protein-100.fasta  -db  ../ccds/CCDS_protein.20130430   -out  Human_kinase_protein-blastp-m5.xml   -evalue 1  -outfmt 5   -num_alignments 10   -num_threads 8

../bin/ncbi-blast-2.2.28+/blastp  -query Human_kinase_protein-100.fasta  -db  ../ccds/CCDS_protein.20130430   -out  Human_kinase_protein-blastp-m7.tbl   -evalue 1  -outfmt  “7  qseqid qlen slen qcovhsp sseqid bitscore score evalue pident qstart qend sstart send”   -num_alignments 10   -num_threads 8

../bin/ncbi-blast-2.2.28+/blastp  -query Human_kinase_protein-100.fasta  -db  ../ccds/CCDS_protein.20130430   -out  Human_kinase_protein-blastp-m0.txt   -evalue 1  -outfmt 0   -num_descriptions 10  -num_alignments 10   -num_threads 8

3.3  BLASTN

../bin/ncbi-blast-2.2.28+/blastn  -query  Human_kinase_rna-100.fasta  -db  ../ccds/CCDS_nucleotide.20130430   -out  Human_kinase-rna-blastn-m5.xml   -evalue 1  -outfmt  5   -num_alignments 10   -num_threads 8

../bin/ncbi-blast-2.2.28+/blastn  -query  Human_kinase_rna-100.fasta  -db  ../ccds/CCDS_nucleotide.20130430   -out  Human_kinase-rna-blastn-m7.tbl   -evalue 1  -outfmt “7  qseqid qlen slen qcovhsp sseqid bitscore score evalue pident qstart qend sstart send”    -num_alignments 10   -num_threads 8

../bin/ncbi-blast-2.2.28+/blastn  -query  Human_kinase_rna-100.fasta  -db  ../ccds/CCDS_nucleotide.20130430   -out  Human_kinase-rna-blastn-m0.txt   -evalue 1  -outfmt 0  -num_descriptions 10  -num_alignments 10   -num_threads 8

3.4  BLASTX

../bin/ncbi-blast-2.2.28+/blastx  -query  Human_kinase_rna-100.fasta  -db  ../ccds/CCDS_protein.20130430   -out   Human_kinase-rna-blastx-m5.xml   -evalue 1   -outfmt  5   -num_alignments 10  -num_threads 8 &

../bin/ncbi-blast-2.2.28+/blastx  -query  Human_kinase_rna-100.fasta  -db  ../ccds/CCDS_protein.20130430   -out   Human_kinase-rna-blastx-m0.txt   -evalue 1   -outfmt  0  -num_descriptions 10   -num_alignments 10  -num_threads 8 &

 

../bin/ncbi-blast-2.2.28+/blastx  -query  Human_kinase_rna-100.fasta  -db  ../ccds/CCDS_protein.20130430   -out  Human_kinase-rna-blastx-m7.tbl   -evalue 1  -outfmt “7  qseqid qlen slen qcovhsp sseqid bitscore score evalue pident qstart qend sstart send”   -num_alignments 10  -num_threads 8 &

四、-m 0 和 5格式的提取

由于-m 0和 –m5 格式的输出文件信息比较全面,所有一般我们会选择这个两个格式,根据需要提取信息,比如我们可以根据BLASTx 结果寻找一条核酸中的可能编码蛋白的区域,并直接给出对应的那段氨基酸序列。

在GiteCafe (链接:https://gitcafe.com/jameslz/BLASTparser)中 提供了两个代码:blastXmlParser.pl 和 blastTxtParser.pl

根据测试:

blastXmlParser.pl 可以提取BLASTx , BLASTp, BLASTn 的 –m 5 格式的结果,返回Table格式。

blastTxtParser.pl 可以提取BLASTx , BLASTp, BLASTn 的 –m 0 格式的结果,返回Table格式。

4.1 代码执行

实例:

XMLparser

perl  ../bin/blastXmlParser.pl   Human_kinase-rna-blastn-m5.xml   Human_kinase-rna-blastn-m5.tbl

perl  ../bin/blastXmlParser.pl   Human_kinase-rna-blastx-m5.xml   Human_kinase-rna-blastx-m5.tbl

perl  ../bin/blastXmlParser.pl  Human_kinase_protein-blastp-m5.xml   Human_kinase_protein-blastp-m5.tbl

TxtParser

perl   ../bin/blastTxtParser.pl   Human_kinase_protein-blastp-m0.txt   Human_kinase_protein-blastp-m0.tbl

perl   ../bin/blastTxtParser.pl   Human_kinase-rna-blastx-m0.txt  Human_kinase-rna-blastx-m0.tbl

perl   ../bin/blastTxtParser.pl   Human_kinase-rna-blastn-m0.txt  Human_kinase-rna-blastn-m0.tbl

4.2 文件

代码Host在了GitCafe, 所有示例文件都在百度云网盘 地址:

http://pan.baidu.com/share/link?shareid=2240069514&uk=1293299598

压缩包内容:

Bin目录为可执行程序: 包括了blastpTxtParser.pl  blastXmlParser.pl 以及blastn  blastp  blastx  makeblastdb

ccds 目录为 ccds 的human 蛋白质和核酸序列数据

example 目录 为实例文件:Human_kinase_rna-100.fasta 和 Human_kinase_protein-100.fasta

解压后可以先进入ccds 目录 进行makeblastdb, 然后可以进入example 目录进行运行BLAST程序以及提取程序。

8 comments to BLASTparser : 转换BLAST -m5/ -m0 格式成Table格式

Leave a Reply

To create code blocks or other preformatted text, indent by four spaces:

    This will be displayed in a monospaced font. The first four 
    spaces will be stripped off, but all other whitespace
    will be preserved.
    
    Markdown is turned off in code blocks:
     [This is not a link](http://example.com)

To create not a block, but an inline code span, use backticks:

Here is some inline `code`.

For more help see http://daringfireball.net/projects/markdown/syntax

You can use these HTML tags

<a href="" title=""> <abbr title=""> <acronym title=""> <b> <blockquote cite=""> <cite> <code> <del datetime=""> <em> <i> <q cite=""> <s> <strike> <strong>