Bio/SeqFeatureI.pm

   1 # $Id$
   2 #
   3 # BioPerl module for Bio::SeqFeatureI
   4 #
   5 # Cared for by Ewan Birney <birney@ebi.ac.uk>
   6 #
   7 # Copyright Ewan Birney
   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::SeqFeatureI - Abstract interface of a Sequence Feature
  16
  17 =head1 SYNOPSIS
  18
  19     # get a seqfeature somehow, eg, from a Sequence with Features attached
  20
  21     foreach $feat ( $seq->get_SeqFeatures() ) {
  22        print "Feature from ", $feat->start, "to ",
  23                $feat->end, " Primary tag  ", $feat->primary_tag,
  24                   ", produced by ", $feat->source_tag(), "\n";
  25
  26        if( $feat->strand == 0 ) {
  27                     print "Feature applicable to either strand\n";
  28        } else {
  29           print "Feature on strand ", $feat->strand,"\n"; # -1,1
  30        }
  31
  32        print "feature location is ",$feat->start, "..",
  33           $feat->end, " on strand ", $feat->strand, "\n";
  34        print "easy utility to print locations in GenBank/EMBL way ",
  35           $feat->location->to_FTstring(), "\n";
  36
  37        foreach $tag ( $feat->get_all_tags() ) {
  38                     print "Feature has tag ", $tag, " with values, ",
  39                       join(' ',$feat->get_tag_values($tag)), "\n";
  40        }
  41             print "new feature\n" if $feat->has_tag('new');
  42             # features can have sub features
  43             my @subfeat = $feat->get_SeqFeatures();
  44          }
  45
  46 =head1 DESCRIPTION
  47
  48 This interface is the functions one can expect for any Sequence
  49 Feature, whatever its implementation or whether it is a more complex
  50 type (eg, a Gene). This object does not actually provide any
  51 implemention, it just provides the definitions of what methods one can
  52 call. See Bio::SeqFeature::Generic for a good standard implementation
  53 of this object
  54
  55 =head1 FEEDBACK
  56
  57 User feedback is an integral part of the evolution of this and other
  58 Bioperl modules. Send your comments and suggestions preferably to one
  59 of the Bioperl mailing lists.  Your participation is much appreciated.
  60
  61   bioperl-l@bioperl.org                  - General discussion
  62   http://bioperl.org/wiki/Mailing_lists  - About the mailing lists
  63
  64 =head2 Reporting Bugs
  65
  66 Report bugs to the Bioperl bug tracking system to help us keep track
  67 the bugs and their resolution.  Bug reports can be submitted via the
  68 web:
  69
  70   http://bugzilla.open-bio.org/
  71
  72 =head1 APPENDIX
  73
  74 The rest of the documentation details each of the object
  75 methods. Internal methods are usually preceded with a _
  76
  77 =cut
  78
  79
  80 # Let the code begin...
  81
  82
  83 package Bio::SeqFeatureI;
  84 use vars qw($HasInMemory);
  85 use strict;
  86 BEGIN {
  87     eval { require Bio::DB::InMemoryCache };
  88     if( $@ ) { $HasInMemory = 0 }
  89     else { $HasInMemory = 1 }
  90 }
  91
  92 use Bio::Seq;
  93
  94 use Carp;
  95
  96 use base qw(Bio::RangeI);
  97
  98 =head1 Bio::SeqFeatureI specific methods
  99
 100 New method interfaces.
 101
 102 =cut
 103
 104 =head2 get_SeqFeatures
 105
 106  Title   : get_SeqFeatures
 107  Usage   : @feats = $feat->get_SeqFeatures();
 108  Function: Returns an array of sub Sequence Features
 109  Returns : An array
 110  Args    : none
 111
 112 =cut
 113
 114 sub get_SeqFeatures{
 115    my ($self,@args) = @_;
 116
 117    $self->throw_not_implemented();
 118 }
 119
 120 =head2 display_name
 121
 122  Title   : display_name
 123  Usage   : $name = $feat->display_name()
 124  Function: Returns the human-readable name of the feature for displays.
 125  Returns : a string
 126  Args    : none
 127
 128 =cut
 129
 130 sub display_name {
 131     shift->throw_not_implemented();
 132 }
 133
 134 =head2 primary_tag
 135
 136  Title   : primary_tag
 137  Usage   : $tag = $feat->primary_tag()
 138  Function: Returns the primary tag for a feature,
 139            eg 'exon'
 140  Returns : a string
 141  Args    : none
 142
 143
 144 =cut
 145
 146 sub primary_tag{
 147    my ($self,@args) = @_;
 148
 149    $self->throw_not_implemented();
 150
 151 }
 152
 153 =head2 source_tag
 154
 155  Title   : source_tag
 156  Usage   : $tag = $feat->source_tag()
 157  Function: Returns the source tag for a feature,
 158            eg, 'genscan'
 159  Returns : a string
 160  Args    : none
 161
 162
 163 =cut
 164
 165 sub source_tag{
 166    my ($self,@args) = @_;
 167
 168    $self->throw_not_implemented();
 169 }
 170
 171 =head2 has_tag
 172
 173  Title   : has_tag
 174  Usage   : $tag_exists = $self->has_tag('some_tag')
 175  Function:
 176  Returns : TRUE if the specified tag exists, and FALSE otherwise
 177  Args    :
 178
 179
 180 =cut
 181
 182 sub has_tag{
 183    my ($self,@args) = @_;
 184
 185    $self->throw_not_implemented();
 186
 187 }
 188
 189 =head2 get_tag_values
 190
 191  Title   : get_tag_values
 192  Usage   : @values = $self->get_tag_values('some_tag')
 193  Function:
 194  Returns : An array comprising the values of the specified tag.
 195  Args    : a string
 196
 197 throws an exception if there is no such tag
 198
 199 =cut
 200
 201 sub get_tag_values {
 202     shift->throw_not_implemented();
 203 }
 204
 205 =head2 get_tagset_values
 206
 207  Title   : get_tagset_values
 208  Usage   : @values = $self->get_tagset_values(qw(label transcript_id product))
 209  Function:
 210  Returns : An array comprising the values of the specified tags, in order of tags
 211  Args    : An array of strings
 212
 213 does NOT throw an exception if none of the tags are not present
 214
 215 this method is useful for getting a human-readable label for a
 216 SeqFeatureI; not all tags can be assumed to be present, so a list of
 217 possible tags in preferential order is provided
 218
 219 =cut
 220
 221 # interface + abstract method
 222 sub get_tagset_values {
 223     my ($self, @args) = @_;
 224     my @vals = ();
 225     foreach my $arg (@args) {
 226         if ($self->has_tag($arg)) {
 227             push(@vals, $self->get_tag_values($arg));
 228         }
 229     }
 230     return @vals;
 231 }
 232
 233 =head2 get_all_tags
 234
 235  Title   : get_all_tags
 236  Usage   : @tags = $feat->get_all_tags()
 237  Function: gives all tags for this feature
 238  Returns : an array of strings
 239  Args    : none
 240
 241
 242 =cut
 243
 244 sub get_all_tags{
 245     shift->throw_not_implemented();
 246 }
 247
 248 =head2 attach_seq
 249
 250  Title   : attach_seq
 251  Usage   : $sf->attach_seq($seq)
 252  Function: Attaches a Bio::Seq object to this feature. This
 253            Bio::Seq object is for the *entire* sequence: ie
 254            from 1 to 10000
 255
 256            Note that it is not guaranteed that if you obtain a feature from
 257            an object in bioperl, it will have a sequence attached. Also,
 258            implementors of this interface can choose to provide an empty
 259            implementation of this method. I.e., there is also no guarantee
 260            that if you do attach a sequence, seq() or entire_seq() will not
 261            return undef.
 262
 263            The reason that this method is here on the interface is to enable
 264            you to call it on every SeqFeatureI compliant object, and
 265            that it will be implemented in a useful way and set to a useful
 266            value for the great majority of use cases. Implementors who choose
 267            to ignore the call are encouraged to specifically state this in
 268            their documentation.
 269
 270  Example :
 271  Returns : TRUE on success
 272  Args    : a Bio::PrimarySeqI compliant object
 273
 274
 275 =cut
 276
 277 sub attach_seq {
 278     shift->throw_not_implemented();
 279 }
 280
 281 =head2 seq
 282
 283  Title   : seq
 284  Usage   : $tseq = $sf->seq()
 285  Function: returns the truncated sequence (if there is a sequence attached)
 286            for this feature
 287  Example :
 288  Returns : sub seq (a Bio::PrimarySeqI compliant object) on attached sequence
 289            bounded by start & end, or undef if there is no sequence attached
 290  Args    : none
 291
 292
 293 =cut
 294
 295 sub seq {
 296     shift->throw_not_implemented();
 297 }
 298
 299 =head2 entire_seq
 300
 301  Title   : entire_seq
 302  Usage   : $whole_seq = $sf->entire_seq()
 303  Function: gives the entire sequence that this seqfeature is attached to
 304  Example :
 305  Returns : a Bio::PrimarySeqI compliant object, or undef if there is no
 306            sequence attached
 307  Args    : none
 308
 309
 310 =cut
 311
 312 sub entire_seq {
 313     shift->throw_not_implemented();
 314 }
 315
 316
 317 =head2 seq_id
 318
 319  Title   : seq_id
 320  Usage   : $obj->seq_id($newval)
 321  Function: There are many cases when you make a feature that you
 322            do know the sequence name, but do not know its actual
 323            sequence. This is an attribute such that you can store
 324            the ID (e.g., display_id) of the sequence.
 325
 326            This attribute should *not* be used in GFF dumping, as
 327            that should come from the collection in which the seq
 328            feature was found.
 329  Returns : value of seq_id
 330  Args    : newvalue (optional)
 331
 332
 333 =cut
 334
 335 sub seq_id {
 336     shift->throw_not_implemented();
 337 }
 338
 339 =head2 gff_string
 340
 341  Title   : gff_string
 342  Usage   : $str = $feat->gff_string;
 343            $str = $feat->gff_string($gff_formatter);
 344  Function: Provides the feature information in GFF format.
 345
 346            The implementation provided here returns GFF2 by default. If you
 347            want a different version, supply an object implementing a method
 348            gff_string() accepting a SeqFeatureI object as argument. E.g., to
 349            obtain GFF1 format, do the following:
 350
 351                 my $gffio = Bio::Tools::GFF->new(-gff_version => 1);
 352                 $gff1str = $feat->gff_string($gff1io);
 353
 354  Returns : A string
 355  Args    : Optionally, an object implementing gff_string().
 356
 357
 358 =cut
 359
 360 sub gff_string{
 361    my ($self,$formatter) = @_;
 362
 363    $formatter = $self->_static_gff_formatter unless $formatter;
 364    return $formatter->gff_string($self);
 365 }
 366
 367 my $static_gff_formatter = undef;
 368
 369 =head2 _static_gff_formatter
 370
 371  Title   : _static_gff_formatter
 372  Usage   :
 373  Function:
 374  Example :
 375  Returns :
 376  Args    :
 377
 378
 379 =cut
 380
 381 sub _static_gff_formatter{
 382    my ($self,@args) = @_;
 383    require Bio::Tools::GFF; # on the fly inclusion -- is this better?
 384    if( !defined $static_gff_formatter ) {
 385        $static_gff_formatter = Bio::Tools::GFF->new('-gff_version' => 2);
 386    }
 387    return $static_gff_formatter;
 388 }
 389
 390
 391 =head1 Decorating methods
 392
 393 These methods have an implementation provided by Bio::SeqFeatureI,
 394 but can be validly overwritten by subclasses
 395
 396 =head2 spliced_seq
 397
 398   Title   : spliced_seq
 399
 400   Usage   : $seq = $feature->spliced_seq()
 401             $seq = $feature_with_remote_locations->spliced_seq($db_for_seqs)
 402
 403   Function: Provides a sequence of the feature which is the most
 404             semantically "relevant" feature for this sequence. A default
 405             implementation is provided which for simple cases returns just
 406             the sequence, but for split cases, loops over the split location
 407             to return the sequence. In the case of split locations with
 408             remote locations, eg
 409
 410             join(AB000123:5567-5589,80..1144)
 411
 412             in the case when a database object is passed in, it will attempt
 413             to retrieve the sequence from the database object, and "Do the right thing",
 414             however if no database object is provided, it will generate the correct
 415             number of N's (DNA) or X's (protein, though this is unlikely).
 416
 417             This function is deliberately "magical" attempting to second guess
 418             what a user wants as "the" sequence for this feature.
 419
 420             Implementing classes are free to override this method with their
 421             own magic if they have a better idea what the user wants.
 422
 423   Args    : [optional]
 424             -db        A L<Bio::DB::RandomAccessI> compliant object if
 425                        one needs to retrieve remote seqs.
 426             -nosort    boolean if the locations should not be sorted
 427                        by start location.  This may occur, for instance,
 428                        in a circular sequence where a gene span starts
 429                        before the end of the sequence and ends after the
 430                        sequence start. Example : join(15685..16260,1..207)
 431                                            (default = if sequence is_circular(), 1, otherwise 0)
 432                         -phase     truncates the returned sequence based on the
 433                                            intron phase (0,1,2).
 434
 435   Returns : A L<Bio::PrimarySeqI> object
 436
 437 =cut
 438
 439 sub spliced_seq {
 440     my $self = shift;
 441         my @args = @_;
 442         my ($db, $nosort, $phase) =
 443            $self->_rearrange([qw(DB NOSORT PHASE)], @args);
 444
 445         # set no_sort based on the parent sequence status
 446         if ($self->entire_seq->is_circular) {
 447                 $nosort = 1;
 448         }
 449
 450         # (added 7/7/06 to allow use old API (with warnings)
 451         my $old_api = (!(grep {$_ =~ /(?:nosort|db|phase)/} @args)) ? 1 : 0;
 452         if (@args && $old_api) {
 453                 $self->warn(q(API has changed; please use '-db' or '-nosort' ).
 454                      qq(for args. See POD for more details.));
 455                 $db = shift @args if @args;
 456                 $nosort = shift @args if @args;
 457                 $phase = shift @args if @args;
 458         };
 459
 460         if (defined($phase) && ($phase < 0 || $phase > 2)) {
 461             $self->warn("Phase must be 0,1, or 2.  Setting phase to 0...");
 462             $phase = 0;
 463         }
 464
 465         if( $db && ref($db) && ! $db->isa('Bio::DB::RandomAccessI') ) {
 466         $self->warn("Must pass in a valid Bio::DB::RandomAccessI object".
 467                     " for access to remote locations for spliced_seq");
 468         $db = undef;
 469     } elsif( defined $db && $HasInMemory &&
 470             $db->isa('Bio::DB::InMemoryCache') ) {
 471         $db = Bio::DB::InMemoryCache->new(-seqdb => $db);
 472     }
 473
 474     if( ! $self->location->isa("Bio::Location::SplitLocationI") ) {
 475         if ($phase) {
 476             $self->debug("Subseq start: ",$phase+1,"\tend: ",$self->end,"\n");
 477             my $seqstr = substr($self->seq->seq, $phase);
 478             my $out = Bio::Seq->new( -id => $self->entire_seq->display_id
 479                                 . "_spliced_feat",
 480                          -seq => $seqstr);
 481             return $out;
 482         } else {
 483             return $self->seq(); # nice and easy!
 484         }
 485     }
 486
 487     # redundant test, but the above ISA is probably not ideal.
 488     if( ! $self->location->isa("Bio::Location::SplitLocationI") ) {
 489         $self->throw("not atomic, not split, yikes, in trouble!");
 490     }
 491
 492     my $seqstr = '';
 493     my $seqid = $self->entire_seq->display_id;
 494     # This is to deal with reverse strand features
 495     # so we are really sorting features 5' -> 3' on their strand
 496     # i.e. rev strand features will be sorted largest to smallest
 497     # as this how revcom CDSes seem to be annotated in genbank.
 498     # Might need to eventually allow this to be programable?
 499     # (can I mention how much fun this is NOT! --jason)
 500
 501     my ($mixed,$mixedloc, $fstrand) = (0);
 502
 503     if( $self->isa('Bio::Das::SegmentI') &&
 504         ! $self->absolute ) {
 505         $self->warn("Calling spliced_seq with a Bio::Das::SegmentI which does have absolute set to 1 -- be warned you may not be getting things on the correct strand");
 506     }
 507
 508     my @locset = $self->location->each_Location;
 509     my @locs;
 510     if( ! $nosort ) {
 511         @locs = map { $_->[0] }
 512         # sort so that most negative is first basically to order
 513         # the features on the opposite strand 5'->3' on their strand
 514         # rather than they way most are input which is on the fwd strand
 515
 516         sort { $a->[1] <=> $b->[1] } # Yes Tim, Schwartzian transformation
 517         map {
 518             $fstrand = $_->strand unless defined $fstrand;
 519             $mixed = 1 if defined $_->strand && $fstrand != $_->strand;
 520             if( defined $_->seq_id ) {
 521                 $mixedloc = 1 if( $_->seq_id ne $seqid );
 522             }
 523             [ $_, $_->start * ($_->strand || 1)];
 524         } @locset;
 525
 526         if ( $mixed ) {
 527             $self->warn("Mixed strand locations, spliced seq using the input order rather than trying to sort");
 528             @locs = @locset;
 529         }
 530     } else {
 531         # use the original order instead of trying to sort
 532         @locs = @locset;
 533         $fstrand = $locs[0]->strand;
 534     }
 535
 536     foreach my $loc ( @locs ) {
 537         if( ! $loc->isa("Bio::Location::Atomic") ) {
 538             $self->throw("Can only deal with one level deep locations");
 539         }
 540         my $called_seq;
 541         if( $fstrand != $loc->strand ) {
 542             $self->warn("feature strand is different from location strand!");
 543         }
 544         # deal with remote sequences
 545
 546         if( defined $loc->seq_id &&
 547             $loc->seq_id ne $seqid ) {
 548             if( defined $db ) {
 549                 my $sid = $loc->seq_id;
 550                 $sid =~ s/\.\d+$//g;
 551                 eval {
 552                     $called_seq = $db->get_Seq_by_acc($sid);
 553                 };
 554                 if( $@ ) {
 555                     $self->warn("In attempting to join a remote location, sequence $sid was not in database. Will provide padding N's. Full exception \n\n$@");
 556                     $called_seq = undef;
 557                 }
 558             } else {
 559                 $self->warn( "cannot get remote location for ".$loc->seq_id ." without a valid Bio::DB::RandomAccessI database handle (like Bio::DB::GenBank)");
 560                 $called_seq = undef;
 561             }
 562             if( !defined $called_seq ) {
 563                 $seqstr .= 'N' x $self->length;
 564                 next;
 565             }
 566         } else {
 567             $called_seq = $self->entire_seq;
 568         }
 569
 570     # does the called sequence make sense? Bug 1780
 571     if ($called_seq->length < $loc->end) {
 572         my $accession = $called_seq->accession;
 573         my $end = $loc->end;
 574         my $length = $called_seq->length;
 575         my $orig_id = $self->seq_id; # originating sequence
 576         my ($locus) = $self->get_tagset_values("locus_tag");
 577         $self->throw("Location end ($end) exceeds length ($length) of ".
 578                      "called sequence $accession.\nCheck sequence version used in ".
 579                      "$locus locus-tagged SeqFeature in $orig_id.");
 580     }
 581
 582         if( $self->isa('Bio::Das::SegmentI') ) {
 583             my ($s,$e) = ($loc->start,$loc->end);
 584                 # $called_seq is Bio::DB::GFF::RelSegment, as well as its subseq();
 585                 # Bio::DB::GFF::RelSegment::seq() returns a Bio::PrimarySeq, and using seq()
 586                 # in turn returns a string.  Confused?
 587             $seqstr .= $called_seq->subseq($s,$e)->seq()->seq();
 588         } else {
 589             # This is dumb, subseq should work on locations...
 590             if( $loc->strand == 1 ) {
 591                 $seqstr .= $called_seq->subseq($loc->start,$loc->end);
 592             } else {
 593                 if( $nosort ) {
 594                     $seqstr = $called_seq->trunc($loc->start,$loc->end)->revcom->seq() . $seqstr;
 595                 } else {
 596                     $seqstr .= $called_seq->trunc($loc->start,$loc->end)->revcom->seq();
 597                 }
 598             }
 599         }
 600     }
 601
 602     if (defined($phase)) {
 603         $seqstr = substr($seqstr, $phase);
 604     }
 605
 606     my $out = Bio::Seq->new( -id => $self->entire_seq->display_id
 607                                     . "_spliced_feat",
 608                              -seq => $seqstr);
 609
 610     return $out;
 611 }
 612
 613 =head2 location
 614
 615  Title   : location
 616  Usage   : my $location = $seqfeature->location()
 617  Function: returns a location object suitable for identifying location
 618            of feature on sequence or parent feature
 619  Returns : Bio::LocationI object
 620  Args    : none
 621
 622
 623 =cut
 624
 625 sub location {
 626    my ($self) = @_;
 627
 628    $self->throw_not_implemented();
 629 }
 630
 631
 632 =head2 primary_id
 633
 634  Title   : primary_id
 635  Usage   : $obj->primary_id($newval)
 636  Function:
 637  Example :
 638  Returns : value of primary_id (a scalar)
 639  Args    : on set, new value (a scalar or undef, optional)
 640
 641 Primary ID is a synonym for the tag 'ID'
 642
 643 =cut
 644
 645 sub primary_id{
 646     my $self = shift;
 647     # note from cjm@fruitfly.org:
 648     # I have commented out the following 2 lines:
 649
 650     #return $self->{'primary_id'} = shift if @_;
 651     #return $self->{'primary_id'};
 652
 653     #... and replaced it with the following; see
 654     # http://bioperl.org/pipermail/bioperl-l/2003-December/014150.html
 655     # for the discussion that lead to this change
 656
 657     if (@_) {
 658         if ($self->has_tag('ID')) {
 659             $self->remove_tag('ID');
 660         }
 661         $self->add_tag_value('ID', shift);
 662     }
 663     my ($id) = $self->get_tagset_values('ID');
 664     return $id;
 665 }
 666
 667 sub generate_unique_persistent_id {
 668     # DEPRECATED - us IDHandler
 669     my $self = shift;
 670     require Bio::SeqFeature::Tools::IDHandler;
 671     Bio::SeqFeature::Tools::IDHandler->new->generate_unique_persistent_id($self);
 672 }
 673
 674 =head1 Bio::RangeI methods
 675
 676 These methods are inherited from RangeI and can be used
 677 directly from a SeqFeatureI interface. Remember that a
 678 SeqFeature is-a RangeI, and so wherever you see RangeI you
 679 can use a feature ($r in the below documentation).
 680
 681 =cut
 682
 683 =head2 start()
 684
 685  See L<Bio::RangeI>
 686
 687 =head2 end()
 688
 689  See L<Bio::RangeI>
 690
 691 =head2 strand()
 692
 693  See L<Bio::RangeI>
 694
 695 =head2 overlaps()
 696
 697  See L<Bio::RangeI>
 698
 699 =head2 contains()
 700
 701  See L<Bio::RangeI>
 702
 703 =head2 equals()
 704
 705  See L<Bio::RangeI>
 706
 707 =head2 intersection()
 708
 709  See L<Bio::RangeI>
 710
 711 =head2 union()
 712
 713  See L<Bio::RangeI>
 714
 715 =head1 Bio::AnnotatableI methods
 716
 717 =cut
 718
 719 =head2 has_tag()
 720
 721  B<Deprecated>.  See L<Bio::AnnotatableI>
 722
 723 =head2 remove_tag()
 724
 725  B<Deprecated>.  See L<Bio::AnnotatableI>
 726
 727 =head2 add_tag_value()
 728
 729  B<Deprecated>.  See L<Bio::AnnotatableI>
 730
 731 =head2 get_tag_values()
 732
 733  B<Deprecated>.  See L<Bio::AnnotatableI>
 734
 735 =head2 get_tagset_values()
 736
 737  B<Deprecated>.  See L<Bio::AnnotatableI>
 738
 739 =head2 get_all_tags()
 740
 741  B<Deprecated>.  See L<Bio::AnnotatableI>
 742
 743 =cut
 744
 745 1;