lib/Bio/Search/HSP/FastaHSP.pm

   1 #
   2 # BioPerl module for Bio::Search::HSP::FastaHSP
   3 #
   4 # Please direct questions and support issues to <bioperl-l@bioperl.org>
   5 #
   6 # Cared for by Jason Stajich <jason@bioperl.org>
   7 #
   8 # Copyright Jason Stajich
   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::Search::HSP::FastaHSP - HSP object for FASTA specific data
  17
  18 =head1 SYNOPSIS
  19
  20   # get a FastaHSP from a SearchIO stream
  21   my $in = Bio::SearchIO->new(-format => 'fasta', -file => 'filename.fasta');
  22
  23   while( my $r = $in->next_result) {
  24       while( my $hit = $r->next_result ) {
  25            while( my $hsp = $hit->next_hsp ) {
  26               print "smith-waterman score (if available): ",
  27                     $hsp->sw_score(),"\n";
  28            }
  29       }
  30   }
  31
  32 =head1 DESCRIPTION
  33
  34 Describe the object here
  35
  36 =head1 FEEDBACK
  37
  38 =head2 Mailing Lists
  39
  40 User feedback is an integral part of the evolution of this and other
  41 Bioperl modules. Send your comments and suggestions preferably to
  42 the Bioperl mailing list.  Your participation is much appreciated.
  43
  44   bioperl-l@bioperl.org                  - General discussion
  45   http://bioperl.org/wiki/Mailing_lists  - About the mailing lists
  46
  47 =head2 Support
  48
  49 Please direct usage questions or support issues to the mailing list:
  50
  51 I<bioperl-l@bioperl.org>
  52
  53 rather than to the module maintainer directly. Many experienced and
  54 reponsive experts will be able look at the problem and quickly
  55 address it. Please include a thorough description of the problem
  56 with code and data examples if at all possible.
  57
  58 =head2 Reporting Bugs
  59
  60 Report bugs to the Bioperl bug tracking system to help us keep track
  61 of the bugs and their resolution. Bug reports can be submitted via the
  62 web:
  63
  64   https://github.com/bioperl/bioperl-live/issues
  65
  66 =head1 AUTHOR - Jason Stajich
  67
  68 Email jason-at-bioperl.org
  69
  70 =head1 APPENDIX
  71
  72 The rest of the documentation details each of the object methods.
  73 Internal methods are usually preceded with a _
  74
  75 =cut
  76
  77
  78 # Let the code begin...
  79
  80
  81 package Bio::Search::HSP::FastaHSP;
  82 use strict;
  83
  84
  85 use base qw(Bio::Search::HSP::GenericHSP);
  86
  87 =head2 new
  88
  89  Title   : new
  90  Usage   : my $obj = Bio::Search::HSP::FastaHSP->new();
  91  Function: Builds a new Bio::Search::HSP::FastaHSP object
  92  Returns : Bio::Search::HSP::FastaHSP
  93  Args    : -swscore => smith-waterman score
  94
  95 =cut
  96
  97 sub new {
  98   my($class,@args) = @_;
  99
 100   my $self = $class->SUPER::new(@args);
 101
 102   my ($swscore, $evalue2) = $self->_rearrange([qw(SWSCORE EVALUE2)], @args);
 103
 104   defined $swscore && $self->sw_score($swscore);
 105
 106   defined $evalue2 && $self->evalue2($evalue2);
 107
 108   return $self;
 109 }
 110
 111
 112 =head2 sw_score
 113
 114  Title   : sw_score
 115  Usage   : $obj->sw_score($newval)
 116  Function: Get/Set Smith-Waterman score
 117  Returns : value of sw_score
 118  Args    : newvalue (optional)
 119
 120
 121 =cut
 122
 123 sub sw_score{
 124     my ($self,$value) = @_;
 125     if( defined $value || ! defined $self->{'_sw_score'} ) {
 126         $value = 0 unless defined $value; # default value
 127         $self->{'_sw_score'} = $value;
 128     }
 129     return $self->{'_sw_score'};
 130 }
 131
 132 =head2 evalue2
 133
 134  Title   : evalue2
 135  Usage   : $obj->evalue2($newval)
 136  Function: Get/Set E2() expectation value
 137  Returns : value of evalue2
 138  Args    : newvalue (optional)
 139
 140
 141 =cut
 142
 143 sub evalue2{
 144     my ($self,$value) = @_;
 145     if( defined $value || ! defined $self->{'_evalue2'} ) {
 146         $value = 0 unless defined $value; # default value
 147         $self->{'_evalue2'} = $value;
 148     }
 149     return $self->{'_evalue2'};
 150 }
 151
 152
 153 sub get_aln {
 154     my ($self) = @_;
 155     require Bio::LocatableSeq;
 156     require Bio::SimpleAlign;
 157     my $aln = Bio::SimpleAlign->new();
 158     my $hs = $self->hit_string();
 159     my $qs = $self->query_string();
 160
 161     # fasta reports some extra 'regional' sequence information
 162     # we need to clear out first
 163     # this seemed a bit insane to me at first, but it appears to
 164     # work --jason
 165
 166     # modified to deal with LocatableSeq's end point verification and to deal
 167     # with frameshifts (which shift the end points in translated sequences).
 168
 169     # we infer the end of the regional sequence where the first
 170     # non space is in the homology string
 171     # then we use the HSP->length to tell us how far to read
 172     # to cut off the end of the sequence
 173
 174     my ($start, $rest) = (0, 0);
 175     if( $self->homology_string() =~ /^(\s+)?(.*?)\s*$/ ) {
 176       ($start, $rest) = ($1 ? CORE::length($1) : 0, CORE::length($2));
 177     }
 178     $self->debug("hs seq is '$hs'\n");
 179     $self->debug("qs seq is '$qs'\n");
 180
 181     $hs = substr($hs, $start,$rest);
 182     $qs = substr($qs, $start,$rest);
 183
 184     my $seqonly = $qs;
 185     $seqonly =~ s/\s+//g;
 186     my ($q_nm,$s_nm) = ($self->query->seq_id(),
 187                         $self->hit->seq_id());
 188     unless( defined $q_nm && CORE::length ($q_nm) ) {
 189         $q_nm = 'query';
 190     }
 191     unless( defined $s_nm && CORE::length ($s_nm) ) {
 192         $s_nm = 'hit';
 193     }
 194     $self->_calculate_seq_positions;
 195     my $query = Bio::LocatableSeq->new('-seq'   => $seqonly,
 196                                       '-id'    => $q_nm,
 197                                       '-start' => $self->query->start,
 198                                       '-end'   => $self->query->end,
 199                       '-frameshifts' => (exists $self->{seqinds}{_frameshiftRes_query}) ?
 200                             $self->{seqinds}{_frameshiftRes_query} : undef,
 201                       '-mapping' => [1, $self->{_query_mapping}],
 202                       -verbose => $self->verbose
 203                                       );
 204     $seqonly = $hs;
 205     $seqonly =~ s/\s+//g;
 206     my $hit =  Bio::LocatableSeq->new('-seq'    => $seqonly,
 207                                       '-id'    => $s_nm,
 208                                       '-start' => $self->hit->start,
 209                                       '-end'   => $self->hit->end,
 210                       '-frameshifts' => exists $self->{seqinds}{_frameshiftRes_sbjct} ?
 211                             $self->{seqinds}{_frameshiftRes_sbjct} : undef,
 212                       '-mapping' => [1, $self->{_hit_mapping}],
 213                       -verbose => $self->verbose
 214                                       );
 215     $aln->add_seq($query);
 216     $aln->add_seq($hit);
 217     return $aln;
 218 }
 219
 220
 221 1;