← Back to branch summary

~ubuntu-branches/ubuntu/raring/bioperl/raring

« back to all changes in this revision

Viewing changes to Bio/Tools/dpAlign.pm

  • Committer: Bazaar Package Importer
  • Author(s): Charles Plessy
  • Date: 2008-03-18 14:44:57 UTC
  • mfrom: (4 hardy)
  • mto: This revision was merged to the branch mainline in revision 6.
  • Revision ID: james.westby@ubuntu.com-20080318144457-1jjoztrvqwf0gruk
http://bugs.debian.org/448890
* debian/control:
  - Removed MIA Matt Hope (dopey) from the Uploaders field.
    Thank you for your work, Matt. I hope you are doing well.
  - Downgraded some recommended package to the 'Suggests' priority,
    according to the following discussion on Upstream's mail list.
    http://bioperl.org/pipermail/bioperl-l/2008-March/027379.html
    (Closes: #448890)
* debian/copyright converted to machine-readable format.

Show diffs side-by-side

added added

removed removed

Lines of Context:
Bio/Tools/dpAlign.pm
1
 
## $Id: dpAlign.pm,v 1.2 2003/06/10 18:18:05 jason Exp $
 
1
## $Id: dpAlign.pm,v 1.18.4.2 2006/11/17 09:32:42 sendu Exp $
2
2
 
3
3
# BioPerl module for Bio::Tools::dpAlign
4
4
#
16
16
 
17
17
=head1 SYNOPSIS
18
18
 
19
 
        use Bio::Tools::dpAlign;
20
 
        use Bio::SeqIO;
21
 
        use Bio::SimpleAlign;
22
 
        use Bio::AlignIO;
23
 
 
24
 
        $seq1 = new Bio::SeqIO(-file => $ARGV[0], -format => 'fasta');
25
 
        $seq2 = new Bio::SeqIO(-file => $ARGV[1], -format => 'fasta');
26
 
 
27
 
        # create a dpAlign object
28
 
        $factory = new dpAlign(-match => 3,
29
 
                           -mismatch => -1,
30
 
                           -gap => 3,
31
 
                           -ext => 1,
32
 
                           -alg => Bio::Tools::dpAlign::DPALIGN_LOCAL_MILLER_MYERS);
33
 
 
34
 
        # actually do the alignment
35
 
        $out = $factory->pairwise_alignment($seq1->next_seq, $seq2->next_seq);
36
 
        $alnout = new Bio::AlignIO(-format => 'pfam', -fh => \*STDOUT);
37
 
        $alnout->write_aln($out);
38
 
 
39
 
        # To do protein alignment, set the sequence type to protein
40
 
        # currently all protein alignments are using BLOSUM62 matrix
41
 
        # the gap opening cost is 10 and gap extension is 2. These
42
 
        # values are from ssearch. They won't be changed even though
43
 
        # you set other values for now. Also DPALIGN_LOCAL_GREEN is not
44
 
        # supported for protein in this version.
45
 
        $seq1->alphabet('protein');
46
 
        $seq2->alphabet('protein');
47
 
        $out = $factory->pairwise_alignment($seq1->next_seq, $seq2->next_seq);
48
 
        $alnout->write_aln($out);
49
 
 
50
 
        # use the factory to make some output
51
 
 
52
 
        $factory->align_and_show($seq1, $seq2, STDOUT);
 
19
  use Bio::Tools::dpAlign;
 
20
  use Bio::SeqIO;
 
21
  use Bio::SimpleAlign;
 
22
  use Bio::AlignIO;
 
23
  use Bio::Matrix::IO;
 
24
 
 
25
  $seq1 = new Bio::SeqIO(-file => $ARGV[0], -format => 'fasta');
 
26
  $seq2 = new Bio::SeqIO(-file => $ARGV[1], -format => 'fasta');
 
27
 
 
28
  # create a dpAlign object
 
29
  # to do global alignment, specify DPALIGN_GLOBAL_MILLER_MYERS
 
30
  # to do ends-free alignment, specify DPALIGN_ENDSFREE_MILLER_MYERS
 
31
  $factory = new dpAlign(-match => 3,
 
32
                     -mismatch => -1,
 
33
                     -gap => 3,
 
34
                     -ext => 1,
 
35
                     -alg => Bio::Tools::dpAlign::DPALIGN_LOCAL_MILLER_MYERS);
 
36
 
 
37
  # actually do the alignment
 
38
  $out = $factory->pairwise_alignment($seq1->next_seq, $seq2->next_seq);
 
39
  $alnout = new Bio::AlignIO(-format => 'pfam', -fh => \*STDOUT);
 
40
  $alnout->write_aln($out);
 
41
 
 
42
  # To do protein alignment, set the sequence type to protein
 
43
  # By default all protein alignments are using BLOSUM62 matrix
 
44
  # the gap opening cost is 7 and gap extension is 1. These
 
45
  # values are from ssearch. To use your own custom substitution 
 
46
  # matrix, you can create a Bio::Matrix::MatrixI object.
 
47
 
 
48
  $parser = new Bio::Matrix::IO(-format => 'scoring', -file => 'blosum50.mat');
 
49
  $matrix = $parser->next_matrix;
 
50
  $factory = new Bio::Tools::dpAlign(-matrix => $matrix, -alg => Bio::Tools::dpAlign::DPALIGN_LOCAL_MILLERMYERS);
 
51
  $seq1->alphabet('protein');
 
52
  $seq2->alphabet('protein');
 
53
  $out = $factory->pairwise_alignment($seq1->next_seq, $seq2->next_seq);
 
54
  $alnout->write_aln($out);
 
55
 
 
56
  # use the factory to make some output
 
57
 
 
58
  $factory->align_and_show($seq1, $seq2, STDOUT);
 
59
 
 
60
  # use Phil Green's algorithm to calculate the optimal local
 
61
  # alignment score between two sequences quickly. It is very
 
62
  # useful when you are searching a query sequence in a database
 
63
  # of sequences. Since finding a alignment is more costly 
 
64
  # than just calculating scores, you can save time if you only 
 
65
  # align sequences that have a high alignment score.
 
66
 
 
67
  # To use this feature, first you call the sequence_profile function
 
68
  # to obtain the profile of the query sequence.
 
69
  $profile = $factory->sequence_profile($query);
 
70
 
 
71
  %scores = ();
 
72
  # Then use a loop to run a database of sequences against the
 
73
  # profile to obtain a table of alignment scores
 
74
  $dbseq = Bio::SeqIO(-file => 'dbseq.fa', -format => 'fasta');
 
75
  while (defined($seq = $dbseq->next_seq)) {
 
76
      $scores{$seq->id} = $factory->pairwise_alignment_score($profile, $seq);
 
77
  }
53
78
 
54
79
=head1 DESCRIPTION
55
80
 
56
 
        Dynamic Programming approach is considered to be the most
57
 
        sensitive way to align two biological sequences. There are
58
 
        currently three major types of dynamic programming algorithms:
59
 
        Global Alignment, Local Alignment and Ends-free Alignment.
60
 
 
61
 
        Global Alignment compares two sequences in their entirety.
62
 
        By inserting gaps in the two sequences, it aligns two
63
 
        sequences to minimize the edit distance as defined by the
64
 
        gap cost function and the substitution matrix. Global Alignment
65
 
        is generally applied to two sequences that are very similar
66
 
        in length and content.
67
 
 
68
 
        Local Alignment instead attempts to find out the subsequences
69
 
        that has the minimal edit distance among all possible subsequences.
70
 
        It is good for sequences that has a stretch of subsequences
71
 
        that are similar to each other.
72
 
 
73
 
        Ends-free Alignment is a special case of Global Alignment. There
74
 
        are no gap penalty imposed for the gaps that extended from
75
 
        the end points of two sequences. Therefore it will be a good
76
 
        application when you think one sequence is contained by the
77
 
        other or when you think two sequences overlap each other.
78
 
 
79
 
        Dynamic Programming was first introduced by Needleman-Wunsch (1970)
80
 
        to globally align two sequences. The idea of local alignment
81
 
        was later introduced by Smith-Waterman (1981). Gotoh (1982)
82
 
        improved both algorithms by introducing auxillary arrays that
83
 
        reduced the time complexity of the algorithms to O(m*n).
84
 
        Miller-Myers (1988) exploits the divide-and-conquer idea
85
 
        introduced by Hirschberg (1975) to solve the affine gap cost
86
 
        dynamic programming using only linear space. At the time of
87
 
        this writing, it is accepted that Miller-Myers is the fastest
88
 
        single CPU implementation and using the least memory that is
89
 
        truly equivalent to original algorithm introduced by
90
 
        Needleman-Wunsch. According to Aaron Mackey, Phil Green's SWAT
91
 
        implemention introduced a heuristic that does not consider
92
 
        paths throught the matrix where the score would be less than
93
 
        the gap opening penalty, yielding a 1.5-2X speedup on most
94
 
        comparisons. to skip the calculation of some cells. However,
95
 
        his approach is only good for calculating the minimum edit
96
 
        distance and find out the corresponding subsequences (aka
97
 
        search phase). Bill Pearson's popular dynamic programming
98
 
        alignment program SSEARCH uses Phil Green's algorithm to
99
 
        find the subsequences and then Miller-Myers's algorithm to
100
 
        find the actual alignment. (aka alignment phase)
101
 
 
102
 
        The current implementation supports local alignment of
103
 
        either DNA sequences or protein sequences. It allows you
104
 
        to specify either the Phil Green (DPALIGN_LOCAL_GREEN)
105
 
        or Miller-Myers (DPALIGN_LOCAL_MILLER_MYERS). For DNA
106
 
        alignment, you can specify the scores for match, mismatch,
107
 
        gap opening cost and gap extension cost. For protein
108
 
        alignment, it is using BLOSUM62 by default. Currently the
109
 
        substitution matrix is not configurable.
 
81
Dynamic Programming approach is considered to be the most sensitive
 
82
way to align two biological sequences. There are currently three major
 
83
types of dynamic programming algorithms: Global Alignment, Local
 
84
Alignment and Ends-free Alignment.
 
85
 
 
86
Global Alignment compares two sequences in their entirety.  By
 
87
inserting gaps in the two sequences, it aligns two sequences to
 
88
minimize the edit distance as defined by the gap cost function and the
 
89
substitution matrix. Global Alignment is generally applied to two
 
90
sequences that are very similar in length and content.
 
91
 
 
92
Local Alignment instead attempts to find out the subsequences that has
 
93
the minimal edit distance among all possible subsequences.  It is good
 
94
for sequences that has a stretch of subsequences that are similar to
 
95
each other.
 
96
 
 
97
Ends-free Alignment is a special case of Global Alignment. There are
 
98
no gap penalty imposed for the gaps that extended from the end points
 
99
of two sequences. Therefore it will be a good application when you
 
100
think one sequence is contained by the other or when you think two
 
101
sequences overlap each other.
 
102
 
 
103
Dynamic Programming was first introduced by Needleman-Wunsch (1970) to
 
104
globally align two sequences. The idea of local alignment was later
 
105
introduced by Smith-Waterman (1981). Gotoh (1982) improved both
 
106
algorithms by introducing auxillary arrays that reduced the time
 
107
complexity of the algorithms to O(m*n).  Miller-Myers (1988) exploits
 
108
the divide-and-conquer idea introduced by Hirschberg (1975) to solve
 
109
the affine gap cost dynamic programming using only linear space. At
 
110
the time of this writing, it is accepted that Miller-Myers is the
 
111
fastest single CPU implementation and using the least memory that is
 
112
truly equivalent to original algorithm introduced by
 
113
Needleman-Wunsch. According to Aaron Mackey, Phil Green's SWAT
 
114
implemention introduced a heuristic that does not consider paths
 
115
throught the matrix where the score would be less than the gap opening
 
116
penalty, yielding a 1.5-2X speedup on most comparisons. to skip the
 
117
calculation of some cells. However, his approach is only good for
 
118
calculating the minimum edit distance and find out the corresponding
 
119
subsequences (aka search phase). Bill Pearson's popular dynamic
 
120
programming alignment program SSEARCH uses Phil Green's algorithm to
 
121
find the subsequences and then Miller-Myers's algorithm to find the
 
122
actual alignment. (aka alignment phase)
 
123
 
 
124
The current implementation supports local alignment of either DNA
 
125
sequences or protein sequences. It allows you to specify either the
 
126
Miller-Myers Global Alignment (DPALIGN_GLOBAL_MILLER_MYERS) or
 
127
Miller-Myers Local Alignment (DPALIGN_LOCAL_MILLER_MYERS). For DNA
 
128
alignment, you can specify the scores for match, mismatch, gap opening
 
129
cost and gap extension cost. For protein alignment, it is using
 
130
BLOSUM62 by default. Currently the substitution matrix is not
 
131
configurable.
110
132
 
111
133
=head1 DEPENDENCIES
112
134
 
113
 
        This package comes with the main bioperl distribution. You
114
 
        also need to install the lastest bioperl-ext package which
115
 
        contains the XS code that implements the algorithms. This 
116
 
        package won't work if you haven't compiled the bioperl-ext
117
 
        package.
 
135
This package comes with the main bioperl distribution. You also need
 
136
to install the lastest bioperl-ext package which contains the XS code
 
137
that implements the algorithms. This package won't work if you haven't
 
138
compiled the bioperl-ext package.
118
139
 
119
140
=head1 TO-DO
120
141
 
121
 
        1) Allow custom substitution matrix.
122
 
 
123
 
        2) Support Global Alignment.
124
 
 
125
 
        3) Support Ends-free Alignment.
126
 
 
127
 
        4) Include a score only calculation based on Phil Green's
128
 
        algorithm. The code will be borrowed from do_work in
129
 
        ssearch.
130
 
 
131
 
        5) Support IUPAC code for DNA sequence
 
142
 
 
143
=over 3
 
144
 
 
145
=item 1.
 
146
 
 
147
Support IUPAC code for DNA sequence
 
148
 
 
149
=item 2.
 
150
 
 
151
Allow custom substitution matrix for DNA. Note that for proteins, you
 
152
can now use your own subsitution matirx.
 
153
 
 
154
=back
132
155
 
133
156
=head1 FEEDBACK
134
157
 
138
161
Bioperl modules.  Send your comments and suggestions preferably to one
139
162
of the Bioperl mailing lists.  Your participation is much appreciated.
140
163
 
141
 
    bioperl-l@bioperl.org              - General discussion
142
 
    http://bioperl.org/MailList.shtml  - About the mailing lists
 
164
  bioperl-l@bioperl.org                  - General discussion
 
165
  http://bioperl.org/wiki/Mailing_lists  - About the mailing lists
143
166
 
144
167
=head2 Reporting Bugs
145
168
 
146
169
Report bugs to the Bioperl bug tracking system to help us keep track
147
 
the bugs and their resolution. Bug reports can be submitted via email
148
 
or the web:
 
170
the bugs and their resolution. Bug reports can be submitted via the
 
171
web:
149
172
 
150
 
    bioperl-bugs@bio.perl.org
151
 
    http://bugzilla.bioperl.org/
 
173
  http://bugzilla.open-bio.org/
152
174
 
153
175
=head1 AUTHOR
154
176
 
164
186
 
165
187
package Bio::Tools::dpAlign;
166
188
 
167
 
use Bio::Tools::AlignFactory;
168
189
use Bio::SimpleAlign;
169
190
 
170
 
@ISA = qw(Bio::Tools::AlignFactory);
171
 
 
172
 
$VERSION = '0.50';
 
191
use base qw(Bio::Tools::AlignFactory);
173
192
 
174
193
# Gotoh algorithm as defined in J. Mol. Biol. (1982) 162, 705-708
175
194
# use constant DSW_GOTOH => 1;
178
197
# This algorithm is used in both the search phase and the
179
198
# alignment phase.
180
199
use constant DPALIGN_LOCAL_MILLER_MYERS => 1;
 
200
use constant DPALIGN_GLOBAL_MILLER_MYERS => 2;
 
201
use constant DPALIGN_ENDSFREE_MILLER_MYERS => 3;
181
202
# my toy algorithm that tries to do SW as fast as possible
182
203
# use constant DSW_FSW => 3; 
183
204
# Phil Green's approximation to Smith-Waterman. It avoid calculations
185
206
# This is the algorithm used by ssearch. Phil Green's algorithm is
186
207
# used in the search phase while Miller-Myers algorithm is used in
187
208
# the alignment phase
188
 
use constant DPALIGN_LOCAL_GREEN => 2; 
 
209
#use constant DPALIGN_LOCAL_GREEN => 2; 
189
210
 
190
211
BEGIN {
191
212
    eval {
198
219
}
199
220
 
200
221
sub new {
201
 
    my ($class, @args) = @_;
202
 
 
203
 
    my $self = $class->SUPER::new(@args);
204
 
 
205
 
    my ($match, $mismatch, $gap, $ext, $alg) = $self->_rearrange([qw(MATCH
206
 
                                                               MISMATCH
207
 
                                                               GAP
208
 
                                                               EXT
209
 
                                                               ALG      
 
222
   my ($class, @args) = @_;
 
223
 
 
224
   my $self = $class->SUPER::new(@args);
 
225
 
 
226
   my ($match, $mismatch, $gap, $ext, $alg, $matrix) = $self->_rearrange([qw(MATCH
 
227
                                                                MISMATCH
 
228
                                                                GAP
 
229
                                                                EXT
 
230
                                                                ALG
 
231
                                                                MATRIX  
210
232
                                                        )], @args);
211
233
 
212
 
    $self->match(3) unless defined $match;
213
 
    $self->mismatch(-1) unless defined $mismatch;
214
 
    $self->gap(3) unless defined $gap;
215
 
    $self->ext(1) unless defined $ext;
216
 
    $self->alg(DPALIGN_LOCAL_MILLER_MYERS) unless defined $alg;
 
234
   $self->match(3) unless defined $match;
 
235
   $self->mismatch(-1) unless defined $mismatch;
 
236
   $self->gap(3) unless defined $gap;
 
237
   $self->ext(1) unless defined $ext;
 
238
   $self->alg(DPALIGN_LOCAL_MILLER_MYERS) unless defined $alg;
217
239
 
218
 
    if (defined $match) {
 
240
   if (defined $match) {
219
241
        if ($match =~ /^\d+$/) {
220
242
            $self->match($match);
221
243
        }
252
274
    }
253
275
 
254
276
    if (defined $alg) {
255
 
        if ($alg == DPALIGN_LOCAL_MILLER_MYERS or $alg == DPALIGN_LOCAL_GREEN) {
 
277
        if ($alg == DPALIGN_LOCAL_MILLER_MYERS or $alg == DPALIGN_GLOBAL_MILLER_MYERS or $alg == DPALIGN_ENDSFREE_MILLER_MYERS) {
256
278
            $self->alg($alg);
257
279
        }
258
280
        else {
259
 
            $self->throw("Algorithm must be either 1 or 2");
 
281
            $self->throw("Algorithm must be either 1, 2 or 3");
260
282
        }
261
283
    }
 
284
 
 
285
    if (defined $matrix and $matrix->isa('Bio::Matrix::MatrixI')) {
 
286
        $self->{'matrix'} = Bio::Ext::Align::ScoringMatrix->new(join("", $matrix->row_names), $self->gap, $self->ext);
 
287
        foreach $rowname ($matrix->row_names) {
 
288
            foreach $colname ($matrix->column_names) {
 
289
                Bio::Ext::Align::ScoringMatrix->set_entry($self->{'matrix'}, $rowname, $colname, $matrix->entry($rowname, $colname));
 
290
            }
 
291
        }
 
292
    }
 
293
    else {
 
294
        $self->{'matrix'} = 0;
 
295
    }
 
296
 
262
297
    return $self;
263
298
}
264
299
 
 
300
=head2 sequence_profile
 
301
 
 
302
 Title   : sequence_profile
 
303
 Usage   : $prof = $factory->sequence_profile($seq1)
 
304
 Function: Makes a dpAlign_SequenceProfile object from one sequence
 
305
 Returns : A dpAlign_SequenceProfile object
 
306
 Args    : The lone argument is a Bio::PrimarySeqI that we want to 
 
307
           build a profile for. Usually, this would be the Query sequence
 
308
 
 
309
=cut
 
310
 
 
311
sub sequence_profile {
 
312
    my ($self, $seq1) = @_;
 
313
 
 
314
    if( ! defined $seq1 || ! $seq1->isa('Bio::PrimarySeqI')) {
 
315
        $self->warn("Cannot call sequence_profilewithout specifing one sequence (Bio::PrimarySeqI object)");
 
316
        return;
 
317
    }
 
318
 
 
319
    # fix Jitterbug #1044
 
320
    if( $seq1->length() < 2) {
 
321
        $self->warn("cannot create sequence profile with length less than 2");
 
322
        return;
 
323
    }
 
324
    # create engine objects
 
325
    $seq1->display_id('seq1') unless ( defined $seq1->id() );
 
326
 
 
327
    if ($seq1->alphabet eq 'dna') {
 
328
        return Bio::Ext::Align::SequenceProfile->dna_new($seq1->seq, $self->{'match'}, $self->{'mismatch'}, $self->{'gap'}, $self->{'ext'});
 
329
    }
 
330
    elsif ($seq1->alphabet eq 'protein') {
 
331
        return Bio::Ext::Align::SequenceProfile->protein_new($seq1->seq, $self->{'matrix'}); 
 
332
    }
 
333
    else {
 
334
        croak("There is currently no support for the types of sequences you want to align!\n");
 
335
        return;
 
336
    }
 
337
}
 
338
 
 
339
=head2 pairwise_alignment_score
 
340
 
 
341
 Title   : pairwise_alignment_score
 
342
 Usage   : $score = $factory->pairwise_alignment_score($prof,$seq2)
 
343
 Function: Makes a SimpleAlign object from two sequences
 
344
 Returns : An integer that is the score of the optimal alignment.
 
345
 Args    : The first argument is the sequence profile obtained from a
 
346
           call to the sequence_profile function. The second argument 
 
347
           is a Bio::PrimarySeqI object to be aligned. The second argument
 
348
           is usually a sequence in the database sequence. Note
 
349
           that this function only uses Phil Green's algorithm and 
 
350
           therefore theoretically may not always give you the optimal
 
351
           score.
 
352
 
 
353
=cut
 
354
 
 
355
sub pairwise_alignment_score {
 
356
    my ($self, $prof, $seq2) = @_;
 
357
 
 
358
    if( ! defined $prof || ! $prof->isa('Bio::Ext::Align::SequenceProfile') || 
 
359
        ! defined $seq2 || ! $seq2->isa('Bio::PrimarySeqI') ) {
 
360
        $self->warn("Cannot call pairwise_alignment_score without specifing 2 sequences (Bio::PrimarySeqI objects)");
 
361
        return;
 
362
    }
 
363
    # fix Jitterbug #1044
 
364
    if( $seq2->length() < 2) {
 
365
        $self->warn("cannot align sequences with length less than 2");
 
366
        return;
 
367
    }
 
368
    $self->set_memory_and_report();
 
369
    # create engine objects
 
370
    $seq2->display_id('seq2') unless ( defined $seq2->id() );
 
371
 
 
372
    if ($prof->alphabet eq 'dna' and $seq2->alphabet eq 'dna') {
 
373
        return Bio::Ext::Align::Score_DNA_Sequences($prof, $seq2->seq);
 
374
    }
 
375
    elsif ($prof->alphabet eq 'protein' and $seq2->alphabet eq 'protein') {
 
376
        return Bio::Ext::Align::Score_Protein_Sequences($prof, $seq2->seq);
 
377
    }
 
378
    else {
 
379
        croak("There is currently no support for the types of sequences you want to align!\n");
 
380
        return;
 
381
    }
 
382
}
 
383
 
265
384
=head2 pairwise_alignment
266
385
 
267
386
 Title   : pairwise_alignment
268
387
 Usage   : $aln = $factory->pairwise_alignment($seq1,$seq2)
269
388
 Function: Makes a SimpleAlign object from two sequences
270
 
 Returns : A SimpleAlign object
271
 
 Args    :
272
 
 
 
389
 Returns : A SimpleAlign object if there is an alignment with positive
 
390
           score. Otherwise, return undef.
 
391
 Args    : The first and second arguments are both Bio::PrimarySeqI
 
392
           objects that are to be aligned.
273
393
 
274
394
=cut
275
395
 
280
400
    if( ! defined $seq1 || ! $seq1->isa('Bio::PrimarySeqI') ||
281
401
        ! defined $seq2 || ! $seq2->isa('Bio::PrimarySeqI') ) {
282
402
        $self->warn("Cannot call pairwise_alignment without specifing 2 sequences (Bio::PrimarySeqI objects)");
283
 
        return undef;
 
403
        return;
284
404
    }
285
405
    # fix Jitterbug #1044
286
406
    if( $seq1->length() < 2 ||
287
407
        $seq2->length() < 2 ) {
288
408
        $self->warn("cannot align sequences with length less than 2");
289
 
        return undef;
 
409
        return;
290
410
    }
291
411
    $self->set_memory_and_report();
292
412
    # create engine objects
297
417
        $aln = Bio::Ext::Align::Align_DNA_Sequences($seq1->seq, $seq2->seq, $self->{'match'}, $self->{'mismatch'}, $self->{'gap'}, $self->{'ext'}, $self->{'alg'});
298
418
    }
299
419
    elsif ($seq1->alphabet eq 'protein' and $seq2->alphabet eq 'protein') {
300
 
        $aln = Bio::Ext::Align::Align_Protein_Sequences($seq1->seq, $seq2->seq, 'blosum62');
 
420
        $aln = Bio::Ext::Align::Align_Protein_Sequences($seq1->seq, $seq2->seq, $self->{'matrix'}, $self->{'alg'});
301
421
    }
302
422
    else {
303
423
        croak("There is currently no support for the types of sequences you want to align!\n");
 
424
        return;
 
425
    }
 
426
 
 
427
    if (not defined $aln or $aln == 0) {
 
428
        return;
304
429
    }
305
430
 
306
431
    $out = Bio::SimpleAlign->new();
314
439
                                         -start => $aln->start2,
315
440
                                         -end => $aln->end2,
316
441
                                         -id => $seq2->id));
317
 
 
 
442
    $out->score($aln->score);
318
443
    return $out;
319
444
}
320
445
 
335
460
    if( ! defined $seq1 || ! $seq1->isa('Bio::PrimarySeqI') ||
336
461
        ! defined $seq2 || ! $seq2->isa('Bio::PrimarySeqI') ) {
337
462
        $self->warn("Cannot call pairwise_alignment without specifing 2 sequences (Bio::PrimarySeqI objects)");
338
 
        return undef;
 
463
        return;
339
464
    }
340
465
    # fix Jitterbug #1044
341
466
    if( $seq1->length() < 2 ||
342
467
        $seq2->length() < 2 ) {
343
468
        $self->warn("cannot align sequences with length less than 2");
344
 
        return undef;
 
469
        return;
345
470
    }
346
471
    $self->set_memory_and_report();
347
472
    # create engine objects
352
477
        $aln = Bio::Ext::Align::Align_DNA_Sequences($seq1->seq, $seq2->seq, $self->{'match'}, $self->{'mismatch'}, $self->{'gap'}, $self->{'ext'}, $self->{'alg'});
353
478
    }
354
479
    elsif ($seq1->alphabet eq 'protein' and $seq2->alphabet eq 'protein') {
355
 
        $aln = Bio::Ext::Align::Align_Protein_Sequences($seq1->seq, $seq2->seq, 'blosum62');
 
480
        $aln = Bio::Ext::Align::Align_Protein_Sequences($seq1->seq, $seq2->seq, $self->{'matrix'}, $self->{'alg'});
356
481
    }
357
482
    else {
358
483
        croak("There is currently no support for the types of sequences you want to align!\n");
359
484
    }
360
485
 
361
 
    $out = Bio::SimpleAlign->new();
362
 
 
363
 
    $out->add_seq(Bio::LocatableSeq->new(-seq => $aln->aln1,
364
 
                                         -start => $aln->start1,
365
 
                                         -end => $aln->end1,
366
 
                                         -id => $seq1->id));
367
 
    
368
 
    $out->add_seq(Bio::LocatableSeq->new(-seq => $aln->aln2,
369
 
                                         -start => $aln->start2,
370
 
                                         -end => $aln->end2,
371
 
                                         -id => $seq2->id));
372
 
 
373
 
    my $t1 = &Bio::Ext::Align::new_Sequence_from_strings($seq1->id, $seq1->seq);
374
 
    my $t2 = &Bio::Ext::Align::new_Sequence_from_strings($seq2->id, $seq2->seq);
375
 
 
376
 
    &Bio::Ext::Align::write_pretty_seq_align($out,$t1,$t2,12,50,$fh);
 
486
    $out = Bio::Ext::Align::AlnBlock->new();
 
487
    my $s1 = Bio::Ext::Align::AlnSequence->new();
 
488
    my $s2 = Bio::Ext::Align::AlnSequence->new();
 
489
    my $a1 = $aln->aln1;
 
490
    my $a2 = $aln->aln2;
 
491
    my $first_col = undef;
 
492
    my $last_col = undef;
 
493
    my $col;
 
494
    my $alu1;
 
495
    my $alu2;
 
496
    my $g1 = 0;
 
497
    my $g2 = 0;
 
498
 
 
499
# construct AlnBlock
 
500
    for (my $i = 0; $i < length($a1); ++$i) {
 
501
        $col = Bio::Ext::Align::AlnColumn->new();
 
502
        $alu1 = Bio::Ext::Align::AlnUnit->new();
 
503
        $alu2 = Bio::Ext::Align::AlnUnit->new();
 
504
        $first_col = $col unless defined $first_col;
 
505
        Bio::Ext::Align::AlnColumn::set_next($last_col, $col) if defined $last_col;
 
506
        
 
507
        if (substr($a1, $i, 1) eq "-") {
 
508
            Bio::Ext::Align::AlnUnit::set_text_label($alu1, "INSERT");
 
509
            Bio::Ext::Align::AlnUnit::set_text_label($alu2, "SEQUENCE");
 
510
            ++$g1;
 
511
        }
 
512
        elsif (substr($a2, $i, 1) eq "-") {
 
513
            Bio::Ext::Align::AlnUnit::set_text_label($alu1, "SEQUENCE");
 
514
            Bio::Ext::Align::AlnUnit::set_text_label($alu2, "INSERT");
 
515
            ++$g2;
 
516
        }
 
517
        else {
 
518
            Bio::Ext::Align::AlnUnit::set_text_label($alu1, "SEQUENCE");
 
519
            Bio::Ext::Align::AlnUnit::set_text_label($alu2, "SEQUENCE");
 
520
        }
 
521
 
 
522
        Bio::Ext::Align::AlnUnit::set_start($alu1, $aln->start1+$i-$g1-2);
 
523
        Bio::Ext::Align::AlnUnit::set_end($alu1, $aln->start1+$i-$g1-2);
 
524
        Bio::Ext::Align::AlnUnit::set_start($alu2, $aln->start2+$i-$g2-2);
 
525
        Bio::Ext::Align::AlnUnit::set_end($alu2, $aln->start2+$i-$g2-2);
 
526
        Bio::Ext::Align::AlnColumn::add_alu($col, $alu1);
 
527
        Bio::Ext::Align::AlnColumn::add_alu($col, $alu2);
 
528
        $last_col = $col;
 
529
    }
 
530
    Bio::Ext::Align::AlnBlock::set_start($out, $first_col);
 
531
    $col = Bio::Ext::Align::AlnColumn->new();
 
532
    $alu1 = Bio::Ext::Align::AlnUnit->new();
 
533
    $alu2 = Bio::Ext::Align::AlnUnit->new();
 
534
    Bio::Ext::Align::AlnUnit::set_start($alu1, $aln->end1);
 
535
    Bio::Ext::Align::AlnUnit::set_end($alu1, $aln->end1);
 
536
    Bio::Ext::Align::AlnUnit::set_text_label($alu1, "END");
 
537
    Bio::Ext::Align::AlnUnit::set_start($alu2, $aln->end2);
 
538
    Bio::Ext::Align::AlnUnit::set_end($alu2, $aln->end2);
 
539
    Bio::Ext::Align::AlnUnit::set_text_label($alu2, "END");
 
540
    Bio::Ext::Align::AlnColumn::add_alu($col, $alu1);
 
541
    Bio::Ext::Align::AlnColumn::add_alu($col, $alu2);
 
542
    Bio::Ext::Align::AlnColumn::set_next($last_col, $col);
 
543
 
 
544
    &Bio::Ext::Align::write_pretty_str_align($out,$seq1->id,$seq1->seq,$seq2->id,$seq2->seq,12,50,$fh);
377
545
}
378
546
 
379
547
=head2 match
492
660
    my ($self,$val) = @_;
493
661
 
494
662
    if( defined $val ) {
495
 
        if( $val != DPALIGN_LOCAL_MILLER_MYERS and $val != DPALIGN_LOCAL_GREEN) {    
496
 
            $self->throw("Can't have an algorithm that is not 1 or 2");
 
663
        if( $val != DPALIGN_LOCAL_MILLER_MYERS and $val != DPALIGN_GLOBAL_MILLER_MYERS and $val != DPALIGN_ENDSFREE_MILLER_MYERS) {    
 
664
            $self->throw("Can't have an algorithm that is not 1, 2 or 3");
497
665
        }
498
666
        $self->{'alg'} = $val;
499
667
    }