← Back to branch summary

~ubuntu-branches/ubuntu/trusty/bioperl/trusty-proposed

« back to all changes in this revision

Viewing changes to Bio/Search/HSP/GenericHSP.pm

  • Committer: Bazaar Package Importer
  • Author(s): Charles Plessy
  • Date: 2009-03-10 07:19:11 UTC
  • mfrom: (1.2.3 upstream)
  • Revision ID: james.westby@ubuntu.com-20090310071911-fukqzw54pyb1f0bd
Tags: 1.6.0-2
https://launchpad.net/bugs/324001
* Removed patch system (not used):
  - removed instuctions in debian/rules;
  - removed quilt from Build-Depends in debian/control.
* Re-enabled tests:
  - uncommented test command in debian/rules;
  - uncommented previously missing build-dependencies in debian/control.
  - Re-enabled tests and uncommented build-dependencies accordingly.
* Removed libmodule-build-perl and libtest-harness-perl from
  Build-Depends-Indep (provided by perl-modules).
* Better cleaning of empty directories using find -type d -empty -delete
  instead of rmdir in debian/rules (LP: #324001).

Show diffs side-by-side

added added

removed removed

Lines of Context:
Bio/Search/HSP/GenericHSP.pm
1
 
# $Id: GenericHSP.pm,v 1.68.4.5 2006/10/02 23:10:24 sendu Exp $
 
1
# $Id: GenericHSP.pm 15084 2008-12-03 22:31:23Z cjfields $
2
2
#
3
3
# BioPerl module for Bio::Search::HSP::GenericHSP
4
4
#
16
16
 
17
17
=head1 SYNOPSIS
18
18
 
19
 
    my $hsp = new Bio::Search::HSP::GenericHSP( -algorithm => 'blastp',
 
19
    my $hsp = Bio::Search::HSP::GenericHSP->new( -algorithm => 'blastp',
20
20
                                                -evalue    => '1e-30',
21
21
                                                );
22
22
 
47
47
# TODO: Describe how to configure a SearchIO stream so that it generates
48
48
#       GenericHSP objects.
49
49
 
50
 
 
51
50
=head1 DESCRIPTION
52
51
 
53
52
This implementation is "Generic", meaning it is is suitable for
103
102
# Let the code begin...
104
103
 
105
104
package Bio::Search::HSP::GenericHSP;
106
 
use vars qw($GAP_SYMBOL);
107
105
use strict;
108
106
 
109
107
use Bio::Root::Root;
111
109
 
112
110
use base qw(Bio::Search::HSP::HSPI);
113
111
 
114
 
BEGIN {
115
 
    $GAP_SYMBOL = '-';
116
 
}
117
112
=head2 new
118
113
 
119
114
 Title   : new
120
 
 Usage   : my $obj = new Bio::Search::HSP::GenericHSP();
 
115
 Usage   : my $obj = Bio::Search::HSP::GenericHSP->new();
121
116
 Function: Builds a new Bio::Search::HSP::GenericHSP object
122
117
 Returns : Bio::Search::HSP::GenericHSP
123
118
 Args    : -algorithm => algorithm used (BLASTP, TBLASTX, FASTX, etc)
153
148
           -rank        => HSP rank
154
149
           -links       => HSP links information (WU-BLAST only)
155
150
           -hsp_group   => HSP Group informat (WU-BLAST only)
 
151
           -gap_symbol  => symbol representing a gap (default = '-')
 
152
           -hit_features=> string of features found in or near HSP hit
 
153
                           region (reported in some BLAST text output,
 
154
                           v. 2.2.13 and up)
 
155
           -stranded    => If the algorithm isn't known (i.e. defaults to
 
156
                           'generic'), setting this will indicate start/end
 
157
                           coordinates are to be used to determine the strand
 
158
                           for 'query', 'subject', 'hit', 'both', or 'none'
 
159
                           (default = 'none')
 
160
 
156
161
=cut
157
162
 
158
163
sub new {
159
 
    my($class,@args) = @_;
 
164
    my($class,%args) = @_;
160
165
 
161
 
    # don't pass anything to SUPER; complex heirarchy results in lots of work
 
166
    # don't pass anything to SUPER; complex hierarchy results in lots of work
162
167
    # for nothing
 
168
    
163
169
    my $self = $class->SUPER::new();
164
170
 
165
171
    # for speed, don't use _rearrange and just store all input data directly
166
172
    # with no method calls and no work done. work can be carried
167
173
    # out just-in-time later if desired
168
 
    my %args = @args;
169
174
    while (my ($arg, $value) = each %args) {
170
175
        $arg =~ tr/a-z\055/A-Z/d;
171
176
        $self->{$arg} = $value;
173
178
    my $bits = $self->{BITS};
174
179
 
175
180
    defined $self->{VERBOSE} && $self->verbose($self->{VERBOSE});
176
 
 
 
181
    if (exists $self->{GAP_SYMBOL}) {
 
182
        # not checking anything else but the length (must be 1 as only one gap
 
183
        # symbol allowed currently); can add in support for symbol checks or
 
184
        # multiple symbols later if needed
 
185
        $self->throw("Gap symbol must be of length 1") if
 
186
            CORE::length($self->{GAP_SYMBOL}) != 1;
 
187
    } else {
 
188
        # dafault
 
189
        $self->{GAP_SYMBOL} = '-';
 
190
    }
177
191
    $self->{ALGORITHM} ||= 'GENERIC';
178
 
 
 
192
    $self->{STRANDED} ||= 'NONE';
 
193
    
179
194
    if (! defined $self->{QUERY_LENGTH} || ! defined $self->{HIT_LENGTH}) {
180
195
        $self->throw("Must define hit and query length");
181
196
    }
186
201
    return $self;
187
202
}
188
203
 
189
 
 
190
204
sub _logical_length {
191
205
    my ($self, $type) = @_;
192
206
    my $algo = $self->algorithm;
257
271
 
258
272
=cut
259
273
 
260
 
sub evalue { shift->significance(@_) }
261
 
 
262
 
# Override significance to return the e-value or, if this is
263
 
# not defined (WU-BLAST), return the p-value.
264
 
sub significance {
265
 
    my $self = shift;
266
 
    my $signif = $self->query->significance(@_);
267
 
    return (defined $signif && $signif ne '') ? $signif : $self->pvalue(@_);
 
274
sub evalue {
 
275
    my ($self,$value) = @_;
 
276
    my $previous = $self->{'EVALUE'};
 
277
    if( defined $value  ) {
 
278
        $self->{'EVALUE'} = $value;
 
279
    }
 
280
    return $previous;
268
281
}
269
282
 
270
 
 
271
283
=head2 frac_identical
272
284
 
273
285
 Title   : frac_identical
277
289
 Args    : arg 1:  'query' = num identical / length of query seq (without gaps)
278
290
                   'hit'   = num identical / length of hit seq (without gaps)
279
291
                             synonyms: 'sbjct', 'subject'
280
 
                   'total' = num conserved / length of alignment (with gaps)
 
292
                   'total' = num identical / length of alignment (with gaps)
281
293
                             synonyms: 'hsp'
282
294
                   default = 'total'
283
295
           arg 2: [optional] frac identical value to set for the type requested
349
361
 
350
362
 Title    : gaps
351
363
 Usage    : my $gaps = $hsp->gaps( ['query'|'hit'|'total'] );
352
 
 Function : Get the number of gaps in the query, hit, or total alignment.
 
364
 Function : Get the number of gap characters in the query, hit, or total alignment.
353
365
 Returns  : Integer, number of gaps or 0 if none
354
 
 Args     : arg 1: 'query' = num gaps in query seq
355
 
                   'hit'   = num gaps in hit seq; synonyms: 'sbjct', 'subject'
356
 
                   'total' = num gaps in whole alignment;  synonyms: 'hsp'
 
366
 Args     : arg 1: 'query' = num gap characters in query seq
 
367
                   'hit'   = num gap characters in hit seq; synonyms: 'sbjct', 'subject'
 
368
                   'total' = num gap characters in whole alignment;  synonyms: 'hsp'
357
369
                   default = 'total'
358
370
            arg 2: [optional] integer gap value to set for the type requested
359
371
 
360
372
=cut
361
373
 
362
374
sub gaps {
363
 
    my ($self, $type,$value) = @_;
 
375
    my ($self, $type, $value) = @_;
364
376
 
365
377
    unless ($self->{_did_pregaps}) {
366
378
        $self->_pre_gaps;
527
539
=head2 frame
528
540
 
529
541
 Title   : frame
530
 
 Usage   : $hsp->frame($queryframe,$subjectframe)
 
542
 Usage   : my ($qframe, $hframe) = $hsp->frame('list',$queryframe,$subjectframe)
531
543
 Function: Set the Frame for both query and subject and insure that
532
544
           they agree.
533
545
           This overrides the frame() method implementation in
534
546
           FeaturePair.
535
 
 Returns : array of query and subjects if return type wants an array
536
 
           or query frame if defined or subject frame
537
 
 Args    : none
 
547
 Returns : array of query and subject frame if return type wants an array
 
548
           or query frame if defined or subject frame if not defined
 
549
 Args    : 'hit' or 'subject' or 'sbjct' to retrieve the frame of the subject (default)
 
550
           'query' to retrieve the query frame 
 
551
           'list' or 'array' to retrieve both query and hit frames together
538
552
 Note    : Frames are stored in the GFF way (0-2) not 1-3
539
553
           as they are in BLAST (negative frames are deduced by checking
540
554
                                 the strand of the query or hit)
541
555
 
542
556
=cut
543
557
 
 
558
# Note: changed 4/19/08 - bug 2485
 
559
#
 
560
# frame() is supposed to be a getter/setter (as is implied by the Function desc
 
561
# above; this is also consistent with Bio::SeqFeature::SimilarityPair).  Also,
 
562
# the API is not consistent with the other HSP/SimilarityPair methods such as
 
563
# strand(), start(), end(), etc. 
 
564
#
 
565
# In order to make this consistent with other methods all work is now done
 
566
# when the features are instantiated and not delayed.  We compromise by
 
567
# defaulting back 'to hit' w/o passed args.  Setting is now allowed.  
 
568
 
544
569
sub frame {
545
570
    my $self = shift;
546
 
 
547
 
    unless (defined $self->{_did_preframe}) {
548
 
        $self->_pre_frame;
549
 
    }
550
 
    my $qframe = $self->{QUERY_FRAME};
551
 
    my $sframe = $self->{HIT_FRAME};
552
 
 
553
 
    if( defined $qframe ) {
554
 
        if( $qframe == 0 ) {
555
 
            $qframe = 0;
556
 
        } elsif( $qframe !~ /^([+-])?([1-3])/ ) {
557
 
            $self->warn("Specifying an invalid query frame ($qframe)");
558
 
            $qframe = undef;
559
 
        } else {
560
 
            my $dir = $1;
561
 
            $dir = '+' unless defined $dir;
562
 
            if( ($dir eq '-' && $self->query->strand >= 0) ||
563
 
                ($dir eq '+' && $self->query->strand <= 0) ) {
564
 
                $self->warn("Query frame ($qframe) did not match strand of query (". $self->query->strand() . ")");
565
 
            }
566
 
            # Set frame to GFF [0-2] -
567
 
            # what if someone tries to put in a GFF frame!
568
 
            $qframe = $2 - 1;
569
 
        }
570
 
        $self->query->frame($qframe);
571
 
    }
572
 
    if( defined $sframe ) {
573
 
          if( $sframe == 0 ) {
574
 
            $sframe = 0;
575
 
          } elsif( $sframe !~ /^([+-])?([1-3])/ ) {
576
 
            $self->warn("Specifying an invalid subject frame ($sframe)");
577
 
            $sframe = undef;
578
 
          } else {
579
 
              my $dir = $1;
580
 
              $dir = '+' unless defined $dir;
581
 
              if( ($dir eq '-' && $self->hit->strand >= 0) ||
582
 
                  ($dir eq '+' && $self->hit->strand <= 0) )
583
 
              {
584
 
                  $self->warn("Subject frame ($sframe) did not match strand of subject (". $self->hit->strand() . ")");
585
 
              }
586
 
 
587
 
              # Set frame to GFF [0-2]
588
 
              $sframe = $2 - 1;
589
 
          }
590
 
          $self->hit->frame($sframe);
591
 
      }
592
 
    if (wantarray() && $self->algorithm =~ /^T(BLAST|FAST)(X|Y|XY)/oi)
593
 
    {
594
 
        return ($self->query->frame(), $self->hit->frame());
595
 
    } elsif (wantarray())  {
596
 
        ($self->query->frame() &&
597
 
         return ($self->query->frame(), undef)) ||
598
 
             ($self->hit->frame() &&
599
 
              return (undef, $self->hit->frame()));
600
 
    } else {
601
 
        ($self->query->frame() &&
602
 
         return $self->query->frame()) ||
603
 
        ($self->hit->frame() &&
604
 
         return $self->hit->frame());
605
 
    }
 
571
    my $val = shift;
 
572
    if (!defined $val) {
 
573
        # unfortunately, w/o args we need to warn about API changes
 
574
        $self->warn("API for frame() has changed.\n".
 
575
                    "Please refer to documentation for Bio::Search::HSP::GenericHSP;\n".
 
576
                    "returning query frame");
 
577
        $val = 'query';
 
578
    }
 
579
    $val =~ s/^\s+//;
 
580
 
 
581
    if( $val =~ /^q/i ) {
 
582
        return $self->query->frame(@_);
 
583
    } elsif( $val =~ /^hi|^s/i ) {
 
584
        return $self->hit->frame(@_);
 
585
    } elsif (  $val =~ /^list|array/i ) {
 
586
        return ($self->query->frame($_[0]), 
 
587
            $self->hit->frame($_[1]) );
 
588
    } elsif ( $val =~ /^\d+$/) {
 
589
        # old API i.e. frame($query_frame, $hit_frame). This catches all but one
 
590
        # case, where no arg is passed (so the hit is wanted).  
 
591
        $self->warn("API for frame() has changed.\n".
 
592
                    "Please refer to documentation for Bio::Search::HSP::GenericHSP");
 
593
        wantarray ? 
 
594
        return ($self->query->frame($val), 
 
595
            $self->hit->frame(@_) ) :
 
596
        return $self->hit->frame($val,@_);
 
597
    } else { 
 
598
        $self->warn("unrecognized component '$val' requested\n");
 
599
    }
 
600
    return 0;
606
601
}
607
602
 
608
 
 
609
603
=head2 get_aln
610
604
 
611
605
 Title   : get_aln
620
614
    my ($self) = @_;
621
615
    require Bio::LocatableSeq;
622
616
    require Bio::SimpleAlign;
623
 
    my $aln = new Bio::SimpleAlign;
 
617
    my $aln = Bio::SimpleAlign->new();
624
618
    my $hs = $self->hit_string();
625
619
    my $qs = $self->query_string();
626
620
    # FASTA specific stuff moved to the FastaHSP object
628
622
    $seqonly =~ s/[\-\s]//g;
629
623
    my ($q_nm,$s_nm) = ($self->query->seq_id(),
630
624
                        $self->hit->seq_id());
631
 
    unless( defined $q_nm && CORE::length ($q_nm) ) {
632
 
        $q_nm = 'query';
633
 
    }
634
 
    unless( defined $s_nm && CORE::length ($s_nm) ) {
635
 
        $s_nm = 'hit';
636
 
    }
637
 
    my $query = new Bio::LocatableSeq('-seq'   => $qs,
 
625
    # Should we silently change the name of the query or hit if it isn't
 
626
    # defined?  May need revisiting... cjfields 2008-12-3 (commented out below)
 
627
    
 
628
    #unless( defined $q_nm && CORE::length ($q_nm) ) {
 
629
    #    $q_nm = 'query';
 
630
    #}
 
631
    #unless( defined $s_nm && CORE::length ($s_nm) ) {
 
632
    #    $s_nm = 'hit';
 
633
    #}
 
634
    
 
635
    # mapping: 1 residues for every x coordinate positions
 
636
    my $query = Bio::LocatableSeq->new('-seq'   => $qs,
638
637
                                      '-id'    => $q_nm,
639
638
                                      '-start' => $self->query->start,
640
639
                                      '-end'   => $self->query->end,
 
640
                                      '-force_nse' => $q_nm ? 0 : 1,
 
641
                                      '-mapping' => [1, $self->{_query_mapping}]
641
642
                                      );
642
643
    $seqonly = $hs;
643
644
    $seqonly =~ s/[\-\s]//g;
644
 
    my $hit =  new Bio::LocatableSeq('-seq'    => $hs,
 
645
    my $hit =  Bio::LocatableSeq->new('-seq'    => $hs,
645
646
                                      '-id'    => $s_nm,
646
647
                                      '-start' => $self->hit->start,
647
648
                                      '-end'   => $self->hit->end,
 
649
                                      '-force_nse' => $s_nm ? 0 : 1,
 
650
                                      '-mapping' => [1, $self->{_hit_mapping}]
648
651
                                      );
649
652
    $aln->add_seq($query);
650
653
    $aln->add_seq($hit);
656
659
 Title   : num_conserved
657
660
 Usage   : $obj->num_conserved($newval)
658
661
 Function: returns the number of conserved residues in the alignment
659
 
 Returns : inetger
 
662
 Returns : integer
660
663
 Args    : integer (optional)
661
664
 
662
665
 
717
720
    return $self->{RANK};
718
721
}
719
722
 
720
 
 
721
723
=head2 seq_inds
722
724
 
723
725
 Title   : seq_inds
731
733
           : May include ranges if collapse is true.
732
734
 Argument  : seq_type  = 'query' or 'hit' or 'sbjct'  (default = query)
733
735
           :  ('sbjct' is synonymous with 'hit')
734
 
           : class     = 'identical' or 'conserved' or 'nomatch' or 'gap'
735
 
           :              (default = identical)
736
 
           :              (can be shortened to 'id' or 'cons')
737
 
           :             or 'conserved-not-identical'
 
736
           : class     = 'identical' - identical positions
 
737
           :             'conserved' - conserved positions
 
738
           :             'nomatch'   - mismatched residue or gap positions
 
739
           :             'mismatch'  - mismatched residue positions (no gaps)
 
740
           :             'gap'       - gap positions only
 
741
           :             'frameshift'- frameshift positions only
 
742
           :             'conserved-not-identical' - conserved positions w/o 
 
743
           :                            identical residues
 
744
           :             The name can be shortened to 'id' or 'cons' unless
 
745
           :             the name is ambiguous.  The default value is
 
746
           :             'identical'
 
747
           :
738
748
           : collapse  = boolean, if true, consecutive positions are merged
739
749
           :             using a range notation, e.g., "1 2 3 4 5 7 9 10 11"
740
750
           :             collapses to "1-5 7 9-11". This is useful for
741
751
           :             consolidating long lists. Default = no collapse.
 
752
           :
742
753
 Throws    : n/a.
743
 
 Comments  :
744
 
 
 
754
 Comments  : For HSPs partially or completely derived from translated sequences
 
755
           : (TBLASTN, BLASTX, TBLASTX, or similar), some positional data
 
756
           : cannot easily be attributed to a single position (i.e. the 
 
757
           : positional data is ambiguous).  In these cases all three codon 
 
758
           : positions are reported.  Under these conditions you can check 
 
759
           : ambiguous_seq_inds() to determine whether the query, subject, 
 
760
           : or both are ambiguous.
 
761
           :
745
762
See Also   : L<Bio::Search::SearchUtils::collapse_nums()|Bio::Search::SearchUtils>,
746
763
             L<Bio::Search::Hit::HitI::seq_inds()|Bio::Search::Hit::HitI>
747
764
 
767
784
       $self->warn("unknown seqtype $seqType using 'query'");
768
785
       $seqType = 'query';
769
786
   }
 
787
   
770
788
   $t = lc(substr($class,0,1));
771
789
 
772
790
   if( $t eq 'c' ) {
779
797
       $class = 'identical';
780
798
   } elsif( $t eq 'n' ) {
781
799
       $class = 'nomatch';
 
800
   } elsif( $t eq 'm' ) {
 
801
       $class = 'mismatch';
782
802
   } elsif( $t eq 'g' ) {
783
803
       $class = 'gap';
 
804
   } elsif( $t eq 'f' ) {
 
805
       $class = 'frameshift';    
784
806
   } else {
785
807
       $self->warn("unknown sequence class $class using 'identical'");
786
808
       $class = 'identical';
792
814
   my @ary;
793
815
 
794
816
   if( $class eq '_gap' ) {
795
 
       # this means that we are remapping the gap length that is stored
796
 
       # in the hash (for example $self->{'_gapRes_query'} )
797
 
       # so we'll return an array which has the values of the position of the
798
 
       # of the gap (the key in the hash) + the gap length (value in the
799
 
       # hash for this key - 1.
800
 
 
801
 
       @ary = map { $_ > 1 ?
802
 
                        $_..($_ + $self->{"${class}Res$seqType"}->{$_} - 1) :
803
 
                        $_ }
804
 
              sort { $a <=> $b } keys %{ $self->{"${class}Res$seqType"}};
 
817
        # this means that we are remapping the gap length that is stored
 
818
        # in the hash (for example $self->{'_gapRes_query'} )
 
819
        # so we'll return an array which has the values of the position of the
 
820
        # of the gap (the key in the hash) + the gap length (value in the
 
821
        # hash for this key - 1.
 
822
        
 
823
        # changing this; since the index is the position prior to the insertion,
 
824
        # repeat the position based on the number of gaps inserted
 
825
        @ary = map { my @tmp;
 
826
                     # position holds number of gaps inserted
 
827
                     for my $g (1..$self->{seqinds}{"${class}Res$seqType"}->{$_}) {
 
828
                        push @tmp, $_ ;
 
829
                     }
 
830
                     @tmp}
 
831
              sort { $a <=> $b } keys %{ $self->{seqinds}{"${class}Res$seqType"}};
805
832
   } elsif( $class eq '_conservedall' ) {
806
833
       @ary = sort { $a <=> $b }
807
 
       keys %{ $self->{"_conservedRes$seqType"}},
808
 
       keys %{ $self->{"_identicalRes$seqType"}},
 
834
       keys %{ $self->{seqinds}{"_conservedRes$seqType"}},
 
835
       keys %{ $self->{seqinds}{"_identicalRes$seqType"}},
809
836
   }  else {
810
 
       @ary = sort { $a <=> $b } keys %{ $self->{"${class}Res$seqType"}};
 
837
       @ary = sort { $a <=> $b } keys %{ $self->{seqinds}{"${class}Res$seqType"}};
811
838
   }
812
839
   require Bio::Search::BlastUtils if $collapse;
813
840
 
814
841
   return $collapse ? &Bio::Search::SearchUtils::collapse_nums(@ary) : @ary;
815
842
}
816
843
 
 
844
=head2 ambiguous_seq_inds
 
845
 
 
846
 Title     : ambiguous_seq_inds
 
847
 Purpose   : returns a string denoting whether sequence indices for query, 
 
848
           : subject, or both are ambiguous
 
849
 Returns   : String; 'query', 'subject', 'query/subject', or empty string ''
 
850
 Argument  : none
 
851
 Comments  : For HSPs partially or completely derived from translated sequences
 
852
           : (TBLASTN, BLASTX, TBLASTX, or similar), some positional data
 
853
           : cannot easily be attributed to a single position (i.e. the 
 
854
           : positional data is ambiguous).  In these cases all three codon 
 
855
           : positions are reported when using seq_inds().  Under these
 
856
           : conditions you can check ambiguous_seq_inds() to determine whether
 
857
           : the query, subject, or both are ambiguous.
 
858
See Also   : L<Bio::Search::Hit::HSPI::seq_inds()>
 
859
 
 
860
=cut
 
861
 
 
862
sub ambiguous_seq_inds {
 
863
    my $self = shift;
 
864
    $self->_calculate_seq_positions();    
 
865
    return $self->{seqinds}{'_warnRes'} if exists $self->{seqinds}{'_warnRes'};
 
866
    return $self->{seqinds}{'_warnRes'} = '';
 
867
}
817
868
 
818
869
=head2 Inherited from L<Bio::SeqFeature::SimilarityPair>
819
870
 
890
941
 Returns : numeric
891
942
 Args    : [optional] new value to set
892
943
 
 
944
=cut 
 
945
 
 
946
# Override significance to return the e-value or, if this is
 
947
# not defined (WU-BLAST), return the p-value.
 
948
sub significance {
 
949
    my ($self, $val) = @_;
 
950
    if (!defined $self->{SIGNIFICANCE} || $val) {
 
951
        $self->query->significance($val) if $val;
 
952
        $self->{SIGNIFICANCE} = $val || $self->evalue || $self->pvalue;
 
953
    }
 
954
    return $self->{SIGNIFICANCE};
 
955
}
893
956
 
894
957
=head2 score
895
958
 
941
1004
    my ($self,@args) = @_;
942
1005
    return unless ( $self->{'_sequenceschanged'} );
943
1006
    $self->{'_sequenceschanged'} = 0;
944
 
    my ($mchar, $schar, $qchar);
945
1007
    my ($seqString, $qseq,$sseq) = ( $self->homology_string(),
946
1008
                                     $self->query_string(),
947
1009
                                     $self->hit_string() );
948
 
 
949
 
    # Using hashes to avoid saving duplicate residue numbers.
950
 
    my %identicalList_query = ();
951
 
    my %identicalList_sbjct = ();
952
 
    my %conservedList_query = ();
953
 
    my %conservedList_sbjct = ();
954
 
 
955
 
    my %gapList_query = ();
956
 
    my %gapList_sbjct = ();
957
 
    my %nomatchList_query = ();
958
 
    my %nomatchList_sbjct = ();
959
 
 
 
1010
    my ($mlen, $qlen, $slen) = (CORE::length($seqString), CORE::length($qseq), CORE::length($sseq));
960
1011
    my $qdir = $self->query->strand || 1;
961
1012
    my $sdir = $self->hit->strand || 1;
962
 
    my $resCount_query = ($qdir >=0) ? $self->query->end : $self->query->start;
963
 
    my $resCount_sbjct = ($sdir >=0) ? $self->hit->end : $self->hit->start;
964
 
 
 
1013
    my ($resCount_query, $endpoint_query) = ($qdir <=0) ? ($self->query->end, $self->query->start)
 
1014
        : ($self->query->start, $self->query->end);
 
1015
    my ($resCount_sbjct, $endpoint_sbjct) = ($sdir <=0) ? ($self->hit->end, $self->hit->start)
 
1016
        : ($self->hit->start, $self->hit->end);
 
1017
    
965
1018
    my $prog = $self->algorithm;
 
1019
    
966
1020
    if( $prog  =~ /FAST|SSEARCH|SMITH-WATERMAN/i ) {
 
1021
    
 
1022
        # we infer the end of the regional sequence where the first and last
 
1023
        # non spaces are in the homology string
 
1024
        # then we use the HSP->length to tell us how far to read
 
1025
        # to cut off the end of the sequence
 
1026
        
 
1027
        my ($start, $rest) = (0,0);
 
1028
        if( $seqString =~ /^(\s+)?(.*?)\s*$/ ) {
 
1029
            ($start, $rest) = ($1 ? CORE::length($1) : 0, CORE::length($2));
 
1030
        }
 
1031
    
 
1032
        $seqString = substr($seqString, $start, $rest);
 
1033
        $qseq = substr($qseq, $start, $rest);
 
1034
        $sseq = substr($sseq, $start, $rest);
 
1035
 
 
1036
        # commented out 10/29/08
 
1037
        # removing frameshift symbols doesn't take into account the following:
 
1038
        # 1) does not remove the same point in the homology string (get
 
1039
        # positional errors)
 
1040
        # 2) adjustments to the overall position in the string due to the
 
1041
        # frameshift must be taken into consideration (get balancing errors)
 
1042
        #
 
1043
        # Frameshifts will be handled directly in the main loop. 
 
1044
        # --chris
 
1045
        
967
1046
        # fasta reports some extra 'regional' sequence information
968
1047
        # we need to clear out first
969
1048
        # this seemed a bit insane to me at first, but it appears to
970
1049
        # work --jason
971
 
 
972
 
        # we infer the end of the regional sequence where the first
973
 
        # non space is in the homology string
974
 
        # then we use the HSP->length to tell us how far to read
975
 
        # to cut off the end of the sequence
976
 
 
977
 
        # one possible problem is the sequence which
978
 
 
979
 
        my ($start) = (0);
980
 
        if( $seqString =~ /^(\s+)/ ) {
981
 
            $start = CORE::length($1);
 
1050
        
 
1051
        #$qseq =~ s![\\\/]!!g;
 
1052
        #$sseq =~ s![\\\/]!!g;
 
1053
    }
 
1054
 
 
1055
    $self->{seqinds}{'_warnRes'} = '';
 
1056
    ($self->{_sbjct_offset}, $self->{_query_offset}) = (1,1);
 
1057
    if($prog =~ /^(?:PSI)?T(BLAST|FAST)(N|X|Y)/oi ) {
 
1058
        $self->{_sbjct_offset} = 3;
 
1059
        if ($1 eq 'BLAST' && $2 eq 'X') { #TBLASTX
 
1060
            $self->{_query_offset} = 3;
 
1061
            $self->{seqinds}{'_warnRes'} = 'query/subject';
 
1062
        } else {
 
1063
            $self->{seqinds}{'_warnRes'} = 'subject';
982
1064
        }
983
 
 
984
 
        $seqString = substr($seqString, $start,$self->length('total'));
985
 
        $qseq = substr($qseq, $start,$self->length('total'));
986
 
        $sseq = substr($sseq, $start,$self->length('total'));
987
 
 
988
 
        $qseq =~ s![\\\/]!!g;
989
 
        $sseq =~ s![\\\/]!!g;
990
 
    }
991
 
 
992
 
    if($prog =~ /^(PSI)?T(BLAST|FAST)N/oi ) {
993
 
        $resCount_sbjct = int($resCount_sbjct / 3);
994
1065
    } elsif($prog =~ /^(BLAST|FAST)(X|Y|XY)/oi  ) {
995
 
        $resCount_query = int($resCount_query / 3);
996
 
    } elsif($prog =~ /^T(BLAST|FAST)(X|Y|XY)/oi ) {
997
 
        $resCount_query = int($resCount_query / 3);
998
 
        $resCount_sbjct = int($resCount_sbjct / 3);
999
 
    }
1000
 
    while( $mchar = chop($seqString) ) {
1001
 
        ($qchar, $schar) = (chop($qseq), chop($sseq));
1002
 
        if( $mchar eq '+' || $mchar eq '.' || $mchar eq ':' ) {
1003
 
            $conservedList_query{ $resCount_query } = 1;
1004
 
            $conservedList_sbjct{ $resCount_sbjct } = 1;
1005
 
        } elsif( $mchar ne ' ' ) {
1006
 
            $identicalList_query{ $resCount_query } = 1;
1007
 
            $identicalList_sbjct{ $resCount_sbjct } = 1;
1008
 
        } elsif( $mchar eq ' ') {
1009
 
            $nomatchList_query{ $resCount_query } = 1;
1010
 
            $nomatchList_sbjct{ $resCount_sbjct } = 1;
1011
 
        }
1012
 
        if( $qchar eq $GAP_SYMBOL ) {
1013
 
            $gapList_query{ $resCount_query } ++;
1014
 
        } else {
1015
 
            $resCount_query -= $qdir;
1016
 
        }
1017
 
        if( $schar eq $GAP_SYMBOL ) {
1018
 
            $gapList_sbjct{ $resCount_query } ++;
1019
 
        } else {
1020
 
            $resCount_sbjct -=$sdir;
1021
 
        }
1022
 
    }
1023
 
    $self->{'_identicalRes_query'} = \%identicalList_query;
1024
 
    $self->{'_conservedRes_query'} = \%conservedList_query;
1025
 
    $self->{'_nomatchRes_query'}   = \%nomatchList_query;
1026
 
    $self->{'_gapRes_query'}       = \%gapList_query;
1027
 
 
1028
 
    $self->{'_identicalRes_sbjct'} = \%identicalList_sbjct;
1029
 
    $self->{'_conservedRes_sbjct'} = \%conservedList_sbjct;
1030
 
    $self->{'_nomatchRes_sbjct'}   = \%nomatchList_sbjct;
1031
 
    $self->{'_gapRes_sbjct'}       = \%gapList_sbjct;
 
1066
        $self->{_query_offset} = 3;
 
1067
        $self->{seqinds}{'_warnRes'} = 'query';
 
1068
    }
 
1069
    my ($qfs, $sfs) = (0,0);
 
1070
    CHAR_LOOP:
 
1071
    for my $pos (0..CORE::length($seqString)-1) {
 
1072
        my @qrange = (0 - $qfs)..($self->{_query_offset} - 1);
 
1073
        my @srange = (0 - $sfs)..($self->{_sbjct_offset} - 1);
 
1074
        # $self->debug("QRange:@qrange SRange:@srange\n") if ($qfs || $sfs);
 
1075
        ($qfs, $sfs) = (0,0);
 
1076
        my ($mchar, $qchar, $schar) = (
 
1077
            unpack("x$pos A1",$seqString) || ' ',
 
1078
            $pos < CORE::length($qseq) ? unpack("x$pos A1",$qseq) : '-',
 
1079
            $pos < CORE::length($sseq) ? unpack("x$pos A1",$sseq) : '-'
 
1080
            );
 
1081
        my $matchtype = '';
 
1082
        my ($qgap, $sgap) = (0,0);
 
1083
        if( $mchar eq '+' || $mchar eq '.') {    # conserved
 
1084
            $self->{seqinds}{_conservedRes_query}{ $resCount_query + ($_ * $qdir) } = 1 for @qrange;
 
1085
            $self->{seqinds}{_conservedRes_sbjct}{ $resCount_sbjct + ($_ * $sdir) } = 1 for @srange;
 
1086
            $matchtype = 'conserved';
 
1087
        } elsif( $mchar eq ':' || $mchar ne ' ' ) { # identical
 
1088
            $self->{seqinds}{_identicalRes_query}{ $resCount_query + ($_ * $qdir) } = 1 for @qrange;
 
1089
            $self->{seqinds}{_identicalRes_sbjct}{ $resCount_sbjct + ($_ * $sdir) } = 1 for @srange;
 
1090
            $matchtype = 'identical';
 
1091
        } elsif( $mchar eq ' ' ) {  # possible mismatch/nomatch/frameshift
 
1092
            $qfs = $qchar eq '/'  ?  -1 : # base inserted to match frame
 
1093
                   $qchar eq '\\' ?   1 : # base deleted to match frame
 
1094
                   0;
 
1095
            $sfs = $schar eq '/'  ?  -1 :
 
1096
                   $schar eq '\\' ?   1 :
 
1097
                   0;
 
1098
            if ($qfs) {
 
1099
                # Frameshifts are tricky.
 
1100
                
 
1101
                # '/' indicates that the next residue is shifted back one
 
1102
                # (-1) frame position and is a deletion of one base (so to
 
1103
                # correctly align, a base is inserted). That residue should only
 
1104
                # occupy two sequence positions instead of three.
 
1105
                
 
1106
                # '\' indicates that the next residue is shifted forward
 
1107
                # one (+1) frame position and is an insertion of one base (so to
 
1108
                # correctly align, a base is removed). That residue should
 
1109
                # occupy four sequence positions instead of three.
 
1110
                
 
1111
                # Note that gaps are not counted across from frameshifts.
 
1112
                # Frameshift indices are a range of positions starting in the
 
1113
                # previous sequence position in which the frameshift occurs;
 
1114
                # they are ambiguous by nature.
 
1115
                $self->{seqinds}{_frameshiftRes_query}{ $resCount_query - ($_ * $qdir * $qfs) } = $qfs for @qrange;
 
1116
                $matchtype = "$qfs frameshift-query";
 
1117
                $sgap = $qgap = 1;
 
1118
            }
 
1119
            elsif ($sfs) {
 
1120
                $self->{seqinds}{_frameshiftRes_sbjct}{ $resCount_sbjct - ($_ * $sdir * $sfs) } = $sfs for @srange;
 
1121
                $matchtype = "$sfs frameshift-sbcjt";
 
1122
                $sgap = $qgap = 1;
 
1123
            }
 
1124
            elsif ($qchar eq $self->{GAP_SYMBOL}) {
 
1125
                # gap are counted as being in the immediately preceeding residue
 
1126
                # position; for translated sequences this is not the start of
 
1127
                # the previous codon, but the third codon position
 
1128
                $self->{seqinds}{_gapRes_query}{ $resCount_query - $qdir }++ for @qrange;
 
1129
                $matchtype = 'gap-query';
 
1130
                $qgap++;
 
1131
            }
 
1132
            elsif ($schar eq $self->{GAP_SYMBOL}) {
 
1133
                $self->{seqinds}{_gapRes_sbjct}{ $resCount_sbjct - $sdir }++ for @srange;
 
1134
                $matchtype = 'gap-sbjct';
 
1135
                $sgap++;
 
1136
            }
 
1137
            else {
 
1138
                # if not a gap or frameshift in either seq, must be mismatch
 
1139
                $self->{seqinds}{_mismatchRes_query}{ $resCount_query + ($_ * $qdir) } = 1 for @qrange;
 
1140
                $self->{seqinds}{_mismatchRes_sbjct}{ $resCount_sbjct + ($_ * $sdir) } = 1 for @srange;
 
1141
                $matchtype = 'mismatch';
 
1142
            }
 
1143
            # always add a nomatch unless the current position in the seq is a gap
 
1144
            if (!$sgap) {
 
1145
                $self->{seqinds}{_nomatchRes_sbjct}{ $resCount_sbjct + ($_ * $sdir) } = 1 for @srange;
 
1146
            }
 
1147
            if (!$qgap) {
 
1148
                $self->{seqinds}{_nomatchRes_query}{ $resCount_query + ($_ * $qdir) } = 1 for @qrange;
 
1149
            }
 
1150
        } else {
 
1151
            $self->warn("Unknown midline character: [$mchar]");
 
1152
        }
 
1153
        # leave in and uncomment for future debugging
 
1154
        #$self->debug(sprintf("%7d %1s[%1s]%1s %-7d Type: %-20s QOff:%-2d SOff:%-2d\n",
 
1155
        #                     $resCount_query,
 
1156
        #                     $qchar,
 
1157
        #                     $mchar,
 
1158
        #                     $schar,
 
1159
        #                     $resCount_sbjct,
 
1160
        #                     $matchtype,
 
1161
        #                     ($self->{_query_offset} * $qdir),
 
1162
        #                     ($self->{_sbjct_offset} * $sdir)));
 
1163
        $resCount_query += ($qdir * (scalar(@qrange) + $qfs)) if !$qgap;
 
1164
        $resCount_sbjct += ($sdir * (scalar(@srange) + $sfs)) if !$sgap;
 
1165
    }
1032
1166
    return 1;
1033
1167
}
1034
1168
 
1040
1174
 
1041
1175
sub n {
1042
1176
    my $self = shift;
1043
 
    if(@_) { $self->{'_n'} = shift; }
1044
 
    defined $self->{'_n'} ? $self->{'_n'} : '';
 
1177
    if(@_) { $self->{'N'} = shift; }
 
1178
    defined $self->{'N'} ? $self->{'N'} : '';
1045
1179
}
1046
1180
 
1047
1181
=head2 range
1107
1241
    return $self->{HSP_GROUP};
1108
1242
}
1109
1243
 
 
1244
=head2 hit_features
 
1245
 
 
1246
 Title   : hit_features
 
1247
 Usage   : $obj->hit_features($newval)
 
1248
 Function: Get/Set the HSP hit feature string (from some BLAST 2.2.13 text
 
1249
           output), which is a string of overlapping or nearby features in HSP
 
1250
           hit
 
1251
 Returns : Value of hit features, if present
 
1252
 Args    : On set, new value (a scalar or undef, optional)
 
1253
 
 
1254
 
 
1255
=cut
 
1256
 
 
1257
sub hit_features {
 
1258
    my $self = shift;
 
1259
 
 
1260
    return $self->{HIT_FEATURES} = shift if @_;
 
1261
    return $self->{HIT_FEATURES};
 
1262
}
 
1263
 
1110
1264
# The cigar string code is written by Juguang Xiao <juguang@fugu-sg.org>
1111
1265
 
1112
1266
=head1 Brief introduction on cigar string
1195
1349
    for(my $i=0; $i <= $#qchars; $i++){
1196
1350
        my $qchar = $qchars[$i];
1197
1351
        my $hchar = $hchars[$i];
1198
 
        if($qchar ne $GAP_SYMBOL && $hchar ne $GAP_SYMBOL){ # Match
 
1352
        if($qchar ne $self->{GAP_SYMBOL} && $hchar ne $self->{GAP_SYMBOL}){ # Match
1199
1353
            $cigar_string .= $self->_sub_cigar_string('M');
1200
 
        }elsif($qchar eq $GAP_SYMBOL){ # Deletion
 
1354
        }elsif($qchar eq $self->{GAP_SYMBOL}){ # Deletion
1201
1355
            $cigar_string .= $self->_sub_cigar_string('D');
1202
 
        }elsif($hchar eq $GAP_SYMBOL){ # Insertion
 
1356
        }elsif($hchar eq $self->{GAP_SYMBOL}){ # Insertion
1203
1357
            $cigar_string .= $self->_sub_cigar_string('I');
1204
1358
        }else{
1205
1359
            $self->throw("Impossible state that 2 gaps on each seq aligned");
1227
1381
    return $sub_cigar_string;
1228
1382
}
1229
1383
 
1230
 
 
1231
 
 
1232
1384
# needed before seqfeatures can be made
1233
1385
sub _pre_seq_feature {
1234
1386
    my $self = shift;
1235
1387
    my $algo = $self->{ALGORITHM};
1236
1388
 
1237
1389
    my ($queryfactor, $hitfactor) = (0,0);
1238
 
    if( $algo =~ /^(PSI)?T(BLAST|FAST|SW)[NY]/oi ) {
 
1390
    my ($hitmap, $querymap) = (1,1);
 
1391
    if( $algo =~ /^(?:PSI)?T(?:BLAST|FAST|SW)[NY]/oi ) {
1239
1392
        $hitfactor = 1;
 
1393
        $hitmap = 3;
1240
1394
    }
1241
 
    elsif ($algo =~ /^(FAST|BLAST)(X|Y|XY)/oi || $algo =~ /^P?GENEWISE/oi ) {
 
1395
    elsif ($algo =~ /^(?:FAST|BLAST)(?:X|Y|XY)/oi || $algo =~ /^P?GENEWISE/oi ) {
1242
1396
        $queryfactor = 1;
 
1397
        $querymap = 3;
1243
1398
    }
1244
1399
    elsif ($algo =~ /^T(BLAST|FAST|SW)(X|Y|XY)/oi || $algo =~ /^(BLAST|FAST|SW)N/oi || $algo =~ /^WABA|AXT|BLAT|BLASTZ|PSL|MEGABLAST|EXONERATE|SW|SMITH\-WATERMAN|SIM4$/){
 
1400
        if ($2) {
 
1401
            $hitmap = $querymap = 3;
 
1402
        }
1245
1403
        $hitfactor = 1;
1246
1404
        $queryfactor = 1;
1247
1405
    }
1248
1406
    elsif ($algo =~ /^RPS-BLAST/) {
1249
 
        $queryfactor = ($algo =~ /^RPS-BLAST\(BLASTX\)/) ? 1 : 0;
 
1407
        if ($algo =~ /^RPS-BLAST\(BLASTX\)/) {
 
1408
            $queryfactor = 1;
 
1409
            $querymap  = 3;
 
1410
        }
1250
1411
        $hitfactor = 0;
1251
1412
    }
 
1413
    else {
 
1414
        my $stranded = lc substr($self->{STRANDED}, 0,1);
 
1415
        $queryfactor = ($stranded eq 'q' || $stranded eq 'b') ? 1 : 0;
 
1416
        $hitfactor = ($stranded eq 'h' || $stranded eq 's' || $stranded eq 'b') ? 1 : 0;
 
1417
    }
1252
1418
    $self->{_query_factor} = $queryfactor;
1253
1419
    $self->{_hit_factor} = $hitfactor;
 
1420
    $self->{_hit_mapping} = $hitmap;
 
1421
    $self->{_query_mapping} = $querymap;
1254
1422
}
1255
1423
 
1256
1424
# make query seq feature
1301
1469
    $sim1->seqlength($self->{QUERY_LENGTH});
1302
1470
    $sim1->source_tag($self->{ALGORITHM});
1303
1471
    $sim1->seqdesc($self->{QUERY_DESC});
1304
 
    
1305
 
    $self->SUPER::feature1($sim1);
1306
 
 
 
1472
    $sim1->add_tag_value('meta', $self->{META}) if $self->can('meta');
1307
1473
    # to determine frame from something like FASTXY which doesn't
1308
1474
    # report the frame
1309
1475
    my $qframe = $self->{QUERY_FRAME};
1310
 
    if (defined $strand && ! defined $qframe && $queryfactor) {
1311
 
        $qframe = ( $self->query->start % 3 ) * $strand;
1312
 
    }
1313
 
    elsif (! defined $strand) {
1314
 
        $qframe = 0;
1315
 
    }
1316
 
    $self->{QUERY_FRAME} = $qframe;
 
1476
    
 
1477
    if (defined $strand && !defined $qframe && $queryfactor) {
 
1478
        $qframe = ( $qs % 3 ) * $strand;
 
1479
    }
 
1480
    elsif (!defined $strand) {
 
1481
        $qframe = 0;
 
1482
    }
 
1483
    
 
1484
    if( $qframe =~ /^([+-])?([0-3])/ ) {
 
1485
        my $dir = $1 || '+';
 
1486
        if($qframe && (($dir eq '-' && $strand >= 0) || ($dir eq '+' && $strand <= 0)) ) {
 
1487
            $self->warn("Query frame ($qframe) did not match strand of query ($strand)");
 
1488
        }
 
1489
        $qframe = $2 != 0 ? $2 - 1 : $2;
 
1490
    }  else {
 
1491
        $self->warn("Unknown query frame ($qframe)");
 
1492
        $qframe = 0;
 
1493
    }
 
1494
    
 
1495
    $sim1->frame($qframe);
 
1496
    $self->SUPER::feature1($sim1);
1317
1497
 
1318
1498
    $self->{_created_qff} = 1;
1319
1499
    $self->{_making_qff} = 0;
1320
 
    $self->_pre_frame;
1321
1500
}
1322
1501
 
1323
1502
# make subject seq feature
1363
1542
    $sim2->seqlength($self->{HIT_LENGTH});
1364
1543
    $sim2->source_tag($self->{ALGORITHM});
1365
1544
    $sim2->seqdesc($self->{HIT_DESC});
1366
 
    $self->SUPER::feature2($sim2);
1367
 
 
 
1545
    $sim2->add_tag_value('meta', $self->{META}) if $self->can('meta');
1368
1546
    my $hframe = $self->{HIT_FRAME};
1369
 
    if (defined $strand && ! defined $hframe && $hitfactor) {
 
1547
    
 
1548
    if (defined $strand && !defined $hframe && $hitfactor) {
1370
1549
        $hframe = ( $hs % 3 ) * $strand;
1371
1550
    }
1372
 
    elsif (! defined $strand) {
1373
 
        $hframe = 0;
1374
 
    }
1375
 
    $self->{HIT_FRAME} = $hframe;
 
1551
    elsif (!defined $strand) {
 
1552
        $hframe = 0;
 
1553
    }
 
1554
    
 
1555
    if( $hframe =~ /^([+-])?([0-3])/ ) {
 
1556
        my $dir = $1 || '+';
 
1557
        if($hframe && (($dir eq '-' && $strand >= 0) || ($dir eq '+' && $strand <= 0)) ) {
 
1558
            $self->warn("Subject frame ($hframe) did not match strand of subject ($strand)");
 
1559
        }
 
1560
        $hframe = $2 != 0 ? $2 - 1 : $2;
 
1561
    }  else {
 
1562
        $self->warn("Unknown subject frame ($hframe)");
 
1563
        $hframe = 0;
 
1564
    }
 
1565
    
 
1566
    $sim2->frame($hframe);
 
1567
    $self->SUPER::feature2($sim2);
1376
1568
 
1377
1569
    $self->{_created_sff} = 1;
1378
1570
    $self->{_making_sff} = 0;
1379
 
    $self->_pre_frame;
1380
 
}
1381
 
 
1382
 
# know the frame following seq feature creation
1383
 
sub _pre_frame {
1384
 
    my $self = shift;
1385
 
    $self->{_created_qff} || $self->_query_seq_feature;
1386
 
    $self->{_created_sff} || $self->_subject_seq_feature;
1387
 
    $self->{_did_preframe} = 1;
1388
 
    $self->frame;
1389
1571
}
1390
1572
 
1391
1573
# before calling the num_* methods
1397
1579
 
1398
1580
    if (! defined $identical) {
1399
1581
        if (! defined $percent_id) {
1400
 
            $self->warn("Did not defined the number of identical matches or overall percent identity in the HSP assuming 0");
 
1582
            $self->warn("Did not defined the number of identical matches or overall percent identity in the HSP; assuming 0");
1401
1583
            $identical = 0;
1402
1584
        }
1403
1585
        else {
1404
 
            $identical = int($percent_id * $self->{HSP_LENGTH});
 
1586
            $identical = sprintf("%.0f",$percent_id * $self->{HSP_LENGTH});
1405
1587
        }
1406
1588
    }
1407
1589
 
1408
1590
    if (! defined $conserved) {
1409
 
        $self->warn("Did not defined the number of conserved matches in the HSP assuming conserved == identical ($identical)")
 
1591
        $self->warn("Did not define the number of conserved matches in the HSP; assuming conserved == identical ($identical)")
1410
1592
            if( $self->{ALGORITHM} !~ /^((FAST|BLAST)N)|EXONERATE|SIM4|AXT|PSL|BLAT|BLASTZ|WABA/oi);
1411
1593
        $conserved = $identical;
1412
1594
    }
1416
1598
}
1417
1599
 
1418
1600
# before calling the frac_* methods
 
1601
 
1419
1602
sub _pre_frac {
1420
1603
    my $self = shift;
1421
1604
    my $hsp_len = $self->{HSP_LENGTH};
1446
1629
}
1447
1630
 
1448
1631
# before calling gaps()
 
1632
# This relies first on passed parameters (parser-dependent), then on gaps
 
1633
# calculated by seq_inds() (if set), then falls back to directly checking
 
1634
# for '-' as a last resort  
 
1635
 
1449
1636
sub _pre_gaps {
1450
1637
    my $self = shift;
1451
1638
    my $query_gaps = $self->{QUERY_GAPS};
1458
1645
    if( defined $query_gaps ) {
1459
1646
        $self->gaps('query', $query_gaps);
1460
1647
    } elsif( defined $query_seq ) {
1461
 
        $self->gaps('query', scalar ( $query_seq =~ tr/\-//));
 
1648
        my $qg = (defined $self->{'_query_offset'}) ? $self->seq_inds('query','gaps') : scalar( $query_seq =~ tr/\-//);
 
1649
        my $offset = $self->{'_query_offset'} || 1;
 
1650
        $self->gaps('query', $qg/$offset);
1462
1651
    }
1463
1652
    if( defined $hit_gaps ) {
1464
1653
        $self->gaps('hit', $hit_gaps);
1465
1654
    } elsif( defined $hit_seq ) {
1466
 
        $self->gaps('hit', scalar ( $hit_seq =~ tr/\-//));
 
1655
        my $hg = (defined $self->{'_sbjct_offset'}) ? $self->seq_inds('hit','gaps') : scalar( $hit_seq =~ tr/\-//);
 
1656
        my $offset = $self->{'_sbjct_offset'} || 1;
 
1657
        $self->gaps('hit', $hg/$offset);
1467
1658
    }
1468
1659
    if( ! defined $gaps ) {
1469
1660
        $gaps = $self->gaps("query") + $self->gaps("hit");