#!/usr/bin/perl
# An example of how to use the different alignment tools in bioperl
# to align some sequences
# All these methods except Bio::Tools::pSW will work for DNA sequence
# (need to use a different matrix however)
use Bio::Factory::EMBOSS;
use Bio::SeqIO;
use Bio::AlignIO;
use Bio::Tools::pSW;
use Bio::PrimarySeq;
use Bio::Tools::Run::Alignment::Clustalw;
use Bio::Tools::Run::Alignment::TCoffee;
use Bio::Tools::Run::StandAloneBlast;
use strict;
# build the sequences since EMBOSS expects seqs to be in files
my $seq = new Bio::PrimarySeq(-seq =>
'MAVNPELAPFTLSRGIPSFDDQALSTIIQLQDCIQQAIQQLNYSTAEFLAELLYAECSILDKSSVYWSDAVYLYALSLFLNKSYHTAFQISKEFKEYHLGIAYIFGRCALQLSQGVNEAILTLLSIINVFSSNSSNTRINMVLNSNLVHIPDLATLNCLLGNLYMKLDHSKEGAFYHSEALAINPYLWESYEAICKMRATVDLKRVFFDIAGKKSNSHNNNAASSFPSTSLSHFEPRSQPSLYSKTNKNGNNNINNNVNTLFQSSNSPPSTSASSFSSIQHFSRSQQQQANTSIRTCQNKNTQTPKNPAINSKTSSALPNNISMNLVSPSSKQPTISSLAKVYNRNKLLTTPPSKLLNNDRNHQNNNNNNNNNNNNNNNNNNNNNNNNIINKTTFKTPRNLYSSTGRLTTSKKNPRSLIISNSILTSDYQITLPEIMYNFALILRSSSQYNSFKAIRLFESQIPSHIKDTMPWCLVQLGKLHFEIINYDMSLKYFNRLKDLQPARVKDMEIFSTLLWHLHDKVKSSNLANGLMDTMPNKPETWCCIGNLLSLQKDHDAAIKAFEKATQLDPNFAYAYTLQGHEHSSNDSSDSAKTCYRKALACDPQHYNAYYGLGTSAMKLGQYEEALLYFEKARSINPVNVVLICCCGGSLEKLGYKEKALQYYELACHLQPTSSLSKYKMGQLLYSMTRYNVALQTFEELVKLVPDDATAHYLLGQTYRIVGRKKDAIKELTVAMNLDPKGNQVIIDELQKCHMQE',
-id => 'seq1'
);
my $seq2 = new Bio::PrimarySeq( -seq =>
'CLIFXRLLLIQMIHPQARRAFTFLQQQEPYRIQSMEQLSTLLWHLADLPALSHLSQSLISISRSSPQAWIAVGNCFSLQKDHDEAMRCFRRATQVDEGCAYAWTLCGYEAVEMEEYERAMAFYRTAIRTDARHYNAWYVLFFFFFFFFVPGDIDSXPKKGMEWGXFISKRIDRGMRSIILKEPSKSIQLIPFFYVALVWXVGVSSYPLETMTNIDFPKKKKALEKSNDVVQALHFYERASKYAPTSAMVQFKRIRALVALQRYDEAISALVPLTHSAPDEANVFFLLGKCLLKKERRQEATMAFTNARELEPK',
-id => 'seq2');
my $out = new Bio::SeqIO(-format => 'fasta',
-file => ">seq1.fa");
$out->write_seq($seq);
$out->close();
$out = new Bio::SeqIO(-format => 'fasta',
-file => ">seq2.fa");
$out->write_seq($seq2);
$out->close();
my $embossfactory = Bio::Factory::EMBOSS->new();
my @alignprogs = qw(water needle stretcher matcher);
my $alignout = new Bio::AlignIO(-format => 'msf');
foreach my $prog ( @alignprogs ) {
my $alignfactory = $embossfactory->program('water');
$alignfactory->run({ '-sequencea' => 'seq1.fa',
'-seqall' => 'seq2.fa',
'-gapext' => 2.0,
'-datafile' => 'EBLOSUM62',
'-gapopen' => 14.0,
'-outfile' => "seq1_vs_seq2.$prog"});
my $alnin = new Bio::AlignIO(-format => 'emboss',
-file => "seq1_vs_seq2.$prog");
my $aln = $alnin->next_aln();
$alignout->write_aln($aln);
}
# this should produce the same alignment as 'water'
my $factory = new Bio::Tools::pSW(-matrix=> 'blosum62.bla',
-gap => 14,
-ext => 2);
my $aln = $factory->pairwise_alignment($seq,$seq2);
$alignout->write_aln($aln);
$factory = new Bio::Tools::Run::Alignment::Clustalw('ktuple' => 2,
'matrix' => 'BLOSUM');
$aln = $factory->align([$seq,$seq2]);
$alignout->write_aln($aln);
$factory = new Bio::Tools::Run::Alignment::TCoffee('ktuple' => 2,
'matrix' => 'BLOSUM');
$aln = $factory->align([$seq,$seq2]);
$alignout->write_aln($aln);
$factory = new Bio::Tools::Run::StandAloneBlast();
$aln = $factory->bl2seq($seq,$seq2);
# this actually returns a Bio::Tools::BPbl2seq object
# it can be transformed to a SimpleAlign object see
# the code in Bio::AlignIO::bl2seq
# A transformer object will be written at some point