view release on metacpan or  search on metacpan

Bio/Search/HSP/GenericHSP.pm  view on Meta::CPAN

 Args    : 'hit' or 'subject' or 'sbjct' to retrieve the strand of the subject,
           'query' to retrieve the query strand (default)

=cut

sub strand {
    my ($self,$type) = @_;

    if( $type =~ /^q/i && defined $self->{'QUERY_STRAND'} ) {
        return $self->{'QUERY_STRAND'};
    } elsif( $type =~ /^(hit|subject|sbjct)/i && defined $self->{'HIT_STRAND'} ) {
        return $self->{'HIT_STRAND'};
    } 

    return $self->SUPER::strand($type)
}

=head2 score

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

=head2 bits

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

=head1 Private methods

=cut

=head2 _calculate_seq_positions

 Title   : _calculate_seq_positions
 Usage   : $self->_calculate_seq_positions
 Function: Internal function
 Returns :
 Args    :


=cut

sub _calculate_seq_positions {
    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;