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

« back to all changes in this revision

Viewing changes to Bio/SeqIO/game/featHandler.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
* 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:
1
 
# $Id: 
 
1
# $Id: featHandler.pm,v 1.15.4.1 2006/10/02 23:10:30 sendu Exp $
2
2
3
3
#
4
4
# Helper module for Bio::SeqIO::game::featHandler
5
5
#
6
 
# Cared for by Sheldon McKay <smckay@bcgsc.bc.ca>
7
 
#
8
 
# Copyright Sheldon McKay
 
6
# Cared for by Sheldon McKay <mckays@cshl.edu>
9
7
#
10
8
# You may distribute this module under the same terms as perl itself
11
9
#
22
20
 
23
21
=head1 DESCRIPTION
24
22
 
25
 
Bio::SeqIO::game::featHandler converts game XML <annotation> elements into
26
 
flattened Bio::SeqFeature::Generic objects to be added to the sequence
 
23
Bio::SeqIO::game::featHandler converts game XML E<lt>annotationE<gt>
 
24
elements into flattened Bio::SeqFeature::Generic objects to be added
 
25
to the sequence
27
26
 
28
27
=head1 FEEDBACK
29
28
 
36
35
Your participation is much appreciated.
37
36
 
38
37
  bioperl-l@bioperl.org                  - General discussion
39
 
  http://bioperl.org/MailList.shtml      - About the mailing lists
 
38
  http://bioperl.org/wiki/Mailing_lists  - About the mailing lists
40
39
 
41
40
=head2 Reporting Bugs
42
41
 
43
42
Report bugs to the Bioperl bug tracking system to help us keep track
44
 
of the bugs and their resolution.
45
 
 
46
 
_Bug reports can be submitted via email or the web:
47
 
 
48
 
  bioperl-bugs@bioperl.org
49
 
  http://bugzilla.bioperl.org/
 
43
of the bugs and their resolution. Bug reports can be submitted via the
 
44
web:
 
45
 
 
46
  http://bugzilla.open-bio.org/
50
47
 
51
48
=head1 AUTHOR - Sheldon McKay
52
49
 
53
 
Email smckay@bcgsc.bc.ca
 
50
Email mckays@cshl.edu
54
51
 
55
52
=head1 APPENDIX
56
53
 
61
58
 
62
59
package Bio::SeqIO::game::featHandler;
63
60
 
64
 
use Bio::SeqIO::game::gameSubs;
65
61
use Bio::SeqFeature::Generic;
66
62
use Bio::Location::Split;
 
63
use Data::Dumper;
67
64
use strict;
68
65
 
69
 
use vars qw { @ISA };                                                                                
 
66
use vars qw {};                                                                                
70
67
 
71
 
@ISA = qw { Bio::SeqIO::game::gameSubs };
 
68
use base qw(Bio::SeqIO::game::gameSubs);
72
69
 
73
70
=head2 new
74
71
 
139
136
sub has_gene {
140
137
    my ($self, $gene, $gname, $id) = @_;
141
138
    
 
139
    # use name preferentially over id. We can't edit IDs in Apollo
 
140
    # AFAIK, and this will create an orphan CDS for newly created 
 
141
    # transcipts -- I think this needs more work
 
142
    #$id = $gname if $id && $gname;
 
143
 
142
144
    unless ( $gene ) {
143
145
        if ( defined $self->{curr_gene} ) {
144
146
            return $self->{curr_gene};
151
153
        if ( $id && !$self->{curr_ltag} ) {
152
154
            $self->{curr_ltag} = $id;
153
155
        }
154
 
        if ( $gname && $gname ne $id && !$self->{curr_gname} ) {
 
156
        if ( $gname && !$self->{curr_gname} ) {
155
157
            $self->{curr_gname} = $gname;
156
158
        }
157
159
            
173
175
        my $feat = Bio::SeqFeature::Generic->new( 
174
176
            -primary => 'gene',
175
177
        );
176
 
 
 
178
        my %seen;
177
179
        for ( keys %{$tags} ) {
178
180
            for my $val ( @{$tags->{$_}} ) {
179
 
                $feat->add_tag_value( $_ => $val );
 
181
                $feat->add_tag_value( $_ => $val ) unless ++$seen{$_.$val} > 1;
180
182
            }
181
183
        }
182
184
            
208
210
    }
209
211
    else {
210
212
        my $tags = $self->{curr_tags};
211
 
 
212
 
        unless ( defined $tags->{product} ){
213
 
            return 0;
214
 
        }
215
 
        
216
213
        $self->{curr_cds} = $self->_add_CDS( $transcript, $tags );
217
214
    }
218
215
}
245
242
 
246
243
    if ( $type eq 'gene' ) {
247
244
        $feat = $self->has_gene;
248
 
        $feat->add_tag_value( standard_name => $id );
249
 
        $feat->add_tag_value( gene => ($self->{curr_gname} || $id) );
250
 
        
 
245
        $feat->add_tag_value( gene => ($self->{curr_gname} || $id) )
 
246
            unless $feat->has_tag('gene');
251
247
    }
252
248
    else {
253
249
        $feat = Bio::SeqFeature::Generic->new;
254
250
        $feat->primary_tag($type);
255
251
        my $gene = $self->has_gene;
256
 
        $gene->add_tag_value( standard_name => $id );
257
 
        $gene->add_tag_value( gene => ($self->{curr_gname} || $id) );
 
252
        $gene->add_tag_value( gene => ($self->{curr_gname} || $id) )
 
253
            unless $gene->has_tag('gene');
 
254
        $feat->add_tag_value( gene => ($self->{curr_gname} || $id) )
 
255
            unless $feat->has_tag('gene');;
258
256
    }
259
257
    for ( keys %{$tags} ) {
260
258
        # or else add simple tag/value pairs
261
259
        if ( $_ eq 'name' && $tags->{type}->[0] eq 'gene' ) {
262
 
            unless ( $feat->has_tag( 'gene' ) ) {
263
 
                $feat->add_tag_value( gene => $tags->{name}->[0] );
264
 
            }
 
260
            $feat->add_tag_value( gene => $tags->{name}->[0] )
 
261
                unless $feat->has_tag( 'gene' );
265
262
            delete $tags->{name};
266
263
        }
267
264
        else {
288
285
        $gene->start( $self->{curr_coords}->[0] );
289
286
        $gene->end( $self->{curr_coords}->[-1] );
290
287
        push @annotations, $gene;
291
 
        $self->{curr_gene} = {};
 
288
        $self->{curr_gene} = '';
292
289
    }
293
290
 
294
291
    # add the subfeatures
310
307
    }
311
308
 
312
309
    # garbage collection
 
310
    $self->{curr_gene}   = '';
313
311
    $self->{curr_ltag}   = '';
314
312
    $self->{curr_gname}  = '';
315
313
    $self->{curr_coords} = [];
336
334
 
337
335
sub _add_generic_annotation {
338
336
    my ($self, $seq, $type, $id, $tags, $feats) = @_;
 
337
    
 
338
    for ( @$feats ) {
 
339
        $_->primary_tag($type);
 
340
    }
339
341
 
340
342
    push @{$self->{ann_l}}, @$feats;
341
343
 
370
372
    $self->{curr_subfeats} = [];
371
373
    $self->{curr_strand}   = 0;
372
374
    my @feats = ();
373
 
    my $generic;
374
375
    my $tags = $self->{curr_tags};
375
376
    my $sname = $set->{_name}->{Characters} ||
376
377
        $set->{Attributes}->{id};
385
386
        $self->_build_feature_set($set, 1);
386
387
        my ($feat) = @{$self->{curr_subfeats}};
387
388
        $feat->primary_tag('transcript') if $feat->primary_tag eq 'exon';
388
 
 
 
389
        if ( $feat->primary_tag eq 'transcript' ) {
 
390
            $feat->add_tag_value( gene => ($gname || $id) )
 
391
                unless $feat->has_tag('gene');
 
392
        }
 
393
        
 
394
        my %seen_tag;
389
395
        for my $tag ( keys %{$tags} ) {
390
396
            for my $val ( @{$tags->{$tag}} ) {
391
 
                $feat->add_tag_value( $tag => $val ) if $val;
 
397
                $feat->add_tag_value( $tag => $val ) 
 
398
                    if $val && ++$seen_tag{$tag.$val} < 2;
392
399
            }
393
400
        }
394
 
        $feat->add_tag_value( gene => ($gname || $id) );
395
401
        @feats = ($feat);
396
402
    }
397
403
    else {
399
405
        $self->{curr_cds}      = '';
400
406
        $gname = $id if $gname eq 'gene';
401
407
        $self->{curr_gname} = $gname;
402
 
    
 
408
 
403
409
        if ( $self->has_gene ) {
404
410
            unless ( $anntype =~/RNA/i ) {
405
411
                $stype =~ s/transcript/mRNA/;
413
419
        my $feat = $self->{curr_feat};
414
420
        $self->_build_feature_set($set);
415
421
        
416
 
        if ( $self->{curr_ltag} || $gname ) {
417
 
            $self->{curr_feat}->add_tag_value( gene => ($gname || $self->{curr_ltag}) );
418
 
        }
 
422
        my $gene = $gname || $self->{curr_ltag};
419
423
        
 
424
        $feat->add_tag_value( gene => $gene )
 
425
            unless $feat->has_tag('gene');
 
426
 
420
427
        # if there is an annotated protein product
421
428
        my $cds = $self->_has_CDS( $feat );
422
429
 
423
 
        my $gene = $gname || $self->{curr_ltag};
424
 
        $feat->add_tag_value( gene => $gene );
425
 
        
426
430
        if ( $cds ) {
427
 
            my $gene = $gname || $self->{curr_ltag};
 
431
            $feat->primary_tag('mRNA');
 
432
 
428
433
            # we really just want one value here
429
434
            $cds->remove_tag('standard_name') if $cds->has_tag('standard_name');
430
435
            $cds->add_tag_value( standard_name => $sname );
431
436
            $cds->remove_tag('gene') if $cds->has_tag('gene');
432
437
            $cds->add_tag_value( gene => $gene );
433
438
            
434
 
            # catch missing/empty protein ids
435
 
            if ( $cds->has_tag('protein_id' )) {
436
 
                if ( !$cds->get_tag_values('protein_id') ) {
437
 
                    $cds->remove_tag('protein_id');
438
 
                    if ( $cds->has_tag('product') ) {
439
 
                        $cds->add_tag_value($cds->get_tag_values('product'));
440
 
                    }
441
 
                }
 
439
            # catch empty protein ids
 
440
            if ( $cds->has_tag('protein_id' ) && !$cds->get_tag_values('protein_id') ) {
 
441
                my $pid = $self->protein_id($cds, $sname);
 
442
                $cds->remove_tag('protein_id');
 
443
                $cds->add_tag_value( protein_id => $pid );
442
444
            }
443
445
 
444
446
            # make sure other subfeats are tied to the transcript
445
447
            # via a 'standard_name' qualifier and the gene via a 'gene' qualifier
446
448
            my @subfeats = @{$self->{curr_subfeats}};
447
449
            for my $sf ( @ subfeats ) {
448
 
                $sf->add_tag_value( standard_name => $sname );
449
 
                $sf->add_tag_value( gene => $gene );
 
450
                $sf->add_tag_value( standard_name => $sname )
 
451
                    unless $sf->has_tag('standard_name');
 
452
                $sf->add_tag_value( gene => $gene )
 
453
                    unless $sf->has_tag('gene');
450
454
            }
451
455
            
452
 
            $feat->add_tag_value( standard_name => $sname );
453
 
 
454
 
            @feats = sort { $a->start <=> $b->start } ($cds, @subfeats);
455
 
            unshift @feats, $feat;
 
456
            $feat->add_tag_value( standard_name => $sname )
 
457
                unless $feat->has_tag('standard_name');
 
458
            $feat->add_tag_value( gene => $gene )
 
459
                unless $feat->has_tag('gene');
 
460
 
 
461
            # if the mRNA and CDS are the same length, the mRNA is redundant
 
462
            # lose the mRNA, steal its tags and give them to the CDS
 
463
            my %seen;
 
464
            if ( $feat->length == $cds->length ) {
 
465
                for my $t ( $feat->all_tags ) {
 
466
                    next if $t =~ /gene|standard_name/;
 
467
                    $cds->add_tag_value( $t => $feat->get_tag_values($t) );
 
468
                }
 
469
                undef $feat;
 
470
            }
 
471
 
 
472
            @feats = sort { $a->start <=> $b->start } ($cds, @subfeats);
 
473
            unshift @feats, $feat if $feat;
456
474
        }
457
475
        else {
458
476
            if ( @{$self->{curr_loc}} > 1 ) {
470
488
            else {
471
489
                $feat->location( $self->{curr_loc}->[0] );
472
490
            }   
473
 
            
 
491
            
 
492
            for ( keys %$tags ) {
 
493
                # expunge duplicate gene attributes
 
494
                next if /gene/ && $feat->has_tag('gene');
 
495
                for my $v ( @{$tags->{$_}} ) {
 
496
                    $feat->add_tag_value( $_ => $v );
 
497
                }
 
498
            }
 
499
 
474
500
            # make sure other subfeats are tied to the transcript
475
501
            my @subfeats = @{$self->{curr_subfeats}};
476
502
            for my $sf ( @ subfeats ) {
477
 
                $sf->add_tag_value( standard_name => $sname );
478
 
                $sf->add_tag_value( gene => $gene );
 
503
                $sf->add_tag_value( standard_name => $sname )
 
504
                    unless $sf->has_tag('standard_name');
 
505
                $sf->add_tag_value( gene => $gene )
 
506
                    unless $sf->has_tag('gene');
479
507
            }
480
508
 
481
509
            @feats = ( $feat, @subfeats );
487
515
    $self->{curr_coords}->[0] ||= 1000000000000;
488
516
    $self->{curr_coords}->[1] ||= -1000000000000;
489
517
    for ( @feats ) {
490
 
        if ( $self->{curr_coords}->[0] > $_->start ) {
 
518
        if ( $self->{curr_coords}->[0] > $_->start ) {
491
519
            $self->{curr_coords}->[0] = $_->start;
492
520
        }
493
521
        if ( $self->{curr_coords}->[1] < $_->end ) {
537
565
            $self->property( $child );
538
566
        }
539
567
 
540
 
        # need to add the db_xref tags to the gene??
 
568
        # need to add the db_xref tags to the gene?
541
569
        # otherwise, simple tag/value pairs
542
570
        elsif ( $name =~ /synonym|author|description/) {
543
571
            $self->{curr_tags}->{$name} = [$child->{Characters}];
551
579
}
552
580
 
553
581
=head2 _add_feature_span
554
 
 
 
582
 
555
583
 Title   : _add_feature_span
556
584
 Usage   : $self->_add_feature_span($el, 1)
557
585
 Function: an internal method to process <feature_span> elements
609
637
                                       -primary => $type );
610
638
    }
611
639
 
612
 
    # identify the translation product                                             
613
 
    if ( $el->{Attributes}->{produces_seq} ) {
614
 
        my $subseq = $self->{seq_h}->{$el->{Attributes}->{produces_seq}};                       
 
640
    # identify the translation product     
 
641
    my $tscript = $el->{Attributes}->{produces_seq};
 
642
    if ( $tscript && $tscript ne 'null') {
 
643
        my $subseq = $self->{seq_h}->{$el->{Attributes}->{produces_seq}};
615
644
        $self->{curr_tags}->{product} = [$el->{Attributes}->{produces_seq}];
616
 
        $self->{curr_tags}->{translation} = [$subseq->seq];
 
645
        $self->{curr_tags}->{translation} = [$subseq->seq] if $subseq;
617
646
    }      
618
647
 
619
648
    $self->flush( $el );
627
656
 Returns : a Bio::SeqFeature::Generic object
628
657
 Args    : $transcript -- a Bio::SeqFeature::Generic object for a transcript
629
658
           $tags       -- ref. to a hash of tag/value attributes
630
 
           
 
659
 
631
660
=cut
632
661
 
633
662
sub _add_CDS {
639
668
        $loc = Bio::Location::Split->new;
640
669
 
641
670
        # sort the exons in ascending start order
642
 
        my @loc = map  { $_->[1] }
643
 
                  sort { $a->[0] <=> $b->[0] }
644
 
                  map  { [$_->start, $_] } @{$self->{curr_loc}};
 
671
        my @loc = sort { $a->start <=> $b->start } @{$self->{curr_loc}};
645
672
 
646
673
        # then add them to the location object
647
674
        for ( @loc ) {
657
684
    my @exons = $single ? $loc : $loc->sub_Location(1);
658
685
 
659
686
    $feat->location($loc);
660
 
 
661
687
    # try to find a peptide
662
 
    my $seq  = $self->{seq_h}->{ $tags->{product}->[0] } ||
663
 
               $self->{seq_h}->{ $tags->{protein_id}->[0] } ||
664
 
               $self->{seq_h}->{ $tags->{gene}->[0] } ||
665
 
               $self->{seq_h}->{ $tags->{standard_name}->[0] };
666
 
    
667
 
    unless ( $seq ) {
668
 
        $self->warn("I did not find a protein sequence for " . 
669
 
                     $feat->display_name . ".  I can't figure out the CDS boundaries");
670
 
        $feat->location($loc);
671
 
        next;
672
 
    }   
673
 
    my $peptide = $seq->display_id;
 
688
    my $seq = $self->{seq_h}->{ $tags->{protein_id}->[0] };
 
689
    $seq  ||= $self->{seq_h}->{ $tags->{product}->[0] } ||
 
690
              $self->{seq_h}->{ $tags->{gene}->[0] } ||
 
691
              $self->{seq_h}->{ $tags->{standard_name}->[0] };
 
692
     
674
693
 
675
694
    # Can we count on the description format being consistent?
676
 
    # Why is critical CDS info saved as unconstrained text not 
677
 
    # specified in the DTD? 
678
 
    # Anyone have a better idea?    
679
 
    my $desc = $seq->description;
680
 
    $desc =~ s/,//g;
681
 
    $desc =~ s/\]\)/\]\) /g;
682
 
    my ($start, $stop) = ();
683
 
    
684
 
    if ( $desc =~ /cds_boundaries:.+?(\d+)\.\.(\d+)/ ) {
685
 
        ($start, $stop) = ($1 - $self->{offset}, $2 - $self->{offset});
 
695
    # Why is CDS coordinate info saved as description text not 
 
696
    # specified in the DTD?  Anyone have a better idea? Aww,
 
697
    # who am I kidding, I'm the only one who will ever read this!
 
698
    my ($start, $stop, $peptide) = ();
 
699
    if ( $seq ) {
 
700
        $peptide = $seq->display_id;
 
701
        my $desc = $seq->description || '';
 
702
        $desc =~ s/,|\n//g;
 
703
        $desc =~ s/\)(\w)/\) $1/g;
 
704
 
 
705
        if ( $desc =~ /cds_boundaries:.+?(\d+)\.\.(\d+)/ ) {
 
706
            ($start, $stop) = ($1 - $self->{offset}, $2 - $self->{offset});
 
707
        }
 
708
        else {
 
709
            # OK, I guess the transcript must be the CDS then
 
710
            $start = $loc->start;
 
711
            $stop  = $loc->end;
 
712
        }
686
713
    }
687
714
    else {
688
 
        $self->warn("I could not find the CDS boundaries");
689
 
        $feat->location($loc);
690
 
        next;
 
715
        $self->warn("I did not find a protein sequence for " . $feat->display_name);
691
716
    }
692
 
    
 
717
 
693
718
    delete $tags->{transcript};
694
719
    
695
720
    # now chop off the UTRs to create a CDS
696
721
    my @exons_to_add = ();
 
722
    #warn scalar(@exons), " exons, $start, $stop\n";
697
723
    for ( @exons ) {
698
724
        my $exon = Bio::Location::Simple->new;
699
725
 
700
726
        if ( $_->end < $start || $_->start > $stop ) {
 
727
            #warn "exon out of range\n";
701
728
            next;
702
729
        }
703
730
        if ( $_->start < $start && $_->end > $start ) {
 
731
            #warn "chopping off left UTR\n";
704
732
            $exon->start( $start );
705
733
        }
706
734
        if ( $_->end > $stop && $_->start < $stop ) {
 
735
            #warn "chopping off right UTR\n";
707
736
            $exon->end( $stop );
708
737
        }
709
738
 
725
754
    }
726
755
 
727
756
    my $parent = $self->{curr_gname} || $self->{curr_ltag};
728
 
    my $cds_tags = $tags; 
729
 
    $tags = {};
730
757
    
731
 
    #delete $feat->{_gsf_tag_hash};
732
 
    for ( $feat->all_tags ) {
733
 
        $feat->remove_tag($_);
 
758
    # try not to steal too many mRNA attributes for the CDS
 
759
    my $cds_tags = {};
 
760
    for my $k ( keys %$tags ) {
 
761
        if ( $k =~ /product|protein|translation|codon_start/ ) {
 
762
            $cds_tags->{$k} = $tags->{$k};
 
763
            delete $tags->{$k};
 
764
        }
 
765
    } 
 
766
 
 
767
    for ( keys %$tags ) {
 
768
        for my $v ( @{$tags->{$_}} ) {
 
769
            $feat->add_tag_value( $_ => $v )
 
770
                unless $feat->has_tag($_);
 
771
        }
734
772
    }    
735
773
 
736
 
    $cds_tags->{product_desc} = [$desc];
737
 
 
738
774
    if ( $self->{curr_gname} ) {
739
775
        $cds_tags->{gene} = [$self->{curr_gname}];      
740
776
    }
749
785
    $cds_tags->{translation} = [$seq->seq];
750
786
  
751
787
    for ( keys %{$cds_tags} ) {
 
788
        my %seen;
752
789
        for my $val (@{$cds_tags->{$_}}) {
 
790
            next if ++$seen{$val} > 1;
753
791
            $cds->add_tag_value( $_ => $val );
754
792
        }        
755
793
    }
756
794
    
757
795
    $cds;
758
 
 
759
796
}
760
797
 
761
798
1;