[Bioinfo] NCBI BLAST command line

1. Download executables

ftp://ftp.ncbi.nlm.nih.gov/blast/executables/blast+/LATEST

2. Add the paths of executable and databse to .bash_profile in home directory

vi ~/.bash_profile

enter i to get into edit mode
copy and paste the text below

PATH=$PATH:$HOME/bin:/home/chlin/bin/ncbi-blast-2.2.29+/bin
BLASTDB=/home/chlin/blast/db

press esc then : then type wq to save and quit editing the file

3. (optional) Load the path setting immidiately

sourse .bash_profile

4. Make BLAST database

makeblastdb -in <input_fa> -dbtype <db_type> -out <db_name> -title <db_title>
Examples:
(for nucleotides)

makeblastdb -in hs_chr.fa -dbtype nucl -out hs_chr -title "Human Chromosome, Ref B37.1"

(for peptides)

makeblastdb -in refseq_protein.fa -dbtype prot -out refseq_protein -title "RefSeq Protein Database"

5. BLAST

blastn -query <query.fa> -db <db_name> -out <output_name>
Example:
blastn -query test.fa -db hs_chr -out hs_chr.out

options:
-evalue 1e-10 (E value threshold)
-num_threads 8 (Number of threads (CPUs) to use in blast search.)
-perc_identity 97 (% identity cutoff)
-num_descriptions 5 (hit # threshold)

6. Parse results

Perl module (http://www.bioperl.org/wiki/HOWTO:SearchIO)
Example:

use Bio::SearchIO;
my $in = new Bio::SearchIO( -format => 'blast', -file => "hs_chr.out" );
my $n=0;
while( my $result = $in->next_result ){
    my $query = $result->query_name;
    while ( my $hit = $result->next_hit ){
        while( my $hsp = $hit->next_hsp ){
            print "Query = ", $result->query_name, "\n",
            "Num_hits = ", $result->num_hits,"\n",
            "Hit = ", $hit->name,"\n",
            "Length = ", $hsp->length('total'), "\n",
            "Percent_id = ", $hsp->percent_identity, "\n",
            "Start = ", $hsp->start('hit'), "\n",
            "End = ", $hsp->end('hit'), "\n",
            "E-Value = ", $hsp->evalue,"\n\n";
        }
    }
}

Referces:

1. NCBI Manual (http://www.ncbi.nlm.nih.gov/books/NBK1763/)