← Back to branch summary

~ubuntu-branches/ubuntu/raring/bioperl/raring

« back to all changes in this revision

Viewing changes to Bio/Tools/Genewise.pm

  • Committer: Bazaar Package Importer
  • Author(s): Charles Plessy
  • Date: 2008-03-18 14:44:57 UTC
  • mfrom: (4 hardy)
  • mto: This revision was merged to the branch mainline in revision 6.
  • Revision ID: james.westby@ubuntu.com-20080318144457-1jjoztrvqwf0gruk
http://bugs.debian.org/448890
* debian/control:
  - Removed MIA Matt Hope (dopey) from the Uploaders field.
    Thank you for your work, Matt. I hope you are doing well.
  - Downgraded some recommended package to the 'Suggests' priority,
    according to the following discussion on Upstream's mail list.
    http://bioperl.org/pipermail/bioperl-l/2008-March/027379.html
    (Closes: #448890)
* debian/copyright converted to machine-readable format.

Show diffs side-by-side

added added

removed removed

Lines of Context:
Bio/Tools/Genewise.pm
1
 
# $Id: Genewise.pm,v 1.18 2003/08/21 14:53:13 jason Exp $
 
1
# $Id: Genewise.pm,v 1.25.4.1 2006/10/02 23:10:32 sendu Exp $
2
2
#
3
3
# BioPerl module for Bio::Tools::Genewise
4
4
#
29
29
 
30
30
=head1 DESCRIPTION
31
31
 
32
 
This is the parser for the output of Genewise. It takes either a file handle or
33
 
a file name and returns a Bio::SeqFeature::Gene::GeneStructure object.
 
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.
34
35
 
35
36
=head1 FEEDBACK
36
37
 
40
41
Bioperl modules. Send your comments and suggestions preferably to one
41
42
of the Bioperl mailing lists.  Your participation is much appreciated.
42
43
 
43
 
  bioperl-l@bioperl.org          - General discussion
44
 
  http://bio.perl.org/MailList.html             - About the mailing lists
 
44
  bioperl-l@bioperl.org                  - General discussion
 
45
  http://bioperl.org/wiki/Mailing_lists  - About the mailing lists
45
46
 
46
47
=head2 Reporting Bugs
47
48
 
48
49
Report bugs to the Bioperl bug tracking system to help us keep track
49
 
the bugs and their resolution.  Bug reports can be submitted via email
50
 
or the web:
51
 
 
52
 
  bioperl-bugs@bio.perl.org
53
 
  http://bugzilla.bioperl.org/
54
 
 
55
 
=head1 AUTHOR - Fugu Team 
 
50
the bugs and their resolution.  Bug reports can be submitted via the
 
51
web:
 
52
 
 
53
  http://bugzilla.open-bio.org/
 
54
 
 
55
=head1 AUTHOR - Fugu Team, Jason Stajich 
56
56
 
57
57
 Email: fugui@worf.fugu-sg.org
 
58
 Email: jason-at-bioperl.org
58
59
 
59
60
=head1 APPENDIX
60
61
 
61
 
The rest of the documentation details each of the object methods. Internal methods are usually preceded with a _
 
62
The rest of the documentation details each of the object
 
63
methods. Internal methods are usually preceded with a _
62
64
 
63
65
=cut
64
66
 
67
69
 
68
70
 
69
71
package Bio::Tools::Genewise;
70
 
use vars qw(@ISA $Srctag);
 
72
use vars qw($Srctag);
71
73
use strict;
72
74
use Symbol;
73
75
 
74
 
use Bio::Root::Root;
75
 
use Bio::Root::IO;
76
76
use Bio::Tools::AnalysisResult;
77
77
use Bio::SeqFeature::Generic;
78
78
use Bio::SeqFeature::Gene::Exon;
80
80
use Bio::SeqFeature::Gene::Transcript;
81
81
use Bio::SeqFeature::Gene::GeneStructure;
82
82
 
83
 
@ISA = qw(Bio::Root::Root Bio::Root::IO);
 
83
use base qw(Bio::Root::Root Bio::Root::IO);
84
84
$Srctag = 'genewise';
85
85
 
86
86
=head2 new
90
90
           $obj->new(-fh=>\*GW);
91
91
 Function: Constructor for genewise wrapper. Takes either a file or filehandle
92
92
 Example :
93
 
 Returns : L<Bio::Tools::Genewise>
 
93
 Returns : Bio::Tools::Genewise object
 
94
 
 
95
See L<Bio::Tools::Genewise>
94
96
 
95
97
=cut
96
98
 
105
107
 
106
108
 Title   : _get_strand
107
109
 Usage   : $obj->_get_strand
108
 
 Function: takes start and end values, swap them if start>end and returns end
 
110
 Function: takes start and end values, swap them if start>end and 
 
111
           returns end
109
112
 Example :
110
113
 Returns :$start,$end,$strand
111
114
 
113
116
 
114
117
sub _get_strand {
115
118
  my ($self,$start,$end) = @_;
116
 
  $start || $self->throw("Need a start");
117
 
  $end   || $self->throw("Need an end");
 
119
  defined($start) || $self->throw("Need a start");
 
120
  defined($end)   || $self->throw("Need an end");
118
121
  my $strand;
119
122
  if ($start > $end) {
120
123
    my $tmp = $start;
188
191
 Returns : a Bio::SeqFeature::Gene::GeneStructure object
189
192
 Args    :
190
193
 
 
194
See L<Bio::SeqFeature::Gene::GeneStructure>
 
195
 
191
196
=cut
192
197
 
193
198
 
208
213
}
209
214
  
210
215
sub _parse_genes {
211
 
    my ($self) = @_;
212
 
    my @genes;
213
 
    local ($/) = "//";
214
 
    my ($score,$prot_id,$target_id);
215
 
 
216
 
    while ( defined($_ = $self->_readline) ) {
217
 
        $self->debug( $_ ) if( $self->verbose > 0);
218
 
        ($score) = $_=~ m/Score\s+(\d+[\.][\d]+)/o;
219
 
        $self->_score($score) unless defined $self->_score;
220
 
        ($prot_id) = $_=~ m/Query (?:protein|model):\s+(\S+)/o;
221
 
        $self->_prot_id($prot_id) unless defined $self->_prot_id;
222
 
        ($target_id) = $_=~  m/Target Sequence\s+(\S+)/o;       
223
 
        $self->_target_id($target_id) unless defined $self->_target_id;
224
 
        next unless /Gene\s+\d+\n/o;
225
 
 
226
 
        my @genes_txt = split(/Gene\s+\d+\n/);
227
 
        shift @genes_txt; #remove first empty entry
 
216
        my ($self) = @_;
 
217
        my @genes;
 
218
        local ($/) = "//";
 
219
 
 
220
        while ( defined($_ = $self->_readline) ) {
 
221
                $self->debug( $_ ) if( $self->verbose > 0);
 
222
                if( /Score\s+(\-?\d+(\.\d+)?)/ ) {
 
223
              $self->_score($1);# unless defined $self->_score;    
 
224
      } 
 
225
      if( /Query\s+(?:protein|model)\:\s+(\S+)/ ) {
 
226
              $self->_prot_id($1); #unless defined $self->_prot_id;
 
227
           } 
 
228
        
 
229
     if( /Target Sequence\s+(\S+)/ ) {  
 
230
            $self->_target_id($1);# unless defined $self->_target_id;
 
231
          }
 
232
     next unless /Gene\s+\d+\n/;
 
233
 
 
234
     my @genes_txt = split(/Gene\s+\d+\n/,$_);
 
235
     shift @genes_txt; #remove first empty entry
228
236
       
229
 
        foreach my $gene_txt (@genes_txt) {
 
237
     foreach my $gene_txt (@genes_txt) {
230
238
            # If genewise has assigned a strand to the gene as a whole
231
239
            # overall gene start and end
232
 
            my ($g_start, $g_end, 
233
 
                $type) = $gene_txt =~ m/Gene\s+
234
 
                                       (\d+)\s+       # start (1-based)
235
 
                                       (\d+)\s+       # end
236
 
                                       (?:\[(\w+)\])? # 
237
 
                                       /ox;
 
240
            my ($g_start, $g_end, $type) = 
 
241
                        $gene_txt =~ m/Gene\s+
 
242
                                                                (\d+)[\s-]+    # start (1-based)
 
243
                                                                (\d+)\s+       # end
 
244
                                                                (?:\[(\w+)\])? # 
 
245
                                                                /x;
238
246
            my $g_strand;
239
247
            my $source_tag = $type ? "$Srctag". "_$type" : $Srctag;
240
248
            my $genes = new Bio::SeqFeature::Gene::GeneStructure
241
 
                (-source => $source_tag);
 
249
                 (-source => $source_tag);
242
250
            my $transcript = new Bio::SeqFeature::Gene::Transcript
243
 
                (-source => $source_tag,
 
251
                 (-source => $source_tag,
244
252
                 -score  => $self->_score);
245
253
            ($g_start, $g_end, $g_strand) = $self->_get_strand($g_start,$g_end);
246
254
            $genes->strand($g_strand);
247
255
 
248
 
            #grab exon + supporting feature info
 
256
            # grab exon + supporting feature info
249
257
            my @exons;
250
258
            unless ( @exons = $gene_txt =~ m/(Exon .+\s+Supporting .+)/g ) {
251
 
                @exons = $gene_txt =~ m/(Exon .+\s+)/g;
 
259
                    @exons = $gene_txt =~ m/(Exon .+\s+)/g;
252
260
            }
253
261
            my $nbr = 1;
254
 
            #loop through each exon-supporting feature pair
 
262
            # loop through each exon-supporting feature pair
255
263
            foreach my $e (@exons){
256
 
                my ($e_start,$e_end,
257
 
                    $phase) = $e =~ m/Exon\s+
258
 
                                     (\d+)\s+        # start (1 based)
259
 
                                     (\d+)\s+        # end
260
 
                                     phase\s+(\d+)   # phase
261
 
                                     /ox;
 
264
                   my ($e_start,$e_end,$phase) = 
 
265
                 $e =~ m/Exon\s+
 
266
                                        (\d+)[\s-]+     # start (1 based)
 
267
                                             (\d+)\s+        # end
 
268
                                             phase\s+(\d+)   # phase
 
269
                                             /x;
262
270
                my $e_strand;
263
271
                ($e_start,$e_end,$e_strand) = $self->_get_strand($e_start,
264
272
                                                                 $e_end);
278
286
                    $exon->add_tag_value('Target',"Protein:".$self->_prot_id);
279
287
                }
280
288
                $exon->add_tag_value("Exon",$nbr++);
281
 
                if( $e =~ m/Supporting\s+(\d+)\s+(\d+)\s+(\d+)\s+(\d+)/o) {
 
289
                if( $e =~ m/Supporting\s+(\d+)\s+(\d+)\s+(\d+)\s+(\d+)/) {
282
290
                    my ($geno_start,$geno_end,
283
291
                        $prot_start, $prot_end) = ($1,$2,$3,$4);
284
292
                    my $prot_strand;
315
323
                $transcript->add_exon($exon);
316
324
            }
317
325
            $transcript->seq_id($self->_target_id);
 
326
            $transcript->add_tag_value('Id', $self->_prot_id);
318
327
            $genes->add_transcript($transcript);
319
328
            $genes->seq_id($self->_target_id);
320
329
            push @genes, $genes;
322
331
    }
323
332
    $self->{'_genes'} = \@genes;
324
333
}
 
334
 
325
335
1;