[Bioperl-l] ClustalW alignment and Bioperl::Run module

Kevin Brown Kevin.M.Brown at asu.edu
Thu Sep 28 17:48:15 EDT 2006


I've gotten a very simple script to run using bioperl that creates an
alignment using clustalw of two sequences.  I see that clustal outputs
to stdout information like the score, but I don't see any way to store
that or retrieve that from the alignment object that is returned (unless
I'm just blind).  What follows is my very basic script which used code
found in the Wiki.

print $aln->score() spits out an error about using an uninitialized
value.


#!/usr/bin/perl -w

use strict;
use Bio::SeqIO;
use Bio::Perl;
use Bio::AlignIO;
use Getopt::Long qw(:config no_ignore_case bundling pass_through);
use POSIX;
use Bio::Tools::Run::Alignment::Clustalw;

my $fileName   = "";         # filename(s) to be parsed for information
my $output_dir = "";
my $format     = 'fasta';    # default format for SeqIO module

GetOptions(
                   'file=s'   => \$fileName,
                   'output=s' => \$output_dir,
                  );

# Parse the input file for the needed information
# SeqIO supports several normal formats including <tab>, <fasta> and
<excel>

my @files = split(/\|/, $fileName);
my @seq_array;

my $stream_out =
  Bio::AlignIO->new(-file => '>test.msf', -format => 'msf', -flush =>
0);

foreach my $fileName (@files)
{
        my $file = Bio::SeqIO->new(-format => $format, -file =>
$fileName);
        my $seq;
        while ($seq = $file->next_seq())
        {
                push(@seq_array, $seq);
        }
}

my @params  = ('ktuple' => 2, 'matrix' => 'BLOSUM');
my $factory = Bio::Tools::Run::Alignment::Clustalw->new(@params);
my $ktuple  = 3;
$factory->ktuple($ktuple);    # change the parameter before executing
    # where @seq_array is an array of {{PM|Bio::Seq}} objects

open my $out, ">seq.txt";

for (my $i = 1 ; $i <= $#seq_array ; $i++)
{
        my @seq = ($seq_array[0], $seq_array[$i]);
        my $aln = $factory->align(\@seq);
        $stream_out->write_aln($aln);
        print $aln->score;
        for my $seq ($aln->each_seq) {
                print $out $seq->display_id() ."\t". $seq->seq()."\n";
        }
}



More information about the Bioperl-l mailing list