Perl: make sequencing alignment using local blast

Abstract: Carry out sequence alignment and extract the results using the aligner blast2+ at the local computer

1. Install Blast

The installers and source code of Blast are available from NCBI FTP site:

#wget -ck ftp://ftp.ncbi.nlm.nih.gov/blast/executables/blast+/LATEST/ncbi-blast-2.3.0+-x64-linux.tar.gz

#sudo Tar -xvpf ncbi-blast-2.3.0+-x64-linux.tar.gz

#cd ncbi-blast-2.3.0+-x64-linux/bin

#sudo cp * /usr/bin/

Another quick way is:

#sudo apt-get install blast2

#Blastn -version

2. Test Blast

Create local sequence database:

#makeblastdb -in VirScan_v1.fa -dbtype nucl -parse_seqids -out ../blast/VirScan_v1

Test:

#blastn -query test.fa -db ../blast/VirScan_v1 -outfmt 11 -num_threads 8 -out test.bls

3. One sequence alignment using BioPerl

Here is the code:

#!/usr/bin/perl

use warnings;

#use strict;

use Data::Dumper;

use Bio::SimpleAlign;

use warnings;

use strict;

use Bio::Tools::Run::StandAloneBlast;

my $blast_obj=Bio::Tools::Run::StandAloneBlast->new(-program=>'blastn',

-database=>'/home/yuan/phip/blast/VirScan_v1',

-outfile=>'/home/yuan/phip/blast/test.bls');

#align one sequence

my $seq_obj=Bio::Seq->new(-id=>"test", -seq=>"TTTTATTTGGTTGT");

my $report_obj=$blast_obj->blastall($seq_obj);

while( my $result_obj = $report_obj->next_result ) {

while( my $hit= $result_obj->next_hit ) {

my $hit_name=$hit->name;

while( my $hsp = $hit->next_hsp ) {

printf("Hit:%s\t Length:%s\t Percent_id:%s\n", $hit_name,$hsp->length('total'), $hsp->percent_identity);

}

}

}

3. Multiple sequence alignment

The code:

#!/usr/bin/perl

use warnings;

use strict;

use Bio::Tools::Run::StandAloneBlast;

my $blast_obj=Bio::Tools::Run::StandAloneBlast->new(-program=>'blastn',

-database=>'/home/yuan/phip/blast/VirScan_v1', -expect => 0.01,

-outfile=>'/home/yuan/phip/blast/test.bls');

#input is fasta file

my $seqio=Bio::SeqIO->new(-file=>"/home/yuan/phip/ref_seq/test.fa", -format=>"Fasta");

while(my $seq_obj=$seqio->next_seq()){

my $report_obj=$blast_obj->blastall($seq_obj);

while( my $result_obj = $report_obj->next_result ) {

my $query_name=$result_obj->query_name;

my $query_acc=$result_obj->query_accession;

my $query_desc=$result_obj->query_description;

while( my $hit= $result_obj->next_hit ) {

my $hit_name=$hit->name;

my $hits_num=$hit->num_hsps;

printf("Query name: %s\t Query accession:%s\t Description: %s\t Hits_num:%s\n", $query_name, $query_acc, $query_desc, $hits_num);

while( my $hsp = $hit->next_hsp ) {

my $identity_percent= $hsp->percent_identity;

if ( $identity_percent > 90) {

printf("Hit name:%s\t Length:%s\t Percent_id:%s\t Score:%s\n", $hit_name,$hsp->length('total'), $identity_percent, $hsp->score);

}

}

}

}

}