lib/Bio/Tools/ESTScan.pm

   1 #
   2 # BioPerl module for Bio::Tools::ESTScan
   3 #
   4 # Please direct questions and support issues to <bioperl-l@bioperl.org>
   5 #
   6 # Cared for by Hilmar Lapp <hlapp@gmx.net>
   7 #
   8 # Copyright Hilmar Lapp
   9 #
  10 # You may distribute this module under the same terms as perl itself
  11
  12 # POD documentation - main docs before the code
  13
  14 =head1 NAME
  15
  16 Bio::Tools::ESTScan - Results of one ESTScan run
  17
  18 =head1 SYNOPSIS
  19
  20    $estscan = Bio::Tools::ESTScan->new(-file => 'result.estscan');
  21    # filehandle:
  22    $estscan = Bio::Tools::ESTScan->new( -fh  => \*INPUT );
  23
  24    # parse the results
  25    # note: this class is-a Bio::Tools::AnalysisResult which implements
  26    # Bio::SeqAnalysisParserI, i.e., $genscan->next_feature() is the same
  27    while($gene = $estscan->next_prediction()) {
  28        # $gene is an instance of Bio::Tools::Prediction::Gene
  29        foreach my $orf ($gene->exons()) {
  30            # $orf is an instance of Bio::Tools::Prediction::Exon
  31            $cds_str = $orf->predicted_cds();
  32        }
  33    }
  34
  35    # essential if you gave a filename at initialization (otherwise the file
  36    # will stay open)
  37    $estscan->close();
  38
  39 =head1 DESCRIPTION
  40
  41 The ESTScan module provides a parser for ESTScan coding region prediction
  42 output.
  43
  44 This module inherits off L<Bio::Tools::AnalysisResult> and therefore
  45 implements the L<Bio::SeqAnalysisParserI> interface.
  46 See L<Bio::SeqAnalysisParserI>.
  47
  48 =head1 FEEDBACK
  49
  50 =head2 Mailing Lists
  51
  52 User feedback is an integral part of the evolution of this and other
  53 Bioperl modules. Send your comments and suggestions preferably to one
  54 of the Bioperl mailing lists.  Your participation is much appreciated.
  55
  56   bioperl-l@bioperl.org                  - General discussion
  57   http://bioperl.org/wiki/Mailing_lists  - About the mailing lists
  58
  59 =head2 Support
  60
  61 Please direct usage questions or support issues to the mailing list:
  62
  63 I<bioperl-l@bioperl.org>
  64
  65 rather than to the module maintainer directly. Many experienced and
  66 reponsive experts will be able look at the problem and quickly
  67 address it. Please include a thorough description of the problem
  68 with code and data examples if at all possible.
  69
  70 =head2 Reporting Bugs
  71
  72 Report bugs to the Bioperl bug tracking system to help us keep track
  73 the bugs and their resolution.  Bug reports can be submitted via the
  74 web:
  75
  76   https://github.com/bioperl/bioperl-live/issues
  77
  78 =head1 AUTHOR - Hilmar Lapp
  79
  80 Email hlapp@gmx.net (or hilmar.lapp@pharma.novartis.com)
  81
  82 =head1 APPENDIX
  83
  84 The rest of the documentation details each of the object methods. Internal methods are usually preceded with a _
  85
  86 =cut
  87
  88 # Let the code begin...
  89
  90 package Bio::Tools::ESTScan;
  91
  92 use strict;
  93 use Symbol;
  94
  95 use Bio::Root::Root;
  96 use Bio::Tools::Prediction::Exon;
  97
  98 use base qw(Bio::Tools::AnalysisResult);
  99
 100 sub _initialize_state {
 101     my ($self,@args) = @_;
 102
 103     # first call the inherited method!
 104     my $make = $self->SUPER::_initialize_state(@args);
 105
 106     if(! $self->analysis_method()) {
 107         $self->analysis_method('ESTScan');
 108     }
 109 }
 110
 111 =head2 analysis_method
 112
 113  Usage     : $estscan->analysis_method();
 114  Purpose   : Inherited method. Overridden to ensure that the name matches
 115              /estscan/i.
 116  Returns   : String
 117  Argument  : n/a
 118
 119 =cut
 120
 121 #-------------
 122 sub analysis_method {
 123 #-------------
 124     my ($self, $method) = @_;
 125     if($method && ($method !~ /estscan/i)) {
 126         $self->throw("method $method not supported in " . ref($self));
 127     }
 128     return $self->SUPER::analysis_method($method);
 129 }
 130
 131 =head2 next_feature
 132
 133  Title   : next_feature
 134  Usage   : while($orf = $estscan->next_feature()) {
 135                   # do something
 136            }
 137  Function: Returns the next gene structure prediction of the ESTScan result
 138            file. Call this method repeatedly until FALSE is returned.
 139
 140            The returned object is actually a SeqFeatureI implementing object.
 141            This method is required for classes implementing the
 142            SeqAnalysisParserI interface, and is merely an alias for
 143            next_prediction() at present.
 144
 145  Example :
 146  Returns : A Bio::Tools::Prediction::Gene object.
 147  Args    :
 148
 149 =cut
 150
 151 sub next_feature {
 152     my ($self,@args) = @_;
 153     # even though next_prediction doesn't expect any args (and this method
 154     # does neither), we pass on args in order to be prepared if this changes
 155     # ever
 156     return $self->next_prediction(@args);
 157 }
 158
 159 =head2 next_prediction
 160
 161  Title   : next_prediction
 162  Usage   : while($gene = $estscan->next_prediction()) {
 163                   # do something
 164            }
 165  Function: Returns the next gene structure prediction of the ESTScan result
 166            file. Call this method repeatedly until FALSE is returned.
 167
 168            So far, this method DOES NOT work for reverse strand predictions,
 169            even though the code looks like.
 170  Example :
 171  Returns : A Bio::Tools::Prediction::Gene object.
 172  Args    :
 173
 174 =cut
 175
 176 sub next_prediction {
 177     my ($self) = @_;
 178     my ($gene, $seq, $cds, $predobj);
 179     my $numins = 0;
 180
 181     # predictions are in the format of FASTA sequences and can be parsed one
 182     # at a time
 183     $seq = $self->_fasta_stream()->next_seq();
 184     return unless $seq;
 185     # there is a new prediction
 186     $gene = Bio::Tools::Prediction::Gene->new('-primary' => "ORFprediction",
 187                                               '-source' => "ESTScan");
 188     # score starts the description
 189     $seq->desc() =~ /^([\d.]+)\s*(.*)/ or
 190         $self->throw("unexpected format of description: no score in " .
 191                      $seq->desc());
 192     $gene->score($1);
 193     $seq->desc($2);
 194     # strand may end the description
 195     if($seq->desc() =~ /(.*)minus strand$/) {
 196         my $desc = $1;
 197         $desc =~ s/;\s+$//;
 198         $seq->desc($desc);
 199         $gene->strand(-1);
 200     } else {
 201         $gene->strand(1);
 202     }
 203     # check for the format: default or 'all-in-one' (option -a)
 204     if($seq->desc() =~ /^(\d+)\s+(\d+)\s+(\d+)\s+(\d+)\s*(.*)/) {
 205         # default format
 206         $seq->desc($5);
 207         $predobj = Bio::Tools::Prediction::Exon->new('-source' => "ESTScan",
 208                                                      '-start' => $3,
 209                                                      '-end' => $4);
 210         $predobj->strand($gene->strand());
 211         $predobj->score($gene->score()); # FIXME or $1, or $2 ?
 212         $predobj->primary_tag("InternalExon");
 213         $predobj->seq_id($seq->display_id());
 214         # add to gene structure object
 215         $gene->add_exon($predobj);
 216         # add predicted CDS
 217         $cds = $seq->seq();
 218         $cds =~ s/[a-z]//g; # remove the deletions, but keep the insertions
 219         $cds = Bio::PrimarySeq->new('-seq' => $cds,
 220                                     '-display_id' => $seq->display_id(),
 221                                     '-desc' => $seq->desc(),
 222                                     '-alphabet' => "dna");
 223         $gene->predicted_cds($cds);
 224         $predobj->predicted_cds($cds);
 225         if($gene->strand() == -1) {
 226             $self->warn("reverse strand ORF, but unable to reverse coordinates!");
 227         }
 228     } else {
 229         #
 230         # All-in-one format (hopefully). This encodes the following information
 231         # into the sequence:
 232         # 1) untranslated regions: stretches of lower-case letters
 233         # 2) translated regions: stretches of upper-case letters
 234         # 3) insertions in the translated regions: capital X
 235         # 4) deletions in the translated regions: a single lower-case letter
 236         #
 237         # if reverse strand ORF, save a lot of hassle by reversing the sequence
 238         if($gene->strand() == -1) {
 239             $seq = $seq->revcom();
 240         }
 241         my $seqstr = $seq->seq();
 242         while($seqstr =~ /^([a-z]*)([A-Z].*)$/) {
 243             # leading 5'UTR
 244             my $utr5 = $1;
 245             # exon + 3'UTR
 246             my $exonseq = $2;
 247             # strip 3'UTR and following exons
 248             if($exonseq =~ s/([a-z]{2,}.*)$//) {
 249                 $seqstr = $1;
 250             } else {
 251                 $seqstr = "";
 252             }
 253             # start: take care of yielding the absolute coordinate
 254             my $start = CORE::length($utr5) + 1;
 255             if($predobj) {
 256                 $start += $predobj->end() + $numins;
 257             }
 258             # for the end coordinate, we need to subtract the insertions
 259             $cds = $exonseq;
 260             $cds =~ s/[X]//g;
 261             my $end = $start + CORE::length($cds) - 1;
 262             # construct next exon object
 263             $predobj = Bio::Tools::Prediction::Exon->new('-start' => $start,
 264                                                          '-end' => $end);
 265             $predobj->source_tag("ESTScan");
 266             $predobj->primary_tag("InternalExon");
 267             $predobj->seq_id($seq->display_id());
 268             $predobj->strand($gene->strand());
 269             $predobj->score($gene->score());
 270             # add the exon to the gene structure object
 271             $gene->add_exon($predobj);
 272             # add the predicted CDS
 273             $cds = $exonseq;
 274             $cds =~ s/[a-z]//g; # remove the deletions, but keep the insertions
 275             $cds = Bio::PrimarySeq->new('-seq' => $cds,
 276                                         '-display_id' => $seq->display_id(),
 277                                         '-desc' => $seq->desc(),
 278                                         '-alphabet' => "dna");
 279             # only store the first one in the overall prediction
 280             $gene->predicted_cds($cds) unless $gene->predicted_cds();
 281             $predobj->predicted_cds($cds);
 282             # add the predicted insertions and deletions as subfeatures
 283             # of the exon
 284             my $fea = undef;
 285             while($exonseq =~ /([a-zX])/g) {
 286                 my $indel = $1;
 287                 # start and end: start looking at the position after the
 288                 # previous feature
 289                 if($fea) {
 290                     $start = $fea->start()+$numins;
 291                     $start -= 1 if($fea->primary_tag() eq 'insertion');
 292                 } else {
 293                     $start = $predobj->start()+$numins-1;
 294                 }
 295                 #print "# numins = $numins, indel = $indel, start = $start\n";
 296                 $start = index($seq->seq(), $indel, $start) + 1 - $numins;
 297                 $fea = Bio::SeqFeature::Generic->new('-start' => $start,
 298                                                      '-end' => $start);
 299                 $fea->source_tag("ESTScan");
 300                 $fea->seq_id($seq->display_id());
 301                 $fea->strand($predobj->strand());
 302                 if($indel eq 'X') {
 303                     # an insertion (depends on viewpoint: to get the 'real'
 304                     # CDS, a base has to be inserted, i.e., the HMMER model
 305                     # inserted a base; however, the sequencing process deleted
 306                     # a base that was there).
 307                     $fea->primary_tag("insertion");
 308                     # we need to count insertions because these are left out
 309                     # of any coordinates saved in the objects (which is correct
 310                     # because insertions change the original sequence, so
 311                     # coordinates wouldn't match)
 312                     $numins++;
 313                 } else {
 314                     # a deletion (depends on viewpoint: to get the 'real'
 315                     # CDS, a base has to be deleted, i.e., the HMMER model
 316                     # deleted a base; however, the sequencing process inserted
 317                     # a base that wasn't there).
 318                     $fea->primary_tag("deletion");
 319                     $fea->add_tag_value('base', $indel);
 320                 }
 321                 $predobj->add_sub_SeqFeature($fea);
 322             }
 323         }
 324     }
 325
 326     return $gene;
 327 }
 328
 329 =head2 close
 330
 331  Title   : close
 332  Usage   : $result->close()
 333  Function: Closes the file handle associated with this result file.
 334            Inherited method, overridden.
 335  Example :
 336  Returns :
 337  Args    :
 338
 339 =cut
 340
 341 sub close {
 342    my ($self, @args) = @_;
 343
 344    delete($self->{'_fastastream'});
 345    $self->SUPER::close(@args);
 346 }
 347
 348 =head2 _fasta_stream
 349
 350  Title   : _fasta_stream
 351  Usage   : $result->_fasta_stream()
 352  Function: Gets/Sets the FASTA sequence IO stream for reading the contents of
 353            the file associated with this MZEF result object.
 354
 355            If called for the first time, creates the stream from the filehandle
 356            if necessary.
 357  Example :
 358  Returns :
 359  Args    :
 360
 361 =cut
 362
 363 sub _fasta_stream {
 364     my ($self, $stream) = @_;
 365
 366     if($stream || (! exists($self->{'_fastastream'}))) {
 367         if(! $stream) {
 368             $stream = Bio::SeqIO->new('-fh' => $self->_fh(),
 369                                       '-format' => "fasta");
 370         }
 371         $self->{'_fastastream'} = $stream;
 372     }
 373     return $self->{'_fastastream'};
 374 }
 375
 376 1;
 377