my $hsp = Bio::Search::HSP::GenericHSP->new( -algorithm => 'blastp',
                                                -evalue    => '1e-30',
                                                );

    $r_type = $hsp->algorithm;

    $pvalue = $hsp->p();

    $evalue = $hsp->evalue();

    $frac_id = $hsp->frac_identical( ['query'|'hit'|'total'] );

    $frac_cons = $hsp->frac_conserved( ['query'|'hit'|'total'] );

    $gaps = $hsp->gaps( ['query'|'hit'|'total'] );

    $qseq = $hsp->query_string;

    $hseq = $hsp->hit_string;

    $homo_string = $hsp->homology_string;

    $len = $hsp->length( ['query'|'hit'|'total'] );

    $len = $hsp->length( ['query'|'hit'|'total'] );

    $rank = $hsp->rank;

    my($class,%args) = @_;

    # don't pass anything to SUPER; complex hierarchy results in lots of work
    # for nothing
    
    my $self = $class->SUPER::new();

    # for speed, don't use _rearrange and just store all input data directly
    # with no method calls and no work done. work can be carried
    # out just-in-time later if desired
    while (my ($arg, $value) = each %args) {
        $arg =~ tr/a-z\055/A-Z/d;
        $self->{$arg} = $value;
    }
    my $bits = $self->{BITS};

    defined $self->{VERBOSE} && $self->verbose($self->{VERBOSE});
    if (exists $self->{GAP_SYMBOL}) {
        # not checking anything else but the length (must be 1 as only one gap
        # symbol allowed currently); can add in support for symbol checks or
        # multiple symbols later if needed
        $self->throw("Gap symbol must be of length 1") if
            CORE::length($self->{GAP_SYMBOL}) != 1;
    } else {
        # dafault
        $self->{GAP_SYMBOL} = '-';
    }
    $self->{ALGORITHM} ||= 'GENERIC';
    $self->{STRANDED} ||= 'NONE';
    
    if (! defined $self->{QUERY_LENGTH} || ! defined $self->{HIT_LENGTH}) {
        $self->throw("Must define hit and query length");
    }

    $self->{'_sequenceschanged'} = 1;
    
    $self->{_finished_new} = 1;
    return $self;

    my ($self, $type) = @_;
    if (!defined($self->{_sbjct_offset}) || !defined($self->{_query_offset})) {
        $self->_calculate_seq_offsets();
    }    
    my $key = $type =~ /sbjct|hit|tot/i ? 'sbjct' : 'query';
    
    my $offset = $self->{"_${key}_offset"};
    return $self->length($type) / $offset ;

    my ($self,$value) = @_;
    my $previous = $self->{'ALGORITHM'};
    if( defined $value || ! defined $previous ) {
        $value = $previous = '' unless defined $value;
        $self->{'ALGORITHM'} = $value;
    }

    return $previous;

    my ($self,$value) = @_;
    my $previous = $self->{'PVALUE'};
    if( defined $value  ) {
        $self->{'PVALUE'} = $value;
    }
    return $previous;

    my ($self,$value) = @_;
    my $previous = $self->{'EVALUE'};
    if( defined $value  ) {
        $self->{'EVALUE'} = $value;
    }
    return $previous;

   my ($self, $type,$value) = @_;

    unless ($self->{_did_prefrac}) {
        $self->_pre_frac;
    }

   $type = lc $type if defined $type;
   $type = 'hit' if( defined $type &&
		     $type =~ /subject|sbjct/);
   $type = 'total' if( ! defined $type || $type eq 'hsp' ||
                        $type !~ /query|hit|subject|sbjct|total/);
   my $previous = $self->{'_frac_identical'}->{$type};
   if( defined $value || ! defined $previous ) {
       $value = $previous = '' unless defined $value;
       if( $type eq 'hit' || $type eq 'query' ) {
           $self->$type()->frac_identical( $value);
       }
       $self->{'_frac_identical'}->{$type} = $value;
   }
   return $previous;

    my ($self, $type,$value) = @_;

    unless ($self->{_did_prefrac}) {
        $self->_pre_frac;
    }

    $type = lc $type if defined $type;
    $type = 'hit' if( defined $type && $type =~ /subject|sbjct/);
    $type = 'total' if( ! defined $type || $type eq 'hsp' ||
                        $type !~ /query|hit|subject|sbjct|total/);
    my $previous = $self->{'_frac_conserved'}->{$type};
    if( defined $value || ! defined $previous ) {
        $value = $previous = '' unless defined $value;
        $self->{'_frac_conserved'}->{$type} = $value;
    }
    return $previous;

    my ($self, $type, $value) = @_;

    unless ($self->{_did_pregaps}) {
        $self->_pre_gaps;
    }

    $type = lc $type if defined $type;
    $type = 'total' if( ! defined $type || $type eq 'hsp' ||
                        $type !~ /query|hit|subject|sbjct|total/);
    $type = 'hit' if $type =~ /sbjct|subject/;
    my $previous = $self->{'_gaps'}->{$type};
    if( defined $value || ! defined $previous ) {
        $value = $previous = '' unless defined $value;
        $self->{'_gaps'}->{$type} = $value;
    }
    return $previous || 0;

    my ($self,$value) = @_;
    my $previous = $self->{QUERY_SEQ};
    if( defined $value || ! defined $previous ) {
        $value = $previous = '' unless defined $value;
        $self->{QUERY_SEQ} = $value;
        # do some housekeeping so we know when to
        # re-run _calculate_seq_positions
        $self->{'_sequenceschanged'} = 1;
    }
    return $previous;

    my ($self,$value) = @_;
    my $previous = $self->{HIT_SEQ};
    if( defined $value || ! defined $previous ) {
        $value = $previous = '' unless defined $value;
        $self->{HIT_SEQ} = $value;
        # do some housekeeping so we know when to
        # re-run _calculate_seq_positions
        $self->{'_sequenceschanged'} = 1;
    }
    return $previous;

    my ($self,$value) = @_;
    my $previous = $self->{HOMOLOGY_SEQ};
    if( defined $value || ! defined $previous ) {
        $value = $previous = '' unless defined $value;
        $self->{HOMOLOGY_SEQ} = $value;
        # do some housekeeping so we know when to
        # re-run _calculate_seq_positions
        $self->{'_sequenceschanged'} = 1;
    }
    return $previous;

    my $self = shift;
    my $type = shift;

    $type = 'total' unless defined $type;
    $type = lc $type;

    if( $type =~ /^q/i ) {
        return $self->query()->length(shift);
    } elsif( $type =~ /^(hit|subject|sbjct)/ ) {
        return $self->hit()->length(shift);
    } else {
        my $v = shift;
        if( defined $v ) {
            $self->{HSP_LENGTH} = $v;
        }
        return $self->{HSP_LENGTH};
   }
    return 0; # should never get here

 return shift->length('hsp', shift);

    my $self = shift;

    unless ($self->{_did_prepi}) {
        $self->_pre_pi;
    }

    return $self->SUPER::percent_identity(@_);

    my $self = shift;
    my $val = shift;
    if (!defined $val) {
        # unfortunately, w/o args we need to warn about API changes
        $self->warn("API for frame() has changed.\n".
                    "Please refer to documentation for Bio::Search::HSP::GenericHSP;\n".
                    "returning query frame");
        $val = 'query';
    }
    $val =~ s/^\s+//;

    if( $val =~ /^q/i ) {
        return $self->query->frame(@_);
    } elsif( $val =~ /^hi|^s/i ) {
        return $self->hit->frame(@_);
    } elsif (  $val =~ /^list|array/i ) {
        return ($self->query->frame($_[0]), 
            $self->hit->frame($_[1]) );
    } elsif ( $val =~ /^\d+$/) {
        # old API i.e. frame($query_frame, $hit_frame). This catches all but one
        # case, where no arg is passed (so the hit is wanted).  
        $self->warn("API for frame() has changed.\n".
                    "Please refer to documentation for Bio::Search::HSP::GenericHSP");
        wantarray ? 
        return ($self->query->frame($val), 
            $self->hit->frame(@_) ) :
        return $self->hit->frame($val,@_);
    } else { 
        $self->warn("unrecognized component '$val' requested\n");
    }
    return 0;

    my ($self) = @_;
    require Bio::LocatableSeq;
    require Bio::SimpleAlign;
    my $aln = Bio::SimpleAlign->new();
    my $hs = $self->hit_string();
    my $qs = $self->query_string();
    # FASTA specific stuff moved to the FastaHSP object
    my $seqonly = $qs;
    $seqonly =~ s/[\-\s]//g;
    my ($q_nm,$s_nm) = ($self->query->seq_id(),
                        $self->hit->seq_id());
    # Should we silently change the name of the query or hit if it isn't
    # defined?  May need revisiting... cjfields 2008-12-3 (commented out below)
    
    #unless( defined $q_nm && CORE::length ($q_nm) ) {
    #    $q_nm = 'query';
    #}
    #unless( defined $s_nm && CORE::length ($s_nm) ) {
    #    $s_nm = 'hit';
    #}
    
    # mapping: 1 residues for every x coordinate positions
    my $query = Bio::LocatableSeq->new('-seq'   => $qs,
                                      '-id'    => $q_nm,
                                      '-start' => $self->query->start,
                                      '-end'   => $self->query->end,
                                      '-strand'   => $self->query->strand,
                                      '-force_nse' => $q_nm ? 0 : 1,
                                      '-mapping' => [1, $self->{_query_mapping}]
                                      );
    $seqonly = $hs;
    $seqonly =~ s/[\-\s]//g;
    my $hit =  Bio::LocatableSeq->new('-seq'    => $hs,
                                      '-id'    => $s_nm,
                                      '-start' => $self->hit->start,
                                      '-end'   => $self->hit->end,
                                      '-strand'   => $self->hit->strand,
                                      '-force_nse' => $s_nm ? 0 : 1,
                                      '-mapping' => [1, $self->{_hit_mapping}]
                                      );
    $aln->add_seq($query);
    $aln->add_seq($hit);
    return $aln;

    my ($self,$value) = @_;

    unless ($self->{_did_presimilar}) {
        $self->_pre_similar_stats;
    }

    if (defined $value) {
        $self->{CONSERVED} = $value;
    }
    return $self->{CONSERVED};

   my ($self,$value) = @_;

   unless ($self->{_did_presimilar}) {
        $self->_pre_similar_stats;
    }

   if( defined $value) {
       $self->{IDENTICAL} = $value;
   }
   return $self->{IDENTICAL};

    my ($self,$value) = @_;
    if( defined $value) {
        $self->{RANK} = $value;
    }
    return $self->{RANK};

   my ($self, $seqType, $class, $collapse) = @_;

   # prepare the internal structures - this is cached so
   # if the strings have not changed we're okay
   $self->_calculate_seq_positions();

   $seqType  ||= 'query';
   $class ||= 'identical';
   $collapse ||= 0;
   $seqType = 'sbjct' if $seqType eq 'hit';
   my $t = lc(substr($seqType,0,1));
   if( $t eq 'q' ) {
       $seqType = 'query';
   } elsif ( $t eq 's' || $t eq 'h' ) {
       $seqType = 'sbjct';
   } else {
       $self->warn("unknown seqtype $seqType using 'query'");
       $seqType = 'query';
   }
   
   $t = lc(substr($class,0,1));

   if( $t eq 'c' ) {
       if( $class =~ /conserved\-not\-identical/ ) {
	   $class = 'conserved';
       } else {
	   $class = 'conservedall';
       }
   } elsif( $t eq 'i' ) {
       $class = 'identical';
   } elsif( $t eq 'n' ) {
       $class = 'nomatch';
   } elsif( $t eq 'm' ) {
       $class = 'mismatch';
   } elsif( $t eq 'g' ) {
       $class = 'gap';
   } elsif( $t eq 'f' ) {
       $class = 'frameshift';    
   } else {
       $self->warn("unknown sequence class $class using 'identical'");
       $class = 'identical';
   }

   ## Sensitive to member name changes.
   $seqType  = "_\L$seqType\E";
   $class = "_\L$class\E";
   my @ary;

   if( $class eq '_gap' ) {
        # this means that we are remapping the gap length that is stored
        # in the hash (for example $self->{'_gapRes_query'} )
        # so we'll return an array which has the values of the position of the
        # of the gap (the key in the hash) + the gap length (value in the
        # hash for this key - 1.
        
        # changing this; since the index is the position prior to the insertion,
        # repeat the position based on the number of gaps inserted
        @ary = map { my @tmp;
                     # position holds number of gaps inserted
                     for my $g (1..$self->{seqinds}{"${class}Res$seqType"}->{$_}) {
                        push @tmp, $_ ;
                     }
                     @tmp}
              sort { $a <=> $b } keys %{ $self->{seqinds}{"${class}Res$seqType"}};
   } elsif( $class eq '_conservedall' ) {
       @ary = sort { $a <=> $b }
       keys %{ $self->{seqinds}{"_conservedRes$seqType"}},
       keys %{ $self->{seqinds}{"_identicalRes$seqType"}},
   }  else {
       @ary = sort { $a <=> $b } keys %{ $self->{seqinds}{"${class}Res$seqType"}};
   }
   require Bio::Search::BlastUtils if $collapse;

   return $collapse ? &Bio::Search::SearchUtils::collapse_nums(@ary) : @ary;

    my $self = shift;
    $self->_calculate_seq_positions();
    my $type = ($self->{_query_offset} == 3 && $self->{_sbjct_offset} == 3) ?
        'query/subject' :
        ($self->{_query_offset} == 3) ? 'query' :
        ($self->{_sbjct_offset} == 3) ? 'subject' : '';
    return $type;

    my $self = shift;
    unless ($self->{_created_qff}) {
        $self->_query_seq_feature;
    }
    return $self->SUPER::query(@_);

    my $self = shift;
    if (! $self->{_finished_new} || $self->{_making_qff}) {
        return $self->{_sim1} if $self->{_sim1};
        $self->{_sim1} = Bio::SeqFeature::Similarity->new();
        return $self->{_sim1};
    }
    unless ($self->{_created_qff}) {
        $self->_query_seq_feature;
    }
    return $self->SUPER::feature1(@_);

    my $self = shift;
    unless ($self->{_created_sff}) {
        $self->_subject_seq_feature;
    }
    return $self->SUPER::hit(@_);

    my $self = shift;
    if (! $self->{_finished_new} || $self->{_making_sff}) {
        return $self->{_sim2} if $self->{_sim2};
        $self->{_sim2} = Bio::SeqFeature::Similarity->new();
        return $self->{_sim2};
    }
    unless ($self->{_created_sff}) {
        $self->_subject_seq_feature;
    }
    return $self->SUPER::feature2(@_);

    my ($self, $val) = @_;
    if (!defined $self->{SIGNIFICANCE} || defined $val) {
        $self->{SIGNIFICANCE} = defined $val ? $val :
                                defined $self->evalue ?  $self->evalue :
                                defined $self->pvalue ? $$self->pvalue :
                                undef;        
        $self->query->significance($self->{SIGNIFICANCE});
    }
    return $self->{SIGNIFICANCE};

    my ($self,@args) = @_;
    return unless ( $self->{'_sequenceschanged'} );
    $self->{'_sequenceschanged'} = 0;
    my ($seqString, $qseq,$sseq) = ( $self->homology_string(),
                                     $self->query_string(),
                                     $self->hit_string() );
    my ($mlen, $qlen, $slen) = (CORE::length($seqString), CORE::length($qseq), CORE::length($sseq));
    my $qdir = $self->query->strand || 1;
    my $sdir = $self->hit->strand || 1;
    my ($resCount_query, $endpoint_query) = ($qdir <=0) ? ($self->query->end, $self->query->start)
        : ($self->query->start, $self->query->end);
    my ($resCount_sbjct, $endpoint_sbjct) = ($sdir <=0) ? ($self->hit->end, $self->hit->start)
        : ($self->hit->start, $self->hit->end);
    
    my $prog = $self->algorithm;
    
    if( $prog  =~ /FAST|SSEARCH|SMITH-WATERMAN/i ) {
    
        # we infer the end of the regional sequence where the first and last
        # non spaces are in the homology string
        # then we use the HSP->length to tell us how far to read
        # to cut off the end of the sequence
        
        my ($start, $rest) = (0,0);
        if( $seqString =~ /^(\s+)?(.*?)\s*$/ ) {
            ($start, $rest) = ($1 ? CORE::length($1) : 0, CORE::length($2));
        }
    
        $seqString = substr($seqString, $start, $rest);
        $qseq = substr($qseq, $start, $rest);
        $sseq = substr($sseq, $start, $rest);

        # commented out 10/29/08
        # removing frameshift symbols doesn't take into account the following:
        # 1) does not remove the same point in the homology string (get
        # positional errors)
        # 2) adjustments to the overall position in the string due to the
        # frameshift must be taken into consideration (get balancing errors)
        #
        # Frameshifts will be handled directly in the main loop. 
        # --chris
        
        # fasta reports some extra 'regional' sequence information
        # we need to clear out first
        # this seemed a bit insane to me at first, but it appears to
        # work --jason
        
        #$qseq =~ s![\\\/]!!g;
        #$sseq =~ s![\\\/]!!g;
    }

    if (!defined($self->{_sbjct_offset}) || !defined($self->{_query_offset})) {
        $self->_calculate_seq_offsets();
    }
    
    my ($qfs, $sfs) = (0,0);
    CHAR_LOOP:
    for my $pos (0..CORE::length($seqString)-1) {
        my @qrange = (0 - $qfs)..($self->{_query_offset} - 1);
        my @srange = (0 - $sfs)..($self->{_sbjct_offset} - 1);
        # $self->debug("QRange:@qrange SRange:@srange\n") if ($qfs || $sfs);
        ($qfs, $sfs) = (0,0);
        my ($mchar, $qchar, $schar) = (
            unpack("x$pos A1",$seqString) || ' ',
            $pos < CORE::length($qseq) ? unpack("x$pos A1",$qseq) : '-',
            $pos < CORE::length($sseq) ? unpack("x$pos A1",$sseq) : '-'
            );
        my $matchtype = '';
        my ($qgap, $sgap) = (0,0);
        if( $mchar eq '+' || $mchar eq '.') {    # conserved
            $self->{seqinds}{_conservedRes_query}{ $resCount_query + ($_ * $qdir) } = 1 for @qrange;
            $self->{seqinds}{_conservedRes_sbjct}{ $resCount_sbjct + ($_ * $sdir) } = 1 for @srange;
            $matchtype = 'conserved';
        } elsif( $mchar eq ':' || $mchar ne ' ' ) { # identical
            $self->{seqinds}{_identicalRes_query}{ $resCount_query + ($_ * $qdir) } = 1 for @qrange;
            $self->{seqinds}{_identicalRes_sbjct}{ $resCount_sbjct + ($_ * $sdir) } = 1 for @srange;
            $matchtype = 'identical';
        } elsif( $mchar eq ' ' ) {  # possible mismatch/nomatch/frameshift
            $qfs = $qchar eq '/'  ?  -1 : # base inserted to match frame
                   $qchar eq '\\' ?   1 : # base deleted to match frame
                   0;
            $sfs = $schar eq '/'  ?  -1 :
                   $schar eq '\\' ?   1 :
                   0;
            if ($qfs) {
                # Frameshifts are tricky.
                
                # '/' indicates that the next residue is shifted back one
                # (-1) frame position and is a deletion of one base (so to
                # correctly align, a base is inserted). That residue should only
                # occupy two sequence positions instead of three.
                
                # '\' indicates that the next residue is shifted forward
                # one (+1) frame position and is an insertion of one base (so to
                # correctly align, a base is removed). That residue should
                # occupy four sequence positions instead of three.
                
                # Note that gaps are not counted across from frameshifts.
                # Frameshift indices are a range of positions starting in the
                # previous sequence position in which the frameshift occurs;
                # they are ambiguous by nature.
                $self->{seqinds}{_frameshiftRes_query}{ $resCount_query - ($_ * $qdir * $qfs) } = $qfs for @qrange;
                $matchtype = "$qfs frameshift-query";
                $sgap = $qgap = 1;
            }
            elsif ($sfs) {
                $self->{seqinds}{_frameshiftRes_sbjct}{ $resCount_sbjct - ($_ * $sdir * $sfs) } = $sfs for @srange;
                $matchtype = "$sfs frameshift-sbcjt";
                $sgap = $qgap = 1;
            }
            elsif ($qchar eq $self->{GAP_SYMBOL}) {
                # gap are counted as being in the immediately preceeding residue
                # position; for translated sequences this is not the start of
                # the previous codon, but the third codon position
                $self->{seqinds}{_gapRes_query}{ $resCount_query - $qdir }++ for @qrange;
                $matchtype = 'gap-query';
                $qgap++;
            }
            elsif ($schar eq $self->{GAP_SYMBOL}) {
                $self->{seqinds}{_gapRes_sbjct}{ $resCount_sbjct - $sdir }++ for @srange;
                $matchtype = 'gap-sbjct';
                $sgap++;
            }
            else {
                # if not a gap or frameshift in either seq, must be mismatch
                $self->{seqinds}{_mismatchRes_query}{ $resCount_query + ($_ * $qdir) } = 1 for @qrange;
                $self->{seqinds}{_mismatchRes_sbjct}{ $resCount_sbjct + ($_ * $sdir) } = 1 for @srange;
                $matchtype = 'mismatch';
            }
            # always add a nomatch unless the current position in the seq is a gap
            if (!$sgap) {
                $self->{seqinds}{_nomatchRes_sbjct}{ $resCount_sbjct + ($_ * $sdir) } = 1 for @srange;
            }
            if (!$qgap) {
                $self->{seqinds}{_nomatchRes_query}{ $resCount_query + ($_ * $qdir) } = 1 for @qrange;
            }
        } else {
            $self->warn("Unknown midline character: [$mchar]");
        }
        # leave in and uncomment for future debugging
        #$self->debug(sprintf("%7d %1s[%1s]%1s %-7d Type: %-20s QOff:%-2d SOff:%-2d\n",
        #                     $resCount_query,
        #                     $qchar,
        #                     $mchar,
        #                     $schar,
        #                     $resCount_sbjct,
        #                     $matchtype,
        #                     ($self->{_query_offset} * $qdir),
        #                     ($self->{_sbjct_offset} * $sdir)));
        $resCount_query += ($qdir * (scalar(@qrange) + $qfs)) if !$qgap;
        $resCount_sbjct += ($sdir * (scalar(@srange) + $sfs)) if !$sgap;
    }
    return 1;

    my $self = shift;
    my $prog = $self->algorithm;
    ($self->{_sbjct_offset}, $self->{_query_offset}) = (1,1);
    if($prog =~ /^(?:PSI)?T(BLAST|FAST)(N|X|Y)/oi ) {
        $self->{_sbjct_offset} = 3;
        if ($1 eq 'BLAST' && $2 eq 'X') { #TBLASTX
            $self->{_query_offset} = 3;
        } 
    } elsif($prog =~ /^(BLAST|FAST)(X|Y|XY)/oi  ) {
        $self->{_query_offset} = 3;
    }
    1;

    my $self = shift;
    if(@_) { $self->{'N'} = shift; }
    # note that returning 1 is completely an assumption
    defined $self->{'N'} ? $self->{'N'} : 1;

    my ($self, $seqType) = @_;

    $seqType ||= 'query';
    $seqType = 'sbjct' if $seqType eq 'hit';

    my ($start, $end);
    if( $seqType eq 'query' ) {
        $start = $self->query->start;
        $end = $self->query->end;
    }
    else {
        $start = $self->hit->start;
        $end = $self->hit->end;
    }
    return ($start, $end);

    my $self = shift;

    return $self->{LINKS} = shift if @_;
    return $self->{LINKS};

    my $self = shift;

    return $self->{HSP_GROUP} = shift if @_;
    return $self->{HSP_GROUP};

    my $self = shift;

    return $self->{HIT_FEATURES} = shift if @_;
    return $self->{HIT_FEATURES};
}

# The cigar string code is written by Juguang Xiao <juguang@fugu-sg.org>

    my ($self, $arg) = @_;
    $self->warn("this is not a setter") if(defined $arg);

    unless(defined $self->{_cigar_string}){ # generate cigar string
        my $cigar_string = $self->generate_cigar_string($self->query_string, $self->hit_string);
        $self->{_cigar_string} = $cigar_string;
    } # end of unless

    return $self->{_cigar_string};

    my ($self, $qstr, $hstr) = @_;
    my @qchars = split //, $qstr;
    my @hchars = split //, $hstr;

    unless(scalar(@qchars) == scalar(@hchars)){
        $self->throw("two sequences are not equal in lengths");
    }

    $self->{_count_for_cigar_string} = 0;
    $self->{_state_for_cigar_string} = 'M';

    my $cigar_string = '';
    for(my $i=0; $i <= $#qchars; $i++){
        my $qchar = $qchars[$i];
        my $hchar = $hchars[$i];
        if($qchar ne $self->{GAP_SYMBOL} && $hchar ne $self->{GAP_SYMBOL}){ # Match
            $cigar_string .= $self->_sub_cigar_string('M');
        }elsif($qchar eq $self->{GAP_SYMBOL}){ # Deletion
            $cigar_string .= $self->_sub_cigar_string('D');
        }elsif($hchar eq $self->{GAP_SYMBOL}){ # Insertion
            $cigar_string .= $self->_sub_cigar_string('I');
        }else{
            $self->throw("Impossible state that 2 gaps on each seq aligned");
        }
    }
    $cigar_string .= $self->_sub_cigar_string('X'); # not forget the tail.
    return $cigar_string;
}

# an internal method to help generate cigar string

    my ($self, $new_state) = @_;

    my $sub_cigar_string = '';
    if($self->{_state_for_cigar_string} eq $new_state){
        $self->{_count_for_cigar_string} += 1; # Remain the state and increase the counter
    }else{
        $sub_cigar_string .= $self->{_count_for_cigar_string}
            unless $self->{_count_for_cigar_string} == 1;
        $sub_cigar_string .= $self->{_state_for_cigar_string};
        $self->{_count_for_cigar_string} = 1;
        $self->{_state_for_cigar_string} = $new_state;
    }
    return $sub_cigar_string;
}

# needed before seqfeatures can be made

    my $self = shift;
    my $algo = $self->{ALGORITHM};

    my ($queryfactor, $hitfactor) = (0,0);
    my ($hitmap, $querymap) = (1,1);
    if( $algo =~ /^(?:PSI)?T(?:BLAST|FAST|SW)[NY]/oi ) {
        $hitfactor = 1;
        $hitmap = 3;
    }
    elsif ($algo =~ /^(?:FAST|BLAST)(?:X|Y|XY)/oi || $algo =~ /^P?GENEWISE/oi ) {
        $queryfactor = 1;
        $querymap = 3;
    }
    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|SSEARCH|SMITH\-WATERMAN|SIM4$/){
        if ($2) {
            $hitmap = $querymap = 3;
        }
        $hitfactor = 1;
        $queryfactor = 1;
    }
    elsif ($algo =~ /^RPS-BLAST/) {
        if ($algo =~ /^RPS-BLAST\(BLASTX\)/) {
            $queryfactor = 1;
            $querymap  = 3;
        }
        $hitfactor = 0;
    }
    else {
        my $stranded = lc substr($self->{STRANDED}, 0,1);
        $queryfactor = ($stranded eq 'q' || $stranded eq 'b') ? 1 : 0;
        $hitfactor = ($stranded eq 'h' || $stranded eq 's' || $stranded eq 'b') ? 1 : 0;
    }
    $self->{_query_factor} = $queryfactor;
    $self->{_hit_factor} = $hitfactor;
    $self->{_hit_mapping} = $hitmap;
    $self->{_query_mapping} = $querymap;
}

# make query seq feature

    my $self = shift;
    $self->{_making_qff} = 1;
    my $qs = $self->{QUERY_START};
    my $qe = $self->{QUERY_END};
    unless (defined $self->{_query_factor}) {
        $self->_pre_seq_feature;
    }
    my $queryfactor = $self->{_query_factor};

    unless( defined $qe && defined $qs ) { $self->throw("Did not specify a Query End or Query Begin"); }

    my $strand;
    if ($qe > $qs) {  # normal query: start < end
        if ($queryfactor) {
            $strand = 1;
        }
        else {
            $strand = undef;
        }
    }
    else {
        if ($queryfactor) {
            $strand = -1;
        }
        else {
            $strand = undef;
        }
        ($qs,$qe) = ($qe,$qs);
    }

    # Note: many of these data are not query- and hit-specific.
    # Only start, end, name, length are.
    # We could be more efficient by only storing this info once.
    # steve chervitz --- Sat Apr  5 00:55:07 2003

    my $sim1 = $self->{_sim1} || Bio::SeqFeature::Similarity->new();
    $sim1->start($qs);
    $sim1->end($qe);
    $sim1->significance($self->{EVALUE});
    $sim1->bits($self->{BITS});
    $sim1->score($self->{SCORE});
    $sim1->strand($strand);
    $sim1->seq_id($self->{QUERY_NAME});
    $sim1->seqlength($self->{QUERY_LENGTH});
    $sim1->source_tag($self->{ALGORITHM});
    $sim1->seqdesc($self->{QUERY_DESC});
    $sim1->add_tag_value('meta', $self->{META}) if $self->can('meta');
    # to determine frame from something like FASTXY which doesn't
    # report the frame
    my $qframe = $self->{QUERY_FRAME};
    
    if (defined $strand && !defined $qframe && $queryfactor) {
        $qframe = ( $qs % 3 ) * $strand;
    }
    elsif (!defined $strand) {
        $qframe = 0;
    }
    
    if( $qframe =~ /^([+-])?([0-3])/ ) {
        my $dir = $1 || '+';
        if($qframe && (($dir eq '-' && $strand >= 0) || ($dir eq '+' && $strand <= 0)) ) {
            $self->warn("Query frame ($qframe) did not match strand of query ($strand)");
        }
        $qframe = $2 != 0 ? $2 - 1 : $2;
    }  else {
        $self->warn("Unknown query frame ($qframe)");
        $qframe = 0;
    }
    
    $sim1->frame($qframe);
    $self->SUPER::feature1($sim1);

    $self->{_created_qff} = 1;
    $self->{_making_qff} = 0;
}

# make subject seq feature

    my $self = shift;
    $self->{_making_sff} = 1;
    my $hs = $self->{HIT_START};
    my $he = $self->{HIT_END};
    unless (defined $self->{_hit_factor}) {
        $self->_pre_seq_feature;
    }
    my $hitfactor = $self->{_hit_factor};

    unless( defined $he && defined $hs ) { $self->throw("Did not specify a Hit End or Hit Begin"); }

    my $strand;
    if ($he > $hs) { # normal subject
        if ($hitfactor) {
            $strand = 1;
        }
        else {
            $strand = undef;
        }
    }
    else {
        if ($hitfactor) {
            $strand = -1;
        }
        else {
            $strand = undef;
        }
        ($hs,$he) = ( $he,$hs); # reverse subject: start bigger than end
    }

    my $sim2 = $self->{_sim2} || Bio::SeqFeature::Similarity->new();
    $sim2->start($hs);
    $sim2->end($he);
    $sim2->significance($self->{EVALUE});
    $sim2->bits($self->{BITS});
    $sim2->score($self->{SCORE});
    $sim2->strand($strand);
    $sim2->seq_id($self->{HIT_NAME});
    $sim2->seqlength($self->{HIT_LENGTH});
    $sim2->source_tag($self->{ALGORITHM});
    $sim2->seqdesc($self->{HIT_DESC});
    $sim2->add_tag_value('meta', $self->{META}) if $self->can('meta');
    my $hframe = $self->{HIT_FRAME};
    
    if (defined $strand && !defined $hframe && $hitfactor) {
        $hframe = ( $hs % 3 ) * $strand;
    }
    elsif (!defined $strand) {
        $hframe = 0;
    }
    
    if( $hframe =~ /^([+-])?([0-3])/ ) {
        my $dir = $1 || '+';
        if($hframe && (($dir eq '-' && $strand >= 0) || ($dir eq '+' && $strand <= 0)) ) {
            $self->warn("Subject frame ($hframe) did not match strand of subject ($strand)");
        }
        $hframe = $2 != 0 ? $2 - 1 : $2;
    }  else {
        $self->warn("Unknown subject frame ($hframe)");
        $hframe = 0;
    }
    
    $sim2->frame($hframe);
    $self->SUPER::feature2($sim2);

    $self->{_created_sff} = 1;
    $self->{_making_sff} = 0;
}

# before calling the num_* methods

    my $self = shift;
    my $identical = $self->{IDENTICAL};
    my $conserved = $self->{CONSERVED};
    my $percent_id = $self->{PERCENT_IDENTITY};

    if (! defined $identical) {
        if (! defined $percent_id) {
            $self->warn("Did not defined the number of identical matches or overall percent identity in the HSP; assuming 0");
            $identical = 0;
        }
        else {
            $identical = sprintf("%.0f",$percent_id * $self->{HSP_LENGTH});
        }
    }

    if (! defined $conserved) {
        $self->warn("Did not define the number of conserved matches in the HSP; assuming conserved == identical ($identical)")
            if( $self->{ALGORITHM} !~ /^((FAST|BLAST)N)|EXONERATE|SIM4|AXT|PSL|BLAT|BLASTZ|WABA/oi);
        $conserved = $identical;
    }
    $self->{IDENTICAL} = $identical;
    $self->{CONSERVED} = $conserved;
    $self->{_did_presimilar} = 1;
}

# before calling the frac_* methods

    my $self = shift;
    
    my $hsp_len = $self->{HSP_LENGTH};
    my $hit_len = $self->{HIT_LENGTH};
    my $query_len = $self->{QUERY_LENGTH};

    my $identical = $self->num_identical;
    my $conserved = $self->num_conserved;

    $self->{_did_prefrac} = 1;
    my $logical;
    if( $hsp_len ) {
        $self->length('total', $hsp_len);
        $logical = $self->_logical_length('total');
        $self->frac_identical( 'total', $identical / $hsp_len);
        $self->frac_conserved( 'total', $conserved / $hsp_len);
    }
    if( $hit_len ) {
        $logical = $self->_logical_length('hit');
        $self->frac_identical( 'hit', $identical / $logical);
        $self->frac_conserved( 'hit', $conserved / $logical);
    }
    if( $query_len ) {
        $logical = $self->_logical_length('query');
        $self->frac_identical( 'query', $identical / $logical) ;
        $self->frac_conserved( 'query', $conserved / $logical);
    }
}

# before calling gaps()
# This relies first on passed parameters (parser-dependent), then on gaps
# calculated by seq_inds() (if set), then falls back to directly checking
# for '-' as a last resort  

    my $self = shift;
    my $query_gaps = $self->{QUERY_GAPS};
    my $query_seq = $self->{QUERY_SEQ};
    my $hit_gaps = $self->{HIT_GAPS};
    my $hit_seq = $self->{HIT_SEQ};
    my $gaps = $self->{HSP_GAPS};

    $self->{_did_pregaps} = 1; # well, we're in the process; avoid recursion
    if( defined $query_gaps ) {
        $self->gaps('query', $query_gaps);
    } elsif( defined $query_seq ) {
        my $qg = (defined $self->{'_query_offset'}) ? $self->seq_inds('query','gaps') : scalar( $query_seq =~ tr/\-//);
        my $offset = $self->{'_query_offset'} || 1;
        $self->gaps('query', $qg/$offset);
    }
    if( defined $hit_gaps ) {
        $self->gaps('hit', $hit_gaps);
    } elsif( defined $hit_seq ) {
        my $hg = (defined $self->{'_sbjct_offset'}) ? $self->seq_inds('hit','gaps') : scalar( $hit_seq =~ tr/\-//);
        my $offset = $self->{'_sbjct_offset'} || 1;
        $self->gaps('hit', $hg/$offset);
    }
    if( ! defined $gaps ) {
        $gaps = $self->gaps("query") + $self->gaps("hit");
    }
    $self->gaps('total', $gaps);
}

# before percent_identity

    my $self = shift;
    $self->{_did_prepi} = 1;
    $self->percent_identity($self->{PERCENT_IDENTITY} || $self->frac_identical('total')*100) if( $self->{HSP_LENGTH} > 0 );
}

1;

  bioperl-l@bioperl.org                  - General discussion
  http://bioperl.org/wiki/Mailing_lists  - About the mailing lists

  https://redmine.open-bio.org/projects/bioperl/

 Title   : score
 Usage   : $score = $obj->score();
           $obj->score($value);
 Function: Get/Set the score value
 Returns : numeric
 Args    : [optional] new value to set

 Title   : bits
 Usage   : $bits = $obj->bits();
           $obj->bits($value);
 Function: Get/Set the bits value
 Returns : numeric
 Args    : [optional] new value to set

 Title   : strand
 Usage   : $hsp->strand('query')
 Function: Retrieves the strand for the HSP component requested
 Returns : +1 or -1 (0 if unknown)
 Args    : 'hit' or 'subject' or 'sbjct' to retrieve the strand of the subject
           'query' to retrieve the query strand (default)

            PG    P    G     GP   R      PLGP