Bio/Tools/Genomewise.pm

   1 #
   2 # BioPerl module for Bio::Tools::Genomewise
   3 #
   4 # Copyright Jason Stajich <jason-at-bioperl.org>
   5 #
   6 # You may distribute this module under the same terms as perl itself
   7 #
   8 # POD documentation - main docs before the code
   9
  10 =head1 NAME
  11
  12 Bio::Tools::Genomewise - Results of one Genomewise run
  13
  14 =head1 SYNOPSIS
  15
  16   use Bio::Tools::Genomewise;
  17   my $gw = Bio::Tools::Genomewise(-file=>"genomewise.out");
  18
  19   while (my $gene = $gw->next_prediction){
  20       my @transcripts = $gene->transcripts;
  21       foreach my $t(@transcripts){
  22         my @exons =  $t->exons;
  23         foreach my $e(@exons){
  24             print $e->start." ".$e->end."\n";
  25         }
  26       }
  27   }
  28
  29 =head1 DESCRIPTION
  30
  31 This is the parser for the output of Genewise. It takes either a file
  32 handle or a file name and returns a
  33 Bio::SeqFeature::Gene::GeneStructure object.  You will need to specify
  34 the proper target sequence id on the object with the
  35 $feature-E<gt>seq_id($seqid).
  36
  37 =head1 FEEDBACK
  38
  39 =head2 Mailing Lists
  40
  41 User feedback is an integral part of the evolution of this and other
  42 Bioperl modules. Send your comments and suggestions preferably to one
  43 of the Bioperl mailing lists.  Your participation is much appreciated.
  44
  45   bioperl-l@bioperl.org                  - General discussion
  46   http://bioperl.org/wiki/Mailing_lists  - About the mailing lists
  47
  48 =head2 Support
  49
  50 Please direct usage questions or support issues to the mailing list:
  51
  52 I<bioperl-l@bioperl.org>
  53
  54 rather than to the module maintainer directly. Many experienced and
  55 reponsive experts will be able look at the problem and quickly
  56 address it. Please include a thorough description of the problem
  57 with code and data examples if at all possible.
  58
  59 =head2 Reporting Bugs
  60
  61 Report bugs to the Bioperl bug tracking system to help us keep track
  62 the bugs and their resolution.  Bug reports can be submitted via the
  63 web:
  64
  65   https://github.com/bioperl/bioperl-live/issues
  66
  67 =head1 AUTHOR - Fugu Team, Jason Stajich
  68
  69  Email: fugui-at-worf.fugu-sg.org
  70         jason-at-bioperl-dot-org
  71
  72 =head1 APPENDIX
  73
  74 The rest of the documentation details each of the object methods. Internal methods are usually preceded with a _
  75
  76 =cut
  77
  78
  79 # Let the code begin...
  80
  81
  82 package Bio::Tools::Genomewise;
  83 use vars qw($Srctag);
  84 use strict;
  85
  86 use Bio::Tools::AnalysisResult;
  87 use Bio::SeqFeature::Generic;
  88 use Bio::SeqFeature::Gene::Exon;
  89 use Bio::SeqFeature::FeaturePair;
  90 use Bio::SeqFeature::Gene::Transcript;
  91 use Bio::SeqFeature::Gene::GeneStructure;
  92
  93 use base qw(Bio::Tools::Genewise);
  94
  95 $Srctag = 'genomewise';
  96
  97 =head2 new
  98
  99  Title   : new
 100  Usage   : $obj->new(-file=>"genewise.out");
 101            $obj->new(-fh=>\*GW);
 102  Function: Constructor for genomewise wrapper. Takes either a file or filehandle
 103  Example :
 104  Returns : L<Bio::Tools::Genomewise>
 105
 106 =cut
 107
 108 sub new {
 109   my($class,@args) = @_;
 110   my $self = $class->SUPER::new(@args);
 111   return $self;
 112 }
 113
 114 =head2 _get_strand
 115
 116  Title   : _get_strand
 117  Usage   : $obj->_get_strand
 118  Function: takes start and end values, swap them if start>end and returns end
 119  Example :
 120  Returns :$start,$end,$strand
 121
 122 =cut
 123
 124 =head2 score
 125
 126  Title   : score
 127  Usage   : $obj->score
 128  Function: get/set for score info
 129  Example :
 130  Returns : a score value
 131
 132 =cut
 133
 134 =head2 _prot_id
 135
 136  Title   : _prot_id
 137  Usage   : $obj->_prot_id
 138  Function: get/set for protein id
 139  Example :
 140  Returns :a protein id
 141
 142 =cut
 143
 144 =head2 _target_id
 145
 146  Title   : _target_id
 147  Usage   : $obj->_target_id
 148  Function: get/set for genomic sequence id
 149  Example :
 150  Returns :a target id
 151
 152 =cut
 153
 154
 155 =head2 next_prediction
 156
 157  Title   : next_prediction
 158  Usage   : while($gene = $genewise->next_prediction()) {
 159                   # do something
 160            }
 161  Function: Returns the gene structure prediction of the Genomewise result
 162            file. Call this method repeatedly until FALSE is returned.
 163
 164  Example :
 165  Returns : a Bio::SeqFeature::Gene::GeneStructure object
 166  Args    :
 167
 168 =cut
 169
 170
 171 sub next_prediction {
 172     my ($self) = @_;
 173
 174     my $genes;
 175     while ($_ = $self->_readline) {
 176         $self->debug( $_ );
 177         last if m{^//};
 178
 179         if( /^Gene\s+\d+\s*$/ ) {
 180             $genes = Bio::SeqFeature::Gene::GeneStructure->new
 181                 (-source => $Srctag,
 182                  -seq_id => $self->_target_id, # if this had been specified
 183                  );
 184             $_ = $self->_readline;
 185             $self->debug( $_ );
 186
 187             unless ( /^Gene\s+(\d+)\s+(\d+)\s*$/ ) {
 188                 $self->warn("Unparseable genomewise output");
 189                 last;
 190             }
 191             my $transcript = Bio::SeqFeature::Gene::Transcript->new
 192                 (-source => $Srctag,
 193                  -seq_id => $self->_target_id, # if this had been specified
 194                  -start  => $1,
 195                  -end    => $2,
 196                  );
 197             my $nbr = 1;
 198             while( $_ = $self->_readline ) {
 199                 $self->debug( $_ );
 200
 201                 unless( m/^\s+Exon\s+(\d+)\s+(\d+)\s+phase\s+(\d+)/ ){
 202                     $self->_pushback($_);
 203                     last;
 204                 }
 205                 my ($e_start,$e_end,$phase,$e_strand) = ($1,$2,$3);
 206
 207                 ($e_start,$e_end,$e_strand) = $self->_get_strand($e_start,
 208                                                                  $e_end);
 209                 $transcript->strand($e_strand) unless $transcript->strand != 0;
 210
 211                 my $exon = Bio::SeqFeature::Gene::Exon->new
 212                     (-seq_id=>$self->_target_id,
 213                      -source => $Srctag,
 214                      -start=>$e_start,
 215                      -end=>$e_end,
 216                      -frame => $phase,
 217                      -strand=>$e_strand);
 218                 $exon->add_tag_value("Exon",$nbr++);
 219                 $exon->add_tag_value('phase',$phase);
 220                 $transcript->add_exon($exon);
 221             }
 222             $genes->add_transcript($transcript);
 223             last; # only process a single gene at a time
 224         }
 225     }
 226     return $genes;
 227 }
 228
 229 1;