lib/Bio/Tools/Grail.pm

   1 #
   2 # BioPerl module for Bio::Tools::Grail
   3 #
   4 # Please direct questions and support issues to <bioperl-l@bioperl.org>
   5 #
   6 # Cared for by Jason Stajich <jason-at-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::Tools::Grail - Results of one Grail run
  17
  18 =head1 SYNOPSIS
  19
  20    $grail = Bio::Tools::Grail->new(-file => 'result.grail');
  21    # filehandle:
  22    $grail = Bio::Tools::Grail->new( -fh  => \*INPUT );
  23
  24    # parse the results
  25    while($gene = $grail->next_prediction()) {
  26        # $gene is an instance of Bio::Tools::Prediction::Gene
  27
  28        # $gene->exons() returns an array of
  29        # Bio::Tools::Prediction::Exon objects
  30        # all exons:
  31        @exon_arr = $gene->exons();
  32
  33        # initial exons only
  34        @init_exons = $gene->exons('Initial');
  35        # internal exons only
  36        @intrl_exons = $gene->exons('Internal');
  37        # terminal exons only
  38        @term_exons = $gene->exons('Terminal');
  39        # singleton exons only -- should be same as $gene->exons() because
  40        # there are no other exons supposed to exist in this structure
  41        @single_exons = $gene->exons('Single');
  42    }
  43
  44    # essential if you gave a filename at initialization (otherwise the file
  45    # will stay open)
  46    $genscan->close();
  47
  48 =head1 DESCRIPTION
  49
  50 The Grail module provides a parser for Grail gene structure prediction
  51 output.
  52
  53
  54 =head1 FEEDBACK
  55
  56 =head2 Mailing Lists
  57
  58 User feedback is an integral part of the evolution of this and other
  59 Bioperl modules. Send your comments and suggestions preferably to one
  60 of the Bioperl mailing lists.  Your participation is much appreciated.
  61
  62   bioperl-l@bioperl.org                  - General discussion
  63   http://bioperl.org/wiki/Mailing_lists  - About the mailing lists
  64
  65 =head2 Support
  66
  67 Please direct usage questions or support issues to the mailing list:
  68
  69 I<bioperl-l@bioperl.org>
  70
  71 rather than to the module maintainer directly. Many experienced and
  72 reponsive experts will be able look at the problem and quickly
  73 address it. Please include a thorough description of the problem
  74 with code and data examples if at all possible.
  75
  76 =head2 Reporting Bugs
  77
  78 Report bugs to the Bioperl bug tracking system to help us keep track
  79 the bugs and their resolution.  Bug reports can be submitted via the
  80 web:
  81
  82   https://github.com/bioperl/bioperl-live/issues
  83
  84 =head1 AUTHOR - Jason Stajich
  85
  86 Email jason-at-bioperl.org
  87
  88 =head1 APPENDIX
  89
  90 The rest of the documentation details each of the object
  91 methods. Internal methods are usually preceded with a _
  92
  93 =cut
  94
  95 # Let the code begin...
  96
  97 package Bio::Tools::Grail;
  98 use strict;
  99
 100 use Bio::Tools::Prediction::Gene;
 101 use Bio::Tools::Prediction::Exon;
 102 use Symbol;
 103
 104 use base qw(Bio::Root::IO Bio::Root::Root);
 105
 106 sub new {
 107   my($class,@args) = @_;
 108
 109   my $self = $class->SUPER::new(@args);
 110   $self->_initialize_io(@args);
 111
 112   return $self;
 113 }
 114
 115 =head2 next_prediction
 116
 117  Title   : next_prediction
 118  Usage   : while($gene = $grail->next_prediction()) {
 119                   # do something
 120            }
 121  Function: Returns the next gene structure prediction of the Grail result
 122            file. Call this method repeatedly until FALSE is returned.
 123
 124  Example :
 125  Returns : A Bio::Tools::Prediction::Gene object.
 126  Args    :
 127
 128 =cut
 129
 130 sub next_prediction {
 131     my ($self) = @_;
 132
 133     # get next gene structure
 134     my $gene = $self->_prediction();
 135
 136     if($gene) {
 137         # fill in predicted protein, and if available the predicted CDS
 138         #
 139         my ($id, $seq);
 140         # use the seq stack if there's a seq on it
 141         my $seqobj = pop(@{$self->{'_seqstack'}});
 142         if(! $seqobj) {
 143             # otherwise read from input stream
 144             ($id, $seq) = $self->_read_fasta_seq();
 145             $seqobj = Bio::PrimarySeq->new('-seq' => $seq,
 146                                            '-display_id' => $id,
 147                                            '-alphabet' => "protein");
 148         }
 149         # check that prediction number matches the prediction number
 150         # indicated in the sequence id (there may be incomplete gene
 151         # predictions that contain only signals with no associated protein
 152         # and CDS, like promoters, poly-A sites etc)
 153         $gene->primary_tag() =~ /[^0-9]([0-9]+)$/;
 154         my $prednr = $1;
 155         if($seqobj->display_id() !~ /_predicted_\w+_$prednr\|/) {
 156             # this is not our sequence, so push back for the next prediction
 157             push(@{$self->{'_seqstack'}}, $seqobj);
 158         } else {
 159             $gene->predicted_protein($seqobj);
 160             # CDS prediction, too?
 161             if($self->_has_cds()) {
 162                 ($id, $seq) = $self->_read_fasta_seq();
 163                 $seqobj = Bio::PrimarySeq->new('-seq' => $seq,
 164                                                '-display_id' => $id,
 165                                                '-alphabet' => "dna");
 166                 $gene->predicted_cds($seqobj);
 167             }
 168         }
 169     }
 170     return $gene;
 171 }
 172
 173 =head2 _parse_predictions
 174
 175  Title   : _parse_predictions()
 176  Usage   : $obj->_parse_predictions()
 177  Function: Parses the prediction section. Automatically called by
 178            next_prediction() if not yet done.
 179  Example :
 180  Returns :
 181
 182 =cut
 183
 184 sub _parse_predictions {
 185     my ($self) = @_;
 186
 187     # code needs to go here
 188
 189     $self->_predictions_parsed(1);
 190 }
 191
 192 =head2 _prediction
 193
 194  Title   : _prediction()
 195  Usage   : $gene = $obj->_prediction()
 196  Function: internal
 197  Example :
 198  Returns :
 199
 200 =cut
 201
 202 sub _prediction {
 203     my ($self) = @_;
 204
 205     return unless(exists($self->{'_preds'}) && @{$self->{'_preds'}});
 206     return shift(@{$self->{'_preds'}});
 207 }
 208
 209 =head2 _add_prediction
 210
 211  Title   : _add_prediction()
 212  Usage   : $obj->_add_prediction($gene)
 213  Function: internal
 214  Example :
 215  Returns :
 216
 217 =cut
 218
 219 sub _add_prediction {
 220     my ($self, $gene) = @_;
 221
 222     if(! exists($self->{'_preds'})) {
 223         $self->{'_preds'} = [];
 224     }
 225     push(@{$self->{'_preds'}}, $gene);
 226 }
 227
 228 =head2 _predictions_parsed
 229
 230  Title   : _predictions_parsed
 231  Usage   : $obj->_predictions_parsed
 232  Function: internal
 233  Example :
 234  Returns : TRUE or FALSE
 235
 236 =cut
 237
 238 sub _predictions_parsed {
 239     my ($self, $val) = @_;
 240
 241     $self->{'_preds_parsed'} = $val if $val;
 242     if(! exists($self->{'_preds_parsed'})) {
 243         $self->{'_preds_parsed'} = 0;
 244     }
 245     return $self->{'_preds_parsed'};
 246 }
 247
 248 =head2 _has_cds
 249
 250  Title   : _has_cds()
 251  Usage   : $obj->_has_cds()
 252  Function: Whether or not the result contains the predicted CDSs, too.
 253  Example :
 254  Returns : TRUE or FALSE
 255
 256 =cut
 257
 258 sub _has_cds {
 259     my ($self, $val) = @_;
 260
 261     $self->{'_has_cds'} = $val if $val;
 262     if(! exists($self->{'_has_cds'})) {
 263         $self->{'_has_cds'} = 0;
 264     }
 265     return $self->{'_has_cds'};
 266 }
 267
 268 1;