[Bioperl-l] Alignments from psiblast

[Arne Elofsson]

Hi

I just made my first commitment to the bioperl CVS entry. I hope
everything was OK. I do not know if you want a record of what I did
for the mailing list or for anywhere else. However this is what I did

The goal was to be able to produce an alignment (a simple alignment
object) from a psiblast output. The only (easy) way to do this is to
run psiblast with the "-m 6" flag, therefore I had to make som other
changes to. However I think it all works

This is more details
1) A small fix of ID in Bio/SeqIO/fasta.pm

2) Changes in Tools/BPlite/Iteration.pm to handle a different psiblast
output format Some bugs still remains to handle it perfectly.

3) Changed one error message in BPlite.pm so that it is not identical
to one in Iteration.pm

4) Changed so that you can give the full path to the database in
Bio/Tools/Run/StandAloneBlast

Here follows a small program that uses the alignments.
-----
#!/usr/bin/perl

# This is an example how to use the Align module to achievean alignment
# from psiblast


use strict;

use lib '/home/ae/arne/splicing/';
use Bio::Tools::Blast;
use Bio::Tools::BPlite;
use Bio::Tools::BPpsilite;
use Bio::Tools::Run::StandAloneBlast;
use Bio::SeqIO;
use Bio::AlignIO;
use Bio::Seq;
use Bio::Root::IO;
use Bio::Index::Fasta;

#my $db='/home/ae/arne/structpredict/alignments/data/kind';
my $db='swissprot';
my $dir='/home/ae/arne/splicing/';
my $psidir=$dir."/psi/";
my $slxdir=$dir."/slx/";
my $out=$dir.'temp/psiblast.out.'.$$;
my $index=$dir.'temp/bptest.index.'.$$;
my ($seq);
foreach  (@ARGV){
  chomp;
  my $stream = Bio::SeqIO->new(-format => 'Fasta',-file => $_ );
  foreach  ( $seq = $stream->next_seq ) {
    my $id=$seq->display_id;
    my $psiout = $psidir.$id.".psi";
# We need to use the -m 6 flag
    my @psiparams = ('database' => $db , 'output' => $out, 'j' => 3, 'm' => 6,
		     'h' => 1.e-3 , 'F' => 'T' , 'Q' => $psiout ); # 'v' => 99999999,
    my $factory = Bio::Tools::Run::StandAloneBlast->new(@psiparams);
    my $report = $factory->blastpgp($seq);
    my $total_iterations = $report->number_of_iterations();
    my $last_iteration = $report->round($total_iterations);
    my $align=$last_iteration->Align;
    my $slxfile=$slxdir.$id.".slx";
    my $slx = Bio::AlignIO->new('-format' => 'selex','-file' => ">".$slxfile );
    $slx->write_aln($align);
  }
}
-------

yours

arne



-- 
------------------------------------------------------------------------
 Arne Elofsson     Stockholm Bioinformatics Center
 Net:      arne@sbc.su.se http://www.sbc.su.se/~arne/
 Tel:+46-8-161553   Stockholm Bioinformatics Center, Stockholm University
 Fax:+46-8-158057   10691 Stockholm, Sweden