← Back to branch summary

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

« back to all changes in this revision

Viewing changes to Bio/Tools/Fgenesh.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/Fgenesh.pm
 
1
# $Id: Fgenesh.pm,v 1.10.4.2 2006/10/02 23:10:32 sendu Exp $
 
2
#
 
3
# BioPerl module for Bio::Tools::Fgenesh
 
4
#
 
5
# Cared for by Christopher Dwan (chris@dwan.org)
 
6
#
 
7
# Copied, lock stock & barrel from Genscan.pm
 
8
#
 
9
# You may distribute this module under the same terms as perl itself
 
10
 
 
11
# POD documentation - main docs before the code
 
12
 
 
13
=head1 NAME
 
14
 
 
15
Bio::Tools::Fgenesh - parse results of one Fgenesh run
 
16
 
 
17
=head1 SYNOPSIS
 
18
 
 
19
   use Bio::Tools::Fgenesh;
 
20
 
 
21
   $fgenesh = Bio::Tools::Fgenesh->new(-file => 'result.fgenesh');
 
22
   # filehandle:
 
23
   $fgenesh = Bio::Tools::Fgenesh->new( -fh  => \*INPUT );
 
24
 
 
25
   # parse the results
 
26
   # note: this class is-a Bio::Tools::AnalysisResult which implements
 
27
   # Bio::SeqAnalysisParserI, i.e., $fgensh->next_feature() is the same
 
28
   while($gene = $fgenesh->next_prediction()) {
 
29
       # $gene is an instance of Bio::Tools::Prediction::Gene, which inherits
 
30
       # off Bio::SeqFeature::Gene::Transcript.
 
31
       #
 
32
       # $gene->exons() returns an array of 
 
33
       # Bio::Tools::Prediction::Exon objects
 
34
       # all exons:
 
35
       @exon_arr = $gene->exons();
 
36
 
 
37
       # initial exons only
 
38
       @init_exons = $gene->exons('Initial');
 
39
       # internal exons only
 
40
       @intrl_exons = $gene->exons('Internal');
 
41
       # terminal exons only
 
42
       @term_exons = $gene->exons('Terminal');
 
43
       # singleton exons: 
 
44
       ($single_exon) = $gene->exons();
 
45
   }
 
46
 
 
47
   # essential if you gave a filename at initialization (otherwise the file
 
48
   # will stay open)
 
49
   $fgenesh->close();
 
50
 
 
51
=head1 DESCRIPTION
 
52
 
 
53
The Fgenesh module provides a parser for Fgenesh (version 2) gene structure 
 
54
prediction output. It parses one gene prediction into a 
 
55
Bio::SeqFeature::Gene::Transcript- derived object.
 
56
 
 
57
This module also implements the L<Bio::SeqAnalysisParserI> interface, and thus
 
58
can be used wherever such an object fits. 
 
59
 
 
60
=head1 FEEDBACK
 
61
 
 
62
=head2 Mailing Lists
 
63
 
 
64
User feedback is an integral part of the evolution of this and other
 
65
Bioperl modules. Send your comments and suggestions preferably to one
 
66
of the Bioperl mailing lists.  Your participation is much appreciated.
 
67
 
 
68
  bioperl-l@bioperl.org                  - General discussion
 
69
  http://bioperl.org/wiki/Mailing_lists  - About the mailing lists
 
70
 
 
71
=head2 Reporting Bugs
 
72
 
 
73
Report bugs to the Bioperl bug tracking system to help us keep track
 
74
the bugs and their resolution.  Bug reports can be submitted via the
 
75
web:
 
76
 
 
77
  http://bugzilla.open-bio.org/
 
78
 
 
79
=head1 AUTHOR - Chris Dwan
 
80
 
 
81
Email chris-at-dwan.org
 
82
 
 
83
=head1 APPENDIX
 
84
 
 
85
The rest of the documentation details each of the object methods. Internal methods are usually preceded with a _
 
86
 
 
87
=cut
 
88
 
 
89
 
 
90
# Let the code begin...
 
91
 
 
92
 
 
93
package Bio::Tools::Fgenesh;
 
94
use strict;
 
95
use Symbol;
 
96
 
 
97
use Bio::Root::Root;
 
98
use Bio::Tools::Prediction::Gene;
 
99
use Bio::Tools::Prediction::Exon;
 
100
 
 
101
use base qw(Bio::Tools::AnalysisResult);
 
102
 
 
103
my %ExonTags = ('CDSf' => 'Initial',
 
104
                'CDSi' => 'Internal',
 
105
                'CDSl' => 'Terminal',
 
106
                'CDSo' => 'Singleton');
 
107
    
 
108
sub _initialize_state {
 
109
    my ($self,@args) = @_;
 
110
    
 
111
    # first call the inherited method!
 
112
    $self->SUPER::_initialize_state(@args);
 
113
 
 
114
    # our private state variables
 
115
    $self->{'_preds_parsed'} = 0;
 
116
    $self->{'_has_cds'} = 0;
 
117
    # array of pre-parsed predictions
 
118
    $self->{'_preds'} = [];
 
119
    # seq stack
 
120
    $self->{'_seqstack'} = [];
 
121
}
 
122
 
 
123
=head2 analysis_method
 
124
 
 
125
 Usage     : $genscan->analysis_method();
 
126
 Purpose   : Inherited method. Overridden to ensure that the name matches
 
127
             /genscan/i.
 
128
 Returns   : String
 
129
 Argument  : n/a
 
130
 
 
131
=cut
 
132
 
 
133
#-------------
 
134
sub analysis_method { 
 
135
#-------------
 
136
    my ($self, $method) = @_;  
 
137
    if($method && ($method !~ /fgenesh/i)) {
 
138
        $self->throw("method $method not supported in " . ref($self));
 
139
    }
 
140
    return $self->SUPER::analysis_method($method);
 
141
}
 
142
 
 
143
=head2 next_feature
 
144
 
 
145
 Title   : next_feature
 
146
 Usage   : while($gene = $fgenesh->next_feature()) {
 
147
                  # do something
 
148
           }
 
149
 Function: Returns the next gene structure prediction of the Fgenesh result
 
150
           file. Call this method repeatedly until FALSE is returned.
 
151
 
 
152
           The returned object is actually a SeqFeatureI implementing object.
 
153
           This method is required for classes implementing the
 
154
           SeqAnalysisParserI interface, and is merely an alias for 
 
155
           next_prediction() at present.
 
156
 
 
157
 Example :
 
158
 Returns : A Bio::Tools::Prediction::Gene object.
 
159
 Args    :
 
160
 
 
161
=cut
 
162
 
 
163
sub next_feature {
 
164
    my ($self,@args) = @_;
 
165
    # even though next_prediction doesn't expect any args (and this method
 
166
    # does neither), we pass on args in order to be prepared if this changes
 
167
    # ever
 
168
    return $self->next_prediction(@args);
 
169
}
 
170
 
 
171
=head2 next_prediction
 
172
 
 
173
 Title   : next_prediction
 
174
 Usage   : while($gene = $fgenesh->next_prediction()) { ... }
 
175
 Function: Returns the next gene structure prediction of the Genscan result
 
176
           file. Call this method repeatedly until FALSE is returned.
 
177
 Example :
 
178
 Returns : A Bio::Tools::Prediction::Gene object.
 
179
 Args    :
 
180
 
 
181
=cut
 
182
 
 
183
sub next_prediction {
 
184
    my ($self) = @_;
 
185
    my $gene;
 
186
 
 
187
    # if the prediction section hasn't been parsed yet, we do this now
 
188
    $self->_parse_predictions() unless $self->_predictions_parsed();
 
189
 
 
190
    # get next gene structure
 
191
    $gene = $self->_prediction();
 
192
 
 
193
    if($gene) {
 
194
        # fill in predicted protein, and if available the predicted CDS
 
195
        #
 
196
 
 
197
        # use the seq stack if there's a seq on it
 
198
        my $seqobj = pop(@{$self->{'_seqstack'}});
 
199
        my ($id, $seq);
 
200
        unless ($seqobj) {
 
201
           ($id, $seq) = $self->_read_fasta_seq();
 
202
           my $alphabet;
 
203
           if (($id =~ /mrna/) || ($id =~ /cds/)) { $alphabet = 'dna'; }
 
204
           else { $alphabet = 'protein'; }
 
205
           $seqobj = Bio::PrimarySeq->new('-seq'        => $seq,
 
206
                                          '-display_id' => $id,
 
207
                                          '-alphabet'   => $alphabet); 
 
208
        }
 
209
        if ($seqobj) {
 
210
 
 
211
            # check that prediction number matches the prediction number
 
212
            # indicated in the sequence id (there may be incomplete gene
 
213
            # predictions that contain only signals with no associated protein
 
214
            # prediction.
 
215
 
 
216
            $gene->primary_tag() =~ /[^0-9]([0-9]+)$/;
 
217
            my $prednr = $1;
 
218
            if ($id !~ /_predicted_(\w+)_$prednr/) {
 
219
                # this is not our sequence, so push back for next prediction
 
220
                push(@{$self->{'_seqstack'}}, $seqobj);
 
221
            } else {
 
222
                if ($1 eq "protein") {
 
223
                  $gene->predicted_protein($seqobj);
 
224
                } elsif (($1 eq "mrna") || ($1 eq "cds")) {
 
225
                  $self->_has_cds(1);
 
226
                  $gene->predicted_cds($seqobj);
 
227
                  
 
228
                  # Have to go back in and get the protein...
 
229
                  ($id, $seq) = $self->_read_fasta_seq();
 
230
                  if ($id =~ /_cds_/) { 
 
231
                    ($id, $seq) = $self->_read_fasta_seq(); 
 
232
                  }
 
233
 
 
234
                  $seqobj = Bio::PrimarySeq->new('-seq' => $seq,
 
235
                                                 '-display_id' => $id,
 
236
                                                 '-alphabet' => "protein");
 
237
                  $gene->predicted_protein($seqobj);
 
238
                }
 
239
            }
 
240
        }
 
241
    }
 
242
 
 
243
    return $gene;
 
244
}
 
245
 
 
246
=head2 _parse_predictions
 
247
 
 
248
 Title   : _parse_predictions()
 
249
 Usage   : $obj->_parse_predictions()
 
250
 Function: Parses the prediction section. Automatically called by
 
251
           next_prediction() if not yet done.
 
252
 Example :
 
253
 Returns : 
 
254
 
 
255
=cut
 
256
 
 
257
sub _parse_predictions {
 
258
    my ($self) = @_;
 
259
    my $gene;
 
260
    my $seqname;
 
261
 
 
262
    while(defined($_ = $self->_readline())) {
 
263
 
 
264
        if(/^\s*(\d+)\s+([+\-])/) {
 
265
            my $line = $_;
 
266
 
 
267
            # exon or signal
 
268
            my $prednr = $1;
 
269
            my $strand = ($2 eq '+') ? 1 : -1;
 
270
 
 
271
            if(! defined($gene)) {
 
272
                $gene = Bio::Tools::Prediction::Gene->new(
 
273
                                       '-primary' => "GenePrediction$prednr",
 
274
                                       '-source' => 'Fgenesh');
 
275
            }
 
276
            # split into fields
 
277
            chomp();
 
278
            my @flds = split(/\s+/, $line);
 
279
 
 
280
            # create the feature object depending on the type of signal
 
281
            my $predobj;
 
282
            my $is_exon = grep {$line =~ $_} keys(%ExonTags);
 
283
            my ($start, $end);
 
284
            if($is_exon) {
 
285
                $predobj = Bio::Tools::Prediction::Exon->new();
 
286
                $predobj->score($flds[8]);
 
287
                $start   = $flds[5];
 
288
                $end     = $flds[7];
 
289
            } else {
 
290
                # PolyA site, or TSS 
 
291
                $predobj = Bio::SeqFeature::Generic->new();
 
292
                $predobj->score($flds[5]);
 
293
                $start   = $flds[4];
 
294
                $end     = $flds[4];
 
295
            }
 
296
            # set common fields
 
297
            $predobj->source_tag('Fgenesh');
 
298
            $predobj->strand($strand);
 
299
 
 
300
# Following tactical commenting-out made by
 
301
# malcolm.cook@stowers-institute.org since coordinate reversal is
 
302
# apparently vestigial copy/paste detritus from Genscan.pm origins of
 
303
# this module and this is NOT needed for fgenesh (at least in v
 
304
# 2.1.4).
 
305
 
 
306
#           if($predobj->strand() == 1) {
 
307
                $predobj->start($start);
 
308
                $predobj->end($end);
 
309
#           } else {
 
310
#               $predobj->end($start);
 
311
#               $predobj->start($end);
 
312
#           }
 
313
 
 
314
            # print STDERR "start $start end $end\n";
 
315
            # add to gene structure (should be done only when start and end
 
316
            # are set, in order to allow for proper expansion of the range)
 
317
            if($is_exon) {
 
318
                # first, set fields unique to exons
 
319
                $predobj->primary_tag($ExonTags{$flds[3]} . 'Exon');
 
320
                $predobj->is_coding(1);
 
321
                my $cod_offset;
 
322
                if($predobj->strand() == 1) {
 
323
                    $cod_offset = ($flds[9] - $predobj->start()) % 3;
 
324
                    # Possible values are -2, -1, 0, 1, 2. -1 and -2 correspond
 
325
                    # to offsets 2 and 1, resp. Offset 3 is the same as 0.
 
326
                    $cod_offset += 3 if($cod_offset < 1);                   
 
327
                } else {
 
328
                    # On the reverse strand the Genscan frame also refers to
 
329
                    # the first base of the first complete codon, but viewed
 
330
                    # from forward, which is the third base viewed from
 
331
                    # reverse.
 
332
                    $cod_offset = ($flds[11] - $predobj->end()) % 3;
 
333
                    # Possible values are -2, -1, 0, 1, 2. Due to the reverse
 
334
                    # situation, {2,-1} and {1,-2} correspond to offsets
 
335
                    # 1 and 2, resp. Offset 3 is the same as 0.
 
336
                    $cod_offset -= 3 if($cod_offset >= 0);
 
337
                    $cod_offset = -$cod_offset;
 
338
                }
 
339
                # Offsets 2 and 1 correspond to frame 1 and 2 (frame of exon
 
340
                # is the frame of the first base relative to the exon, or the
 
341
                # number of bases the first codon is missing).
 
342
                $predobj->frame(3 - $cod_offset);
 
343
                # print STDERR "  frame is " . $predobj->frame() . "\n";
 
344
                # then add to gene structure object
 
345
                $gene->add_exon($predobj, $ExonTags{$flds[1]});         
 
346
            } elsif($flds[3] eq 'PolA') {
 
347
                $predobj->primary_tag("PolyAsite");
 
348
                $gene->poly_A_site($predobj);
 
349
            } elsif($flds[3] eq 'TSS') {
 
350
                $predobj->primary_tag("Promoter"); # (hey! a TSS is NOT a promoter... what's going on here?...
 
351
                $gene->add_promoter($predobj);
 
352
                #I'd like to do this (for now):
 
353
                #$predobj->primary_tag("TSS"); #this is not the right model, but, it IS a feature at least.
 
354
                #but the followg errs out
 
355
                #$gene->add_SeqFeature($predobj); #err: MSG: start is undefined when bounds at Bio::SeqFeature::Generic::add_SeqFeature 671 check since gene has no start yet
 
356
            }
 
357
            else {
 
358
              $self->throw("unrecognized prediction line: " . $line);
 
359
            }
 
360
            next;
 
361
        }
 
362
 
 
363
        if(/^\s*$/ && defined($gene)) {
 
364
            # current gene is completed
 
365
            $gene->seq_id($seqname);
 
366
            $self->_add_prediction($gene);
 
367
            $gene = undef;
 
368
            next;
 
369
        }
 
370
 
 
371
        if(/^(FGENESH)\s+([\d\.]+)/) {
 
372
            $self->analysis_method($1);
 
373
            $self->analysis_method_version($2);
 
374
            if (/\s(\S+)\sgenomic DNA/) {
 
375
              $self->analysis_subject($1);
 
376
            }
 
377
            next;
 
378
        }
 
379
 
 
380
        if(/^\s*Seq name:\s+(\S+)/) {
 
381
            $seqname = $1;
 
382
            next;
 
383
        }
 
384
        
 
385
        /^Predicted protein/ && do {
 
386
            # section of predicted sequences
 
387
            $self->_pushback($_);
 
388
            last;
 
389
        };
 
390
    }
 
391
    $self->_predictions_parsed(1);
 
392
}
 
393
 
 
394
=head2 _prediction
 
395
 
 
396
 Title   : _prediction()
 
397
 Usage   : $gene = $obj->_prediction()
 
398
 Function: internal
 
399
 Example :
 
400
 Returns : 
 
401
 
 
402
=cut
 
403
 
 
404
sub _prediction {
 
405
    my ($self) = @_;
 
406
 
 
407
    return unless(exists($self->{'_preds'}) && @{$self->{'_preds'}});
 
408
    return shift(@{$self->{'_preds'}});
 
409
}
 
410
 
 
411
=head2 _add_prediction
 
412
 
 
413
 Title   : _add_prediction()
 
414
 Usage   : $obj->_add_prediction($gene)
 
415
 Function: internal
 
416
 Example :
 
417
 Returns : 
 
418
 
 
419
=cut
 
420
 
 
421
sub _add_prediction {
 
422
    my ($self, $gene) = @_;
 
423
 
 
424
    if(! exists($self->{'_preds'})) {
 
425
        $self->{'_preds'} = [];
 
426
    }
 
427
    push(@{$self->{'_preds'}}, $gene);
 
428
}
 
429
 
 
430
=head2 _predictions_parsed
 
431
 
 
432
 Title   : _predictions_parsed
 
433
 Usage   : $obj->_predictions_parsed
 
434
 Function: internal
 
435
 Example :
 
436
 Returns : TRUE or FALSE
 
437
 
 
438
=cut
 
439
 
 
440
sub _predictions_parsed {
 
441
    my ($self, $val) = @_;
 
442
 
 
443
    $self->{'_preds_parsed'} = $val if $val;
 
444
    if(! exists($self->{'_preds_parsed'})) {
 
445
        $self->{'_preds_parsed'} = 0;
 
446
    }
 
447
    return $self->{'_preds_parsed'};
 
448
}
 
449
 
 
450
=head2 _has_cds
 
451
 
 
452
 Title   : _has_cds()
 
453
 Usage   : $obj->_has_cds()
 
454
 Function: Whether or not the result contains the predicted CDSs, too.
 
455
 Example :
 
456
 Returns : TRUE or FALSE
 
457
 
 
458
=cut
 
459
 
 
460
sub _has_cds {
 
461
    my ($self, $val) = @_;
 
462
 
 
463
    $self->{'_has_cds'} = $val if $val;
 
464
    if(! exists($self->{'_has_cds'})) {
 
465
        $self->{'_has_cds'} = 0;
 
466
    }
 
467
    return $self->{'_has_cds'};
 
468
}
 
469
 
 
470
=head2 _read_fasta_seq
 
471
 
 
472
 Title   : _read_fasta_seq()
 
473
 Usage   : ($id,$seqstr) = $obj->_read_fasta_seq();
 
474
 Function: Simple but specialised FASTA format sequence reader. Uses
 
475
           $self->_readline() to retrieve input, and is able to strip off
 
476
           the traling description lines.
 
477
 Example :
 
478
 Returns : An array of two elements: fasta_id & sequence
 
479
 
 
480
=cut
 
481
 
 
482
sub _read_fasta_seq {
 
483
    my ($self) = @_;
 
484
    my ($id, $seq);
 
485
    #local $/ = ">";
 
486
    
 
487
    my $entry = $self->_readline();
 
488
    # print " ^^ $entry\n";
 
489
    return unless ($entry);
 
490
    $entry = $self->_readline() if ($entry =~ /^Predicted protein/);
 
491
    # print " ^^ $entry\n";
 
492
 
 
493
    # Pick up the header / id.
 
494
    if ($entry =~ /^>FGENESH:/) {
 
495
      if ($entry =~ /^>FGENESH:\s+(\d+)/) {
 
496
         # print STDERR "  this is a predicted gene\n";
 
497
         $id  = "_predicted_protein_" . $1;
 
498
      } elsif ($entry =~ /^>FGENESH:\[mRNA\]\s+(\d+)/) {
 
499
        # print STDERR "  this is an mRNA\n";
 
500
         $id  = "_predicted_mrna_" . $1;
 
501
      } elsif ($entry =~ /^>FGENESH:\[exon\]\s+Gene:\s+(\d+)/) {
 
502
         $id  = "_predicted_cds_"  . $1;
 
503
      }
 
504
      $seq = "";
 
505
      $entry = $self->_readline();
 
506
    }
 
507
 
 
508
    my $done = 0;
 
509
    while (!$done) {
 
510
       # print "*** $entry\n";
 
511
       if (($entry =~ /^>FGENESH:\[exon\]/) && ($id =~ /^_predicted_cds_/)) {
 
512
         # print STDERR "  -- informed about an exon header...\n";
 
513
         $entry = $self->_readline();
 
514
       } else {
 
515
         $seq .= $entry;
 
516
         # print STDERR "  Added $entry\n";
 
517
       }
 
518
 
 
519
       last unless $entry  = $self->_readline();
 
520
       if (($entry =~ /^>/) && 
 
521
           (!(($entry =~ /^>FGENESH:\[exon\]/) && ($id =~ /^_predicted_cds_/)))) {
 
522
          $self->_pushback($entry); last;    
 
523
       }
 
524
    }
 
525
 
 
526
    # id and sequence
 
527
    $seq =~ s/\s//g; # Remove whitespace
 
528
    return ($id, $seq);
 
529
}
 
530
 
 
531
1;