Bio/Tools/Genomewise.pm

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