lib/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
  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;