Bio/Tools/dpAlign.pm

   1
   2 # BioPerl module for Bio::Tools::dpAlign
   3 #
   4 # Please direct questions and support issues to <bioperl-l@bioperl.org>
   5 #
   6 # Cared for by Yee Man Chan <ymc@yahoo.com>
   7 #
   8 # Copyright Yee Man Chan
   9 #
  10 # You may distribute this module under the same terms as perl itself
  11
  12 # POD documentation - main docs before the code
  13
  14 =head1 NAME
  15
  16 Bio::Tools::dpAlign - Perl extension to do pairwise dynamic programming sequence alignment
  17
  18 =head1 SYNOPSIS
  19
  20   use Bio::Tools::dpAlign;
  21   use Bio::SeqIO;
  22   use Bio::SimpleAlign;
  23   use Bio::AlignIO;
  24   use Bio::Matrix::IO;
  25
  26   $seq1 = Bio::SeqIO->new(-file => $ARGV[0], -format => 'fasta');
  27   $seq2 = Bio::SeqIO->new(-file => $ARGV[1], -format => 'fasta');
  28
  29   # create a dpAlign object
  30   # to do global alignment, specify DPALIGN_GLOBAL_MILLER_MYERS
  31   # to do ends-free alignment, specify DPALIGN_ENDSFREE_MILLER_MYERS
  32   $factory = new dpAlign(-match => 3,
  33                      -mismatch => -1,
  34                      -gap => 3,
  35                      -ext => 1,
  36                      -alg => Bio::Tools::dpAlign::DPALIGN_LOCAL_MILLER_MYERS);
  37
  38   # actually do the alignment
  39   $out = $factory->pairwise_alignment($seq1->next_seq, $seq2->next_seq);
  40   $alnout = Bio::AlignIO->new(-format => 'pfam', -fh => \*STDOUT);
  41   $alnout->write_aln($out);
  42
  43   # To do protein alignment, set the sequence type to protein
  44   # By default all protein alignments are using BLOSUM62 matrix
  45   # the gap opening cost is 7 and gap extension is 1. These
  46   # values are from ssearch. To use your own custom substitution
  47   # matrix, you can create a Bio::Matrix::MatrixI object.
  48
  49   $parser = Bio::Matrix::IO->new(-format => 'scoring', -file => 'blosum50.mat');
  50   $matrix = $parser->next_matrix;
  51   $factory = Bio::Tools::dpAlign->new(-matrix => $matrix, -alg => Bio::Tools::dpAlign::DPALIGN_LOCAL_MILLERMYERS);
  52   $seq1->alphabet('protein');
  53   $seq2->alphabet('protein');
  54   $out = $factory->pairwise_alignment($seq1->next_seq, $seq2->next_seq);
  55   $alnout->write_aln($out);
  56
  57   # use the factory to make some output
  58
  59   $factory->align_and_show($seq1, $seq2, STDOUT);
  60
  61   # use Phil Green's algorithm to calculate the optimal local
  62   # alignment score between two sequences quickly. It is very
  63   # useful when you are searching a query sequence in a database
  64   # of sequences. Since finding a alignment is more costly
  65   # than just calculating scores, you can save time if you only
  66   # align sequences that have a high alignment score.
  67
  68   # To use this feature, first you call the sequence_profile function
  69   # to obtain the profile of the query sequence.
  70   $profile = $factory->sequence_profile($query);
  71
  72   %scores = ();
  73   # Then use a loop to run a database of sequences against the
  74   # profile to obtain a table of alignment scores
  75   $dbseq = Bio::SeqIO(-file => 'dbseq.fa', -format => 'fasta');
  76   while (defined($seq = $dbseq->next_seq)) {
  77       $scores{$seq->id} = $factory->pairwise_alignment_score($profile, $seq);
  78   }
  79
  80 =head1 DESCRIPTION
  81
  82 Dynamic Programming approach is considered to be the most sensitive
  83 way to align two biological sequences. There are currently three major
  84 types of dynamic programming algorithms: Global Alignment, Local
  85 Alignment and Ends-free Alignment.
  86
  87 Global Alignment compares two sequences in their entirety.  By
  88 inserting gaps in the two sequences, it aligns two sequences to
  89 minimize the edit distance as defined by the gap cost function and the
  90 substitution matrix. Global Alignment is generally applied to two
  91 sequences that are very similar in length and content.
  92
  93 Local Alignment instead attempts to find out the subsequences that has
  94 the minimal edit distance among all possible subsequences.  It is good
  95 for sequences that has a stretch of subsequences that are similar to
  96 each other.
  97
  98 Ends-free Alignment is a special case of Global Alignment. There are
  99 no gap penalty imposed for the gaps that extended from the end points
 100 of two sequences. Therefore it will be a good application when you
 101 think one sequence is contained by the other or when you think two
 102 sequences overlap each other.
 103
 104 Dynamic Programming was first introduced by Needleman-Wunsch (1970) to
 105 globally align two sequences. The idea of local alignment was later
 106 introduced by Smith-Waterman (1981). Gotoh (1982) improved both
 107 algorithms by introducing auxillary arrays that reduced the time
 108 complexity of the algorithms to O(m*n).  Miller-Myers (1988) exploits
 109 the divide-and-conquer idea introduced by Hirschberg (1975) to solve
 110 the affine gap cost dynamic programming using only linear space. At
 111 the time of this writing, it is accepted that Miller-Myers is the
 112 fastest single CPU implementation and using the least memory that is
 113 truly equivalent to original algorithm introduced by
 114 Needleman-Wunsch. According to Aaron Mackey, Phil Green's SWAT
 115 implementation introduced a heuristic that does not consider paths
 116 through the matrix where the score would be less than the gap opening
 117 penalty, yielding a 1.5-2X speedup on most comparisons. to skip the
 118 calculation of some cells. However, his approach is only good for
 119 calculating the minimum edit distance and find out the corresponding
 120 subsequences (aka search phase). Bill Pearson's popular dynamic
 121 programming alignment program SSEARCH uses Phil Green's algorithm to
 122 find the subsequences and then Miller-Myers's algorithm to find the
 123 actual alignment. (aka alignment phase)
 124
 125 The current implementation supports local alignment of either DNA
 126 sequences or protein sequences. It allows you to specify either the
 127 Miller-Myers Global Alignment (DPALIGN_GLOBAL_MILLER_MYERS) or
 128 Miller-Myers Local Alignment (DPALIGN_LOCAL_MILLER_MYERS). For DNA
 129 alignment, you can specify the scores for match, mismatch, gap opening
 130 cost and gap extension cost. For protein alignment, it is using
 131 BLOSUM62 by default. Currently the substitution matrix is not
 132 configurable.
 133
 134 Note: If you supply LocatableSeq objects to pairwise_alignment,
 135 pairwise_alignment_score, align_and_show or sequence_profile and
 136 the sequence supplied contains gaps, these functions will treat
 137 these sequences as if they are without gaps.
 138
 139 =head1 DEPENDENCIES
 140
 141 This package comes with the main bioperl distribution. You also need
 142 to install the lastest bioperl-ext package which contains the XS code
 143 that implements the algorithms. This package won't work if you haven't
 144 compiled the bioperl-ext package.
 145
 146 =head1 TO-DO
 147
 148
 149 =over 3
 150
 151 =item 1.
 152
 153 Basic support for IUPAC code for DNA sequence is now implemented.
 154 X will mismatch any character. T will match U. For others, whenever
 155 there is a possibility for match, it is considered a full match, for
 156 example, W will match B.
 157
 158 =item 2.
 159
 160 Allow custom substitution matrix for DNA. Note that for proteins, you
 161 can now use your own subsitution matirx.
 162
 163 =back
 164
 165 =head1 FEEDBACK
 166
 167 =head2 Mailing Lists
 168
 169 User feedback is an integral part of the evolution of this and other
 170 Bioperl modules.  Send your comments and suggestions preferably to one
 171 of the Bioperl mailing lists.  Your participation is much appreciated.
 172
 173   bioperl-l@bioperl.org                  - General discussion
 174   http://bioperl.org/wiki/Mailing_lists  - About the mailing lists
 175
 176 =head2 Support
 177
 178 Please direct usage questions or support issues to the mailing list:
 179
 180 I<bioperl-l@bioperl.org>
 181
 182 rather than to the module maintainer directly. Many experienced and
 183 reponsive experts will be able look at the problem and quickly
 184 address it. Please include a thorough description of the problem
 185 with code and data examples if at all possible.
 186
 187 =head2 Reporting Bugs
 188
 189 Report bugs to the Bioperl bug tracking system to help us keep track
 190 the bugs and their resolution. Bug reports can be submitted via the
 191 web:
 192
 193   https://github.com/bioperl/bioperl-live/issues
 194
 195 =head1 AUTHOR
 196
 197         This implementation was written by Yee Man Chan (ymc@yahoo.com).
 198         Copyright (c) 2003 Yee Man Chan. All rights reserved. This program
 199         is free software; you can redistribute it and/or modify it under
 200         the same terms as Perl itself. Special thanks to Aaron Mackey
 201         and WIlliam Pearson for the helpful discussions. [The portion
 202         of code inside pgreen subdirectory was borrowed from ssearch. It
 203         should be distributed in the same terms as ssearch.]
 204
 205 =cut
 206
 207 package Bio::Tools::dpAlign;
 208
 209 use Bio::SimpleAlign;
 210
 211 use base qw(Bio::Tools::AlignFactory);
 212
 213 # Gotoh algorithm as defined in J. Mol. Biol. (1982) 162, 705-708
 214 # use constant DSW_GOTOH => 1;
 215 # Hirschberg's algorithm as defined in Myers & Miller in
 216 # CABIOS, Vol 4, No. 1, 1988, p 11-17
 217 # This algorithm is used in both the search phase and the
 218 # alignment phase.
 219 use constant DPALIGN_LOCAL_MILLER_MYERS => 1;
 220 use constant DPALIGN_GLOBAL_MILLER_MYERS => 2;
 221 use constant DPALIGN_ENDSFREE_MILLER_MYERS => 3;
 222 # my toy algorithm that tries to do SW as fast as possible
 223 # use constant DSW_FSW => 3;
 224 # Phil Green's approximation to Smith-Waterman. It avoid calculations
 225 # that might result in a score less than the opening gap penalty.
 226 # This is the algorithm used by ssearch. Phil Green's algorithm is
 227 # used in the search phase while Miller-Myers algorithm is used in
 228 # the alignment phase
 229 #use constant DPALIGN_LOCAL_GREEN => 2;
 230
 231 BEGIN {
 232     eval {
 233         require Bio::Ext::Align;
 234     };
 235     if ( $@ ) {
 236         die("\nThe C-compiled engine for Smith Waterman alignments (Align) has not been installed.\n Please read the install the bioperl-ext package\n\n");
 237         exit(1);
 238     }
 239 }
 240
 241 sub new {
 242    my ($class, @args) = @_;
 243
 244    my $self = $class->SUPER::new(@args);
 245
 246    my ($match, $mismatch, $gap, $ext, $alg, $matrix) = $self->_rearrange([qw(MATCH
 247                                                                 MISMATCH
 248                                                                 GAP
 249                                                                 EXT
 250                                                                 ALG
 251                                                                 MATRIX
 252                                                         )], @args);
 253
 254    $self->match(3) unless defined $match;
 255    $self->mismatch(-1) unless defined $mismatch;
 256    $self->gap(3) unless defined $gap;
 257    $self->ext(1) unless defined $ext;
 258    $self->alg(DPALIGN_LOCAL_MILLER_MYERS) unless defined $alg;
 259
 260    if (defined $match) {
 261         if ($match =~ /^\d+$/) {
 262             $self->match($match);
 263         }
 264         else {
 265             $self->throw("Match score must be a number, not [$match]");
 266         }
 267     }
 268
 269     if (defined $mismatch) {
 270         if ($match =~ /^\d+$/) {
 271             $self->mismatch($mismatch);
 272         }
 273         else {
 274             $self->throw("Mismatch penalty must be a number, not [$mismatch]");
 275         }
 276     }
 277
 278     if (defined $gap) {
 279         if ($gap =~ /^\d+$/) {
 280             $self->gap($gap);
 281         }
 282         else {
 283             $self->throw("Gap penalty must be a number, not [$gap]");
 284         }
 285     }
 286
 287     if (defined $ext) {
 288         if ($ext =~ /^\d+$/) {
 289             $self->ext($ext);
 290         }
 291         else {
 292             $self->throw("Extension penalty must be a number, not [$ext]");
 293         }
 294     }
 295
 296     if (defined $alg) {
 297         if ($alg == DPALIGN_LOCAL_MILLER_MYERS or $alg == DPALIGN_GLOBAL_MILLER_MYERS or $alg == DPALIGN_ENDSFREE_MILLER_MYERS) {
 298             $self->alg($alg);
 299         }
 300         else {
 301             $self->throw("Algorithm must be either 1, 2 or 3");
 302         }
 303     }
 304
 305     if (defined $matrix and $matrix->isa('Bio::Matrix::MatrixI')) {
 306         $self->{'matrix'} = Bio::Ext::Align::ScoringMatrix->new(join("", $matrix->row_names), $self->gap, $self->ext);
 307         foreach $rowname ($matrix->row_names) {
 308             foreach $colname ($matrix->column_names) {
 309                 Bio::Ext::Align::ScoringMatrix->set_entry($self->{'matrix'}, $rowname, $colname, $matrix->entry($rowname, $colname));
 310             }
 311         }
 312     }
 313     else {
 314         $self->{'matrix'} = 0;
 315     }
 316
 317     return $self;
 318 }
 319
 320 =head2 sequence_profile
 321
 322  Title   : sequence_profile
 323  Usage   : $prof = $factory->sequence_profile($seq1)
 324  Function: Makes a dpAlign_SequenceProfile object from one sequence
 325  Returns : A dpAlign_SequenceProfile object
 326  Args    : The lone argument is a Bio::PrimarySeqI that we want to
 327            build a profile for. Usually, this would be the Query sequence
 328
 329 =cut
 330
 331 sub sequence_profile {
 332     my ($self, $seq1) = @_;
 333
 334     if( ! defined $seq1 || ! $seq1->isa('Bio::PrimarySeqI')) {
 335         $self->warn("Cannot call sequence_profilewithout specifing one sequence (Bio::PrimarySeqI object)");
 336         return;
 337     }
 338
 339     # fix Jitterbug #1044
 340     if( $seq1->length() < 2) {
 341         $self->warn("cannot create sequence profile with length less than 2");
 342         return;
 343     }
 344     if ($seq1->isa('Bio::LocatableSeq')) {
 345        my $seqstr = $seq1->seq;
 346        $seqstr =~ s/\-//g;
 347        $seq1 = Bio::Seq->new(-id => $seq1->id, -seq => $seqstr, -alphabet => $seq1->alphabet);
 348     }
 349     # create engine objects
 350     $seq1->display_id('seq1') unless ( defined $seq1->id() );
 351
 352     if ($seq1->alphabet eq 'dna') {
 353         return Bio::Ext::Align::SequenceProfile->dna_new($seq1->seq, $self->{'match'}, $self->{'mismatch'}, $self->{'gap'}, $self->{'ext'});
 354     }
 355     elsif ($seq1->alphabet eq 'protein') {
 356         return Bio::Ext::Align::SequenceProfile->protein_new($seq1->seq, $self->{'matrix'});
 357     }
 358     else {
 359         croak("There is currently no support for the types of sequences you want to align!\n");
 360         return;
 361     }
 362 }
 363
 364 =head2 pairwise_alignment_score
 365
 366  Title   : pairwise_alignment_score
 367  Usage   : $score = $factory->pairwise_alignment_score($prof,$seq2)
 368  Function: Makes a SimpleAlign object from two sequences
 369  Returns : An integer that is the score of the optimal alignment.
 370  Args    : The first argument is the sequence profile obtained from a
 371            call to the sequence_profile function. The second argument
 372            is a Bio::PrimarySeqI object to be aligned. The second argument
 373            is usually a sequence in the database sequence. Note
 374            that this function only uses Phil Green's algorithm and
 375            therefore theoretically may not always give you the optimal
 376            score.
 377
 378 =cut
 379
 380 sub pairwise_alignment_score {
 381     my ($self, $prof, $seq2) = @_;
 382
 383     if( ! defined $prof || ! $prof->isa('Bio::Ext::Align::SequenceProfile') ||
 384         ! defined $seq2 || ! $seq2->isa('Bio::PrimarySeqI') ) {
 385         $self->warn("Cannot call pairwise_alignment_score without specifing 2 sequences (Bio::PrimarySeqI objects)");
 386         return;
 387     }
 388     # fix Jitterbug #1044
 389     if( $seq2->length() < 2) {
 390         $self->warn("cannot align sequences with length less than 2");
 391         return;
 392     }
 393     if ($seq2->isa('Bio::LocatableSeq')) {
 394        my $seqstr = $seq2->seq;
 395        $seqstr =~ s/\-//g;
 396        $seq2 = Bio::Seq->new(-id => $seq2->id, -seq => $seqstr, -alphabet => $seq2->alphabet);
 397     }
 398     $self->set_memory_and_report();
 399     # create engine objects
 400     $seq2->display_id('seq2') unless ( defined $seq2->id() );
 401
 402     if ($prof->alphabet eq 'dna' and $seq2->alphabet eq 'dna') {
 403         return Bio::Ext::Align::Score_DNA_Sequences($prof, $seq2->seq);
 404     }
 405     elsif ($prof->alphabet eq 'protein' and $seq2->alphabet eq 'protein') {
 406         return Bio::Ext::Align::Score_Protein_Sequences($prof, $seq2->seq);
 407     }
 408     else {
 409         $self->throw("There is currently no support for the types of sequences you want to align!\n");
 410         return;
 411     }
 412 }
 413
 414 =head2 pairwise_alignment
 415
 416  Title   : pairwise_alignment
 417  Usage   : $aln = $factory->pairwise_alignment($seq1,$seq2)
 418  Function: Makes a SimpleAlign object from two sequences
 419  Returns : A SimpleAlign object if there is an alignment with positive
 420            score. Otherwise, return undef.
 421  Args    : The first and second arguments are both Bio::PrimarySeqI
 422            objects that are to be aligned.
 423
 424 =cut
 425
 426 sub pairwise_alignment {
 427     my ($self, $seq1, $seq2) = @_;
 428     my ($aln, $out);
 429
 430     if( ! defined $seq1 || ! $seq1->isa('Bio::PrimarySeqI') ||
 431         ! defined $seq2 || ! $seq2->isa('Bio::PrimarySeqI') ) {
 432         $self->warn("Cannot call pairwise_alignment without specifing 2 sequences (Bio::PrimarySeqI objects)");
 433         return;
 434     }
 435     # fix Jitterbug #1044
 436     if( $seq1->length() < 2 ||
 437         $seq2->length() < 2 ) {
 438         $self->warn("cannot align sequences with length less than 2");
 439         return;
 440     }
 441     if ($seq1->isa('Bio::LocatableSeq')) {
 442        my $seqstr = $seq1->seq;
 443        $seqstr =~ s/\-//g;
 444        $seq1 = Bio::Seq->new(-id => $seq1->id, -seq => $seqstr, -alphabet => $seq1->alphabet);
 445     }
 446     if ($seq2->isa('Bio::LocatableSeq')) {
 447        my $seqstr = $seq2->seq;
 448        $seqstr =~ s/\-//g;
 449        $seq2 = Bio::Seq->new(-id => $seq2->id, -seq => $seqstr, -alphabet => $seq2->alphabet);
 450     }
 451     $self->set_memory_and_report();
 452     # create engine objects
 453     $seq1->display_id('seq1') unless ( defined $seq1->id() );
 454     $seq2->display_id('seq2') unless ( defined $seq2->id() );
 455
 456     if ($seq1->alphabet eq 'dna' and $seq2->alphabet eq 'dna') {
 457         $aln = Bio::Ext::Align::Align_DNA_Sequences($seq1->seq, $seq2->seq, $self->{'match'}, $self->{'mismatch'}, $self->{'gap'}, $self->{'ext'}, $self->{'alg'});
 458     }
 459     elsif ($seq1->alphabet eq 'protein' and $seq2->alphabet eq 'protein') {
 460         $aln = Bio::Ext::Align::Align_Protein_Sequences($seq1->seq, $seq2->seq, $self->{'matrix'}, $self->{'alg'});
 461     }
 462     else {
 463         croak("There is currently no support for the types of sequences you want to align!\n");
 464         return;
 465     }
 466
 467     if (not defined $aln or $aln == 0) {
 468         return;
 469     }
 470
 471     $out = Bio::SimpleAlign->new();
 472
 473     $out->add_seq(Bio::LocatableSeq->new(-seq => $aln->aln1,
 474                                          -start => $aln->start1,
 475                                          -end => $aln->end1,
 476                                          -id => $seq1->id));
 477
 478     $out->add_seq(Bio::LocatableSeq->new(-seq => $aln->aln2,
 479                                          -start => $aln->start2,
 480                                          -end => $aln->end2,
 481                                          -id => $seq2->id));
 482     $out->score($aln->score);
 483     return $out;
 484 }
 485
 486 =head2 align_and_show
 487
 488  Title   : align_and_show
 489  Usage   : $factory->align_and_show($seq1,$seq2,STDOUT)
 490
 491 =cut
 492
 493 sub align_and_show {
 494     my ($self, $seq1, $seq2, $fh) = @_;
 495     my ($aln, $out);
 496
 497     if (! defined $fh) {
 498         $fh = \*STDOUT;
 499     }
 500     if( ! defined $seq1 || ! $seq1->isa('Bio::PrimarySeqI') ||
 501         ! defined $seq2 || ! $seq2->isa('Bio::PrimarySeqI') ) {
 502         $self->warn("Cannot call pairwise_alignment without specifing 2 sequences (Bio::PrimarySeqI objects)");
 503         return;
 504     }
 505     # fix Jitterbug #1044
 506     if( $seq1->length() < 2 ||
 507         $seq2->length() < 2 ) {
 508         $self->warn("cannot align sequences with length less than 2");
 509         return;
 510     }
 511     if ($seq1->isa('Bio::LocatableSeq')) {
 512        my $seqstr = $seq1->seq;
 513        $seqstr =~ s/\-//g;
 514        $seq1 = Bio::Seq->new(-id => $seq1->id, -seq => $seqstr, -alphabet => $seq1->alphabet);
 515     }
 516     if ($seq2->isa('Bio::LocatableSeq')) {
 517        my $seqstr = $seq2->seq;
 518        $seqstr =~ s/\-//g;
 519        $seq2 = Bio::Seq->new(-id => $seq2->id, -seq => $seqstr, -alphabet => $seq2->alphabet);
 520     }
 521     $self->set_memory_and_report();
 522     # create engine objects
 523     $seq1->display_id('seq1') unless ( defined $seq1->id() );
 524     $seq2->display_id('seq2') unless ( defined $seq2->id() );
 525
 526     if ($seq1->alphabet eq 'dna' and $seq2->alphabet eq 'dna') {
 527         $aln = Bio::Ext::Align::Align_DNA_Sequences($seq1->seq, $seq2->seq, $self->{'match'}, $self->{'mismatch'}, $self->{'gap'}, $self->{'ext'}, $self->{'alg'});
 528     }
 529     elsif ($seq1->alphabet eq 'protein' and $seq2->alphabet eq 'protein') {
 530         $aln = Bio::Ext::Align::Align_Protein_Sequences($seq1->seq, $seq2->seq, $self->{'matrix'}, $self->{'alg'});
 531     }
 532     else {
 533         croak("There is currently no support for the types of sequences you want to align!\n");
 534     }
 535
 536     $out = Bio::Ext::Align::AlnBlock->new();
 537     my $s1 = Bio::Ext::Align::AlnSequence->new();
 538     my $s2 = Bio::Ext::Align::AlnSequence->new();
 539     my $a1 = $aln->aln1;
 540     my $a2 = $aln->aln2;
 541     my $first_col = undef;
 542     my $last_col = undef;
 543     my $col;
 544     my $alu1;
 545     my $alu2;
 546     my $g1 = 0;
 547     my $g2 = 0;
 548
 549 # construct AlnBlock
 550     for (my $i = 0; $i < length($a1); ++$i) {
 551         $col = Bio::Ext::Align::AlnColumn->new();
 552         $alu1 = Bio::Ext::Align::AlnUnit->new();
 553         $alu2 = Bio::Ext::Align::AlnUnit->new();
 554         $first_col = $col unless defined $first_col;
 555         Bio::Ext::Align::AlnColumn::set_next($last_col, $col) if defined $last_col;
 556
 557         if (substr($a1, $i, 1) eq "-") {
 558             Bio::Ext::Align::AlnUnit::set_text_label($alu1, "INSERT");
 559             Bio::Ext::Align::AlnUnit::set_text_label($alu2, "SEQUENCE");
 560             ++$g1;
 561         }
 562         elsif (substr($a2, $i, 1) eq "-") {
 563             Bio::Ext::Align::AlnUnit::set_text_label($alu1, "SEQUENCE");
 564             Bio::Ext::Align::AlnUnit::set_text_label($alu2, "INSERT");
 565             ++$g2;
 566         }
 567         else {
 568             Bio::Ext::Align::AlnUnit::set_text_label($alu1, "SEQUENCE");
 569             Bio::Ext::Align::AlnUnit::set_text_label($alu2, "SEQUENCE");
 570         }
 571
 572         Bio::Ext::Align::AlnUnit::set_start($alu1, $aln->start1+$i-$g1-2);
 573         Bio::Ext::Align::AlnUnit::set_end($alu1, $aln->start1+$i-$g1-2);
 574         Bio::Ext::Align::AlnUnit::set_start($alu2, $aln->start2+$i-$g2-2);
 575         Bio::Ext::Align::AlnUnit::set_end($alu2, $aln->start2+$i-$g2-2);
 576         Bio::Ext::Align::AlnColumn::add_alu($col, $alu1);
 577         Bio::Ext::Align::AlnColumn::add_alu($col, $alu2);
 578         $last_col = $col;
 579     }
 580     Bio::Ext::Align::AlnBlock::set_start($out, $first_col);
 581     $col = Bio::Ext::Align::AlnColumn->new();
 582     $alu1 = Bio::Ext::Align::AlnUnit->new();
 583     $alu2 = Bio::Ext::Align::AlnUnit->new();
 584     Bio::Ext::Align::AlnUnit::set_start($alu1, $aln->end1);
 585     Bio::Ext::Align::AlnUnit::set_end($alu1, $aln->end1);
 586     Bio::Ext::Align::AlnUnit::set_text_label($alu1, "END");
 587     Bio::Ext::Align::AlnUnit::set_start($alu2, $aln->end2);
 588     Bio::Ext::Align::AlnUnit::set_end($alu2, $aln->end2);
 589     Bio::Ext::Align::AlnUnit::set_text_label($alu2, "END");
 590     Bio::Ext::Align::AlnColumn::add_alu($col, $alu1);
 591     Bio::Ext::Align::AlnColumn::add_alu($col, $alu2);
 592     Bio::Ext::Align::AlnColumn::set_next($last_col, $col);
 593
 594     &Bio::Ext::Align::write_pretty_str_align($out,$seq1->id,$seq1->seq,$seq2->id,$seq2->seq,12,50,$fh);
 595 }
 596
 597 =head2 match
 598
 599  Title     : match
 600  Usage     : $match = $factory->match() #get
 601            : $factory->match($value) #set
 602  Function  : the set get for the match score
 603  Example   :
 604  Returns   : match value
 605  Arguments : new value
 606
 607 =cut
 608
 609 sub match {
 610     my ($self,$val) = @_;
 611
 612
 613     if( defined $val ) {
 614         if( $val < 0 ) {    # Fixed so that match==0 is allowed /AE
 615             $self->throw("Can't have a match score less than 0");
 616         }
 617         $self->{'match'} = $val;
 618     }
 619     return $self->{'match'};
 620 }
 621
 622 =head2 mismatch
 623
 624  Title     : mismatch
 625  Usage     : $mismatch = $factory->mismatch() #get
 626            : $factory->mismatch($value) #set
 627  Function  : the set get for the mismatch penalty
 628  Example   :
 629  Returns   : mismatch value
 630  Arguments : new value
 631
 632 =cut
 633
 634 sub mismatch {
 635     my ($self,$val) = @_;
 636
 637
 638     if( defined $val ) {
 639         if( $val > 0 ) {    # Fixed so that mismatch==0 is allowed /AE
 640             $self->throw("Can't have a mismatch penalty greater than 0");
 641         }
 642         $self->{'mismatch'} = $val;
 643     }
 644     return $self->{'mismatch'};
 645 }
 646
 647
 648 =head2 gap
 649
 650  Title     : gap
 651  Usage     : $gap = $factory->gap() #get
 652            : $factory->gap($value) #set
 653  Function  : the set get for the gap penalty
 654  Example   :
 655  Returns   : gap value
 656  Arguments : new value
 657
 658 =cut
 659
 660 sub gap {
 661     my ($self,$val) = @_;
 662
 663
 664     if( defined $val ) {
 665         if( $val < 0 ) {    # Fixed so that gap==0 is allowed /AE
 666             $self->throw("Can't have a gap penalty less than 0");
 667         }
 668         $self->{'gap'} = $val;
 669     }
 670     return $self->{'gap'};
 671 }
 672
 673 =head2 ext
 674
 675  Title     : ext
 676  Usage     : $ext = $factory->ext() #get
 677            : $factory->ext($value) #set
 678  Function  : the set get for the ext penalty
 679  Example   :
 680  Returns   : ext value
 681  Arguments : new value
 682
 683 =cut
 684
 685 sub ext {
 686     my ($self,$val) = @_;
 687
 688     if( defined $val ) {
 689         if( $val < 0 ) {    # Fixed so that ext==0 is allowed /AE
 690             $self->throw("Can't have a extension penalty less than 0");
 691         }
 692         $self->{'ext'} = $val;
 693     }
 694     return $self->{'ext'};
 695 }
 696
 697 =head2 alg
 698
 699  Title     : alg
 700  Usage     : $alg = $factory->alg() #get
 701            : $factory->alg($value) #set
 702  Function  : the set get for the algorithm
 703  Example   :
 704  Returns   : alg value
 705  Arguments : new value
 706
 707 =cut
 708
 709 sub alg {
 710     my ($self,$val) = @_;
 711
 712     if( defined $val ) {
 713         if( $val != DPALIGN_LOCAL_MILLER_MYERS and $val != DPALIGN_GLOBAL_MILLER_MYERS and $val != DPALIGN_ENDSFREE_MILLER_MYERS) {
 714             $self->throw("Can't have an algorithm that is not 1, 2 or 3");
 715         }
 716         $self->{'alg'} = $val;
 717     }
 718     return $self->{'alg'};
 719 }
 720
 721 1;