← Back to branch summary

~ubuntu-branches/ubuntu/lucid/bioperl/lucid

« back to all changes in this revision

Viewing changes to Bio/Tools/Genemark.pm

  • Committer: Bazaar Package Importer
  • Author(s): Matt Hope
  • Date: 2004-04-18 14:24:11 UTC
  • mfrom: (1.2.1 upstream) (2.1.1 warty)
  • Revision ID: james.westby@ubuntu.com-20040418142411-gr92uexquw4w8liq
Tags: 1.4-1
* New upstream release
* Examples and working code are installed by default to usr/bin,
  this has been moved to usr/share/doc/bioperl/bin

Show diffs side-by-side

added added

removed removed

Lines of Context:
Bio/Tools/Genemark.pm
1
 
# $Id: Genemark.pm,v 1.9 2001/11/20 02:09:40 lstein Exp $
 
1
# $Id: Genemark.pm,v 1.12 2003/04/24 08:49:38 heikki Exp $
2
2
#
3
3
# BioPerl module for Bio::Tools::Genemark
4
4
#
27
27
       # $gene is an instance of Bio::Tools::Prediction::Gene, which inherits
28
28
       # off Bio::SeqFeature::Gene::Transcript.
29
29
       #
30
 
       # $gene->exons() returns an array of 
 
30
       # $gene->exons() returns an array of
31
31
       # Bio::Tools::Prediction::Exon objects
32
32
       # all exons:
33
33
       @exon_arr = $gene->exons();
38
38
       @intrl_exons = $gene->exons('Internal');
39
39
       # terminal exons only
40
40
       @term_exons = $gene->exons('Terminal');
41
 
       # singleton exons: 
 
41
       # singleton exons:
42
42
       ($single_exon) = $gene->exons();
43
43
   }
44
44
 
48
48
 
49
49
=head1 DESCRIPTION
50
50
 
51
 
The Genemark module provides a parser for Genemark gene structure prediction
52
 
output. It parses one gene prediction into a Bio::SeqFeature::Gene::Transcript-
53
 
derived object.
54
 
 
55
 
This module has been developed around genemark.hmm for eukaryots v2.2a and will 
56
 
probably not work with other versions.
57
 
 
58
 
 
59
 
This module also implements the Bio::SeqAnalysisParserI interface, and thus
60
 
can be used wherever such an object fits. See L<Bio::SeqAnalysisParserI>.
 
51
The Genemark module provides a parser for Genemark gene structure
 
52
prediction output. It parses one gene prediction into a
 
53
Bio::SeqFeature::Gene::Transcript- derived object.
 
54
 
 
55
This module has been developed around genemark.hmm for eukaryots v2.2a
 
56
and will probably not work with other versions.
 
57
 
 
58
 
 
59
This module also implements the Bio::SeqAnalysisParserI interface, and
 
60
thus can be used wherever such an object fits. See
 
61
L<Bio::SeqAnalysisParserI>.
61
62
 
62
63
=head1 FEEDBACK
63
64
 
67
68
Bioperl modules. Send your comments and suggestions preferably to one
68
69
of the Bioperl mailing lists.  Your participation is much appreciated.
69
70
 
70
 
  bioperl-l@bioperl.org          - General discussion
71
 
  http://bio.perl.org/MailList.html             - About the mailing lists
 
71
  bioperl-l@bioperl.org                   - General discussion
 
72
  http://bio.perl.org/MailList.html       - About the mailing lists
72
73
 
73
74
=head2 Reporting Bugs
74
75
 
77
78
or the web:
78
79
 
79
80
  bioperl-bugs@bio.perl.org
80
 
  http://bio.perl.org/bioperl-bugs/
 
81
  http://bugzilla.bioperl.org/
81
82
 
82
83
=head1 AUTHOR - Hilmar Lapp, Mark Fiers
83
84
 
84
85
Email hlapp@gmx.net
85
86
      m.w.e.j.fiers@plant.wag-ur.nl
86
87
 
87
 
Describe contact details here
88
 
 
89
88
=head1 APPENDIX
90
89
 
91
 
The rest of the documentation details each of the object methods. Internal methods are usually preceded with a _
 
90
The rest of the documentation details each of the object
 
91
methods. Internal methods are usually preceded with a _
92
92
 
93
93
=cut
94
94
 
135
135
=cut
136
136
 
137
137
#-------------
138
 
sub analysis_method { 
 
138
sub analysis_method {
139
139
#-------------
140
 
    my ($self, $method) = @_;  
 
140
    my ($self, $method) = @_;
141
141
    if($method && ($method !~ /Genemark\.hmm/i)) {
142
142
        $self->throw("method $method not supported in " . ref($self));
143
143
    }
155
155
 
156
156
           The returned object is actually a SeqFeatureI implementing object.
157
157
           This method is required for classes implementing the
158
 
           SeqAnalysisParserI interface, and is merely an alias for 
 
158
           SeqAnalysisParserI interface, and is merely an alias for
159
159
           next_prediction() at present.
160
160
 
161
161
 Example :
196
196
 
197
197
    # get next gene structure
198
198
    $gene = $self->_prediction();
199
 
    
 
199
 
200
200
    return $gene;
201
201
}
202
202
 
207
207
 Function: Parses the prediction section. Automatically called by
208
208
           next_prediction() if not yet done.
209
209
 Example :
210
 
 Returns : 
 
210
 Returns :
211
211
 
212
212
=cut
213
213
 
225
225
    my $current_gene_no = -1;
226
226
 
227
227
    while(defined($_ = $self->_readline())) {
228
 
        
 
228
 
229
229
        if( (/^\s*(\d+)\s+(\d+)/) || (/^\s*(\d+)\s+[\+\-]/)) {
230
230
 
231
231
            #  this is an exon, Genemark doesn't predict anything else
235
235
            #exon no:
236
236
            my $signalnr = 0;
237
237
            if ($2) { my $signalnr = $2; } # used in tag: exon_no
238
 
            
 
238
        
239
239
            # split into fields
240
240
            chomp();
241
241
            my @flds = split(' ', $_);
243
243
            # create the feature (an exon) object
244
244
            my $predobj = Bio::Tools::Prediction::Exon->new();
245
245
 
246
 
                 
 
246
                
247
247
            # define info depending on it being eu- or prokaryot
248
248
            my ($start, $end, $orientation, $prediction_source);
249
249
 
253
253
                ($start, $end) = @flds[(2,3)];
254
254
                $exontag = "_na_";
255
255
 
256
 
            } else {               
 
256
            } else {            
257
257
                $prediction_source = "Genemark.hmm.eu";
258
258
                $orientation = ($flds[2] eq '+') ? 1 : -1;
259
259
                ($start, $end) = @flds[(4,5)];
275
275
                
276
276
            # frame calculation as in the genscan module
277
277
            # is to be implemented...
278
 
            
 
278
        
279
279
            #If the $prednr is not equal to the current gene, we
280
280
            #need to make a new gene and close the old one
281
281
            if($prednr != $current_gene_no) {
282
282
                # a new gene, store the old one if it exists
283
283
                if (defined ($gene)) {
284
 
                    $gene->seqname($seqname);
285
 
                    $self->_add_prediction($gene);          
 
284
                    $gene->seq_id($seqname);
286
285
                    $gene = undef ;
287
286
                }
288
287
                #and make a new one
290
289
                    (
291
290
                     '-primary' => "GenePrediction$prednr",
292
291
                     '-source' => $prediction_source);
293
 
                
 
292
                $self->_add_prediction($gene);          
294
293
                $current_gene_no = $prednr;
295
 
            } 
296
 
            
 
294
            }
 
295
        
297
296
            # Add the exon to the gene
298
297
            $gene->add_exon($predobj, ($exontag eq "_na_" ?
299
298
                                       undef : $exontags{$exontag}));
303
302
        if(/^(Genemark\.hmm\s*[PROKARYOTIC]*)\s+\(Version (.*)\)$/i) {
304
303
            $self->analysis_method($1);
305
304
 
306
 
            my $gm_version = $2; 
 
305
            my $gm_version = $2;
307
306
 
308
307
            $self->analysis_method_version($gm_version);
309
308
            next;
313
312
       if (/^Matrices file:\s+(\S+)?/i)  {
314
313
            $self->analysis_subject($1);
315
314
            # since the line after the matrix file is always the date
316
 
            # (in the output file's I have seen!) extract and store this 
 
315
            # (in the output file's I have seen!) extract and store this
317
316
            # here
318
317
             if (defined(my $_date = $self->_readline())) {
319
318
                 chomp ($_date);
320
319
                 $self->analysis_date($_date);
321
320
             }
322
 
        }                          
 
321
        }                       
323
322
        
324
323
        #Matrix file for prokaryot version
325
324
       if (/^Model file name:\s+(\S+)/) {
326
325
            $self->analysis_subject($1);
327
326
            # since the line after the matrix file is always the date
328
 
            # (in the output file's I have seen!) extract and store this 
 
327
            # (in the output file's I have seen!) extract and store this
329
328
            # here
330
329
            my $_date = $self->_readline() ;
331
330
            if (defined($_date = $self->_readline())) {
350
349
            while (1) {
351
350
               my ($aa_id, $seq) = $self->_read_fasta_seq();
352
351
               last unless ($aa_id);
353
 
                        
 
352
 
354
353
               #now parse through the predictions to add the pred. protein
355
354
               FINDPRED: foreach my $gene (@{$self->{'_preds'}}) {
356
355
                    $gene->primary_tag() =~ /[^0-9]([0-9]+)$/;
362
361
                                                      '-alphabet' => "protein");
363
362
                        $gene->predicted_protein($seqobj);
364
363
                        last FINDPRED;
365
 
                    }     
 
364
                    }   
366
365
 
367
366
               }
368
 
           }                                  
 
367
           }                            
369
368
 
370
369
           last;
371
370
        };
372
371
    }
373
 
    
 
372
 
374
373
    # if the analysis query object contains a ref to a Seq of PrimarySeq
375
374
    # object, then extract the predicted sequences and add it to the gene
376
375
    # object.
377
 
    if (defined $self->analysis_query()) { 
 
376
    if (defined $self->analysis_query()) {
378
377
        my $orig_seq = $self->analysis_query();
379
 
        FINDPREDSEQ: foreach my $gene (@{$self->{'_preds'}}) { 
 
378
        FINDPREDSEQ: foreach my $gene (@{$self->{'_preds'}}) {
380
379
           my $predseq = "";
381
380
           foreach my $exon ($gene->exons()) {
382
381
                #print $exon->start() . " " . $exon->end () . "\n";
388
387
           $gene->predicted_cds($seqobj);
389
388
        }
390
389
    }
391
 
    
392
 
    
 
390
 
 
391
 
393
392
    $self->_predictions_parsed(1);
394
393
}
395
394
 
399
398
 Usage   : $gene = $obj->_prediction()
400
399
 Function: internal
401
400
 Example :
402
 
 Returns : 
 
401
 Returns :
403
402
 
404
403
=cut
405
404
 
416
415
 Usage   : $obj->_add_prediction($gene)
417
416
 Function: internal
418
417
 Example :
419
 
 Returns : 
 
418
 Returns :
420
419
 
421
420
=cut
422
421
 
485
484
    my ($self) = @_;
486
485
    my ($id, $seq);
487
486
    local $/ = ">";
488
 
    
 
487
 
489
488
    return 0 unless (my $entry = $self->_readline());
490
 
    
 
489
 
491
490
    $entry =~ s/^>//;
492
491
    # complete the entry if the first line came from a pushback buffer
493
492
    while(! ($entry =~ />$/)) {
498
497
    # delete everything onwards from an new fasta start (>)
499
498
    $entry =~ s/\n>.*$//s;
500
499
    # id and sequence
501
 
    
 
500
 
502
501
    if($entry =~ s/^(.+)\n//) {
503
502
        $id = $1;
504
503
        $id =~ s/ /_/g;
509
508
        $self->throw("Can't parse Genemark predicted sequence entry");
510
509
    }
511
510
    $seq =~ s/\s//g; # Remove whitespace
512
 
    return ($id, $seq);  
 
511
    return ($id, $seq);
513
512
}
514
513
 
515
514
1;