Bio/Seq.pm

   1 #
   2 # BioPerl module for Bio::Seq
   3 #
   4 # Please direct questions and support issues to <bioperl-l@bioperl.org>
   5 #
   6 # Cared for by Ewan Birney <birney@ebi.ac.uk>
   7 #
   8 # Copyright Ewan Birney
   9 #
  10 # You may distribute this module under the same terms as perl itself
  11
  12 # POD documentation - main docs before the code
  13
  14 =head1 NAME
  15
  16 Bio::Seq - Sequence object, with features
  17
  18 =head1 SYNOPSIS
  19
  20     # This is the main sequence object in Bioperl
  21
  22     # gets a sequence from a file
  23     $seqio  = Bio::SeqIO->new( '-format' => 'embl' , -file => 'myfile.dat');
  24     $seqobj = $seqio->next_seq();
  25
  26     # SeqIO can both read and write sequences; see Bio::SeqIO
  27     # for more information and examples
  28
  29     # get from database
  30     $db = Bio::DB::GenBank->new();
  31     $seqobj = $db->get_Seq_by_acc('X78121');
  32
  33     # make from strings in script
  34     $seqobj = Bio::Seq->new( -display_id => 'my_id',
  35                              -seq => $sequence_as_string);
  36
  37     # gets sequence as a string from sequence object
  38     $seqstr   = $seqobj->seq(); # actual sequence as a string
  39     $seqstr   = $seqobj->subseq(10,50); # slice in biological coordinates
  40
  41     # retrieves information from the sequence
  42     # features must implement Bio::SeqFeatureI interface
  43
  44     @features = $seqobj->get_SeqFeatures(); # just top level
  45     foreach my $feat ( @features ) {
  46         print "Feature ",$feat->primary_tag," starts ",$feat->start," ends ",
  47         $feat->end," strand ",$feat->strand,"\n";
  48
  49         # features retain link to underlying sequence object
  50         print "Feature sequence is ",$feat->seq->seq(),"\n"
  51     }
  52
  53     # sequences may have a species
  54
  55     if( defined $seq->species ) {
  56         print "Sequence is from ",$species->binomial," [",$species->common_name,"]\n";
  57     }
  58
  59     # annotation objects are Bio::AnnotationCollectionI's
  60     $ann      = $seqobj->annotation(); # annotation object
  61
  62     # references is one type of annotations to get. Also get
  63     # comment and dblink. Look at Bio::AnnotationCollection for
  64     # more information
  65
  66     foreach my $ref ( $ann->get_Annotations('reference') ) {
  67         print "Reference ",$ref->title,"\n";
  68     }
  69
  70     # you can get truncations, translations and reverse complements, these
  71     # all give back Bio::Seq objects themselves, though currently with no
  72     # features transfered
  73
  74     my $trunc = $seqobj->trunc(100,200);
  75     my $rev   = $seqobj->revcom();
  76
  77     # there are many options to translate - check out the docs
  78     my $trans = $seqobj->translate();
  79
  80     # these functions can be chained together
  81
  82     my $trans_trunc_rev = $seqobj->trunc(100,200)->revcom->translate();
  83
  84
  85
  86 =head1 DESCRIPTION
  87
  88 A Seq object is a sequence with sequence features placed on it. The
  89 Seq object contains a PrimarySeq object for the actual sequence and
  90 also implements its interface.
  91
  92 In Bioperl we have 3 main players that people are going to use frequently
  93
  94   Bio::PrimarySeq  - just the sequence and its names, nothing else.
  95   Bio::SeqFeatureI - a feature on a sequence, potentially with a sequence
  96                      and a location and annotation.
  97   Bio::Seq         - A sequence and a collection of sequence features
  98                      (an aggregate) with its own annotation.
  99
 100 Although Bioperl is not tied heavily to file formats these distinctions do
 101 map to file formats sensibly and for some bioinformaticians this might help
 102
 103   Bio::PrimarySeq  - Fasta file of a sequence
 104   Bio::SeqFeatureI - A single entry in an EMBL/GenBank/DDBJ feature table
 105   Bio::Seq         - A single EMBL/GenBank/DDBJ entry
 106
 107 By having this split we avoid a lot of nasty circular references
 108 (sequence features can hold a reference to a sequence without the sequence
 109 holding a reference to the sequence feature). See L<Bio::PrimarySeq> and
 110 L<Bio::SeqFeatureI> for more information.
 111
 112 Ian Korf really helped in the design of the Seq and SeqFeature system.
 113
 114 =head2 Examples
 115
 116 A simple and fundamental block of code:
 117
 118   use Bio::SeqIO;
 119
 120   my $seqIOobj = Bio::SeqIO->new(-file=>"1.fa"); # create a SeqIO object
 121   my $seqobj = $seqIOobj->next_seq;              # get a Seq object
 122
 123 With the Seq object in hand one has access to a powerful set of Bioperl
 124 methods and related Bioperl objects. This next script will take a file of sequences
 125 in EMBL format and create a file of the reverse-complemented sequences
 126 in Fasta format using Seq objects. It also prints out details about the
 127 exons it finds as sequence features in Genbank Flat File format.
 128
 129   use Bio::Seq;
 130   use Bio::SeqIO;
 131
 132   $seqin = Bio::SeqIO->new( -format => 'EMBL' , -file => 'myfile.dat');
 133   $seqout= Bio::SeqIO->new( -format => 'Fasta', -file => '>output.fa');
 134
 135   while((my $seqobj = $seqin->next_seq())) {
 136       print "Seen sequence ",$seqobj->display_id,", start of seq ",
 137             substr($seqobj->seq,1,10),"\n";
 138       if( $seqobj->alphabet eq 'dna') {
 139             $rev = $seqobj->revcom;
 140             $id  = $seqobj->display_id();
 141             $id  = "$id.rev";
 142             $rev->display_id($id);
 143             $seqout->write_seq($rev);
 144       }
 145
 146       foreach $feat ( $seqobj->get_SeqFeatures() ) {
 147            if( $feat->primary_tag eq 'exon' ) {
 148               print STDOUT "Location ",$feat->start,":",
 149                     $feat->end," GFF[",$feat->gff_string,"]\n";
 150            }
 151       }
 152   }
 153
 154 Let's examine the script. The lines below import the Bioperl modules.
 155 Seq is the main Bioperl sequence object and SeqIO is the Bioperl support
 156 for reading sequences from files and to files
 157
 158   use Bio::Seq;
 159   use Bio::SeqIO;
 160
 161 These two lines create two SeqIO streams: one for reading in sequences
 162 and one for outputting sequences:
 163
 164   $seqin = Bio::SeqIO->new( -format => 'EMBL' , -file => 'myfile.dat');
 165   $seqout= Bio::SeqIO->new( -format => 'Fasta', -file => '>output.fa');
 166
 167 Notice that in the "$seqout" case there is a greater-than sign,
 168 indicating the file is being opened for writing.
 169
 170 Using the
 171
 172   '-argument' => value
 173
 174 syntax is common in Bioperl. The file argument is like an argument
 175 to open() . You can also pass in filehandles or FileHandle objects by
 176 using the -fh argument (see L<Bio::SeqIO> documentation for details).
 177 Many formats in Bioperl are handled, including Fasta, EMBL, GenBank,
 178 Swissprot (swiss), PIR, and GCG.
 179
 180   $seqin = Bio::SeqIO->new( -format => 'EMBL' , -file => 'myfile.dat');
 181   $seqout= Bio::SeqIO->new( -format => 'Fasta', -file => '>output.fa');
 182
 183 This is the main loop which will loop progressively through sequences
 184 in a file, and each call to $seqio-E<gt>next_seq() provides a new Seq
 185 object from the file:
 186
 187   while((my $seqobj = $seqio->next_seq())) {
 188
 189 This print line below accesses fields in the Seq object directly. The
 190 $seqobj-E<gt>display_id is the way to access the display_id attribute
 191 of the Seq object. The $seqobj-E<gt>seq method gets the actual
 192 sequence out as string. Then you can do manipulation of this if
 193 you want to (there are however easy ways of doing truncation,
 194 reverse-complement and translation).
 195
 196   print "Seen sequence ",$seqobj->display_id,", start of seq ",
 197                substr($seqobj->seq,1,10),"\n";
 198
 199 Bioperl has to guess the alphabet of the sequence, being either 'dna',
 200 'rna', or 'protein'. The alphabet attribute is one of these three
 201 possibilities.
 202
 203   if( $seqobj->alphabet eq 'dna') {
 204
 205 The $seqobj-E<gt>revcom method provides the reverse complement of the Seq
 206 object as another Seq object. Thus, the $rev variable is a reference to
 207 another Seq object. For example, one could repeat the above print line
 208 for this Seq object (putting $rev in place of $seqobj). In this
 209 case we are going to output the object into the file stream we built
 210 earlier on.
 211
 212   $rev = $seqobj->revcom;
 213
 214 When we output it, we want the id of the outputted object
 215 to be changed to "$id.rev", ie, with .rev on the end of the name. The
 216 following lines retrieve the id of the sequence object, add .rev
 217 to this and then set the display_id of the rev sequence object to
 218 this. Notice that to set the display_id attribute you just need
 219 call the same method, display_id(), with the new value as an argument.
 220 Getting and setting values with the same method is common in Bioperl.
 221
 222   $id  = $seqobj->display_id();
 223   $id  = "$id.rev";
 224   $rev->display_id($id);
 225
 226 The write_seq method on the SeqIO output object, $seqout, writes the
 227 $rev object to the filestream we built at the top of the script.
 228 The filestream knows that it is outputting in fasta format, and
 229 so it provides fasta output.
 230
 231   $seqout->write_seq($rev);
 232
 233 This block of code loops over sequence features in the sequence
 234 object, trying to find ones who have been tagged as 'exon'.
 235 Features have start and end attributes and can be outputted
 236 in Genbank Flat File format, GFF, a standarized format for sequence
 237 features.
 238
 239   foreach $feat ( $seqobj->get_SeqFeatures() ) {
 240       if( $feat->primary_tag eq 'exon' ) {
 241           print STDOUT "Location ",$feat->start,":",
 242              $feat->end," GFF[",$feat->gff_string,"]\n";
 243       }
 244   }
 245
 246 The code above shows how a few Bio::Seq methods suffice to read, parse,
 247 reformat and analyze sequences from a file. A full list of methods
 248 available to Bio::Seq objects is shown below. Bear in mind that some of
 249 these methods come from PrimarySeq objects, which are simpler
 250 than Seq objects, stripped of features (see L<Bio::PrimarySeq> for
 251 more information).
 252
 253   # these methods return strings, and accept strings in some cases:
 254
 255   $seqobj->seq();              # string of sequence
 256   $seqobj->subseq(5,10);       # part of the sequence as a string
 257   $seqobj->accession_number(); # when there, the accession number
 258   $seqobj->alphabet();         # one of 'dna','rna',or 'protein'
 259   $seqobj->version()           # when there, the version
 260   $seqobj->keywords();         # when there, the Keywords line
 261   $seqobj->length()            # length
 262   $seqobj->desc();             # description
 263   $seqobj->primary_id();       # a unique id for this sequence regardless
 264                                # of its display_id or accession number
 265   $seqobj->display_id();       # the human readable id of the sequence
 266
 267 Some of these values map to fields in common formats. For example, The
 268 display_id() method returns the LOCUS name of a Genbank entry,
 269 the (\S+) following the E<gt> character in a Fasta file, the ID from
 270 a SwissProt file, and so on. The desc() method will return the DEFINITION
 271 line of a Genbank file, the description following the display_id in a
 272 Fasta file, and the DE field in a SwissProt file.
 273
 274   # the following methods return new Seq objects, but
 275   # do not transfer features across to the new object:
 276
 277   $seqobj->trunc(5,10)  # truncation from 5 to 10 as new object
 278   $seqobj->revcom       # reverse complements sequence
 279   $seqobj->translate    # translation of the sequence
 280
 281   # if new() can be called this method returns 1, else 0
 282
 283   $seqobj->can_call_new
 284
 285   # the following method determines if the given string will be accepted
 286   # by the seq() method - if the string is acceptable then validate()
 287   # returns 1, or 0 if not
 288
 289   $seqobj->validate_seq($string)
 290
 291   # the following method returns or accepts a Species object:
 292
 293   $seqobj->species();
 294
 295 Please see L<Bio::Species> for more information on this object.
 296
 297   # the following method returns or accepts an Annotation object
 298   # which in turn allows access to Annotation::Reference
 299   # and Annotation::Comment objects:
 300
 301   $seqobj->annotation();
 302
 303 These annotations typically refer to entire sequences, unlike
 304 features.  See L<Bio::AnnotationCollectionI>,
 305 L<Bio::Annotation::Collection>, L<Bio::Annotation::Reference>, and
 306 L<Bio::Annotation::Comment> for details.
 307
 308 It is also important to be able to describe defined portions of a
 309 sequence. The combination of some description and the corresponding
 310 sub-sequence is called a feature - an exon and its coordinates within
 311 a gene is an example of a feature, or a domain within a protein.
 312
 313   # the following methods return an array of SeqFeatureI objects:
 314
 315   $seqobj->get_SeqFeatures # The 'top level' sequence features
 316   $seqobj->get_all_SeqFeatures # All sequence features, including sub-seq
 317                                # features, such as features in an exon
 318
 319   # to find out the number of features use:
 320
 321   $seqobj->feature_count
 322
 323 Here are just some of the methods available to SeqFeatureI objects:
 324
 325   # these methods return numbers:
 326
 327   $feat->start          # start position (1 is the first base)
 328   $feat->end            # end position (2 is the second base)
 329   $feat->strand         # 1 means forward, -1 reverse, 0 not relevant
 330
 331   # these methods return or accept strings:
 332
 333   $feat->primary_tag    # the name of the sequence feature, eg
 334                         # 'exon', 'glycoslyation site', 'TM domain'
 335   $feat->source_tag     # where the feature comes from, eg, 'EMBL_GenBank',
 336                         # or 'BLAST'
 337
 338   # this method returns the more austere PrimarySeq object, not a
 339   # Seq object - the main difference is that PrimarySeq objects do not
 340   # themselves contain sequence features
 341
 342   $feat->seq            # the sequence between start,end on the
 343                         # correct strand of the sequence
 344
 345 See L<Bio::PrimarySeq> for more details on PrimarySeq objects.
 346
 347   # useful methods for feature comparisons, for start/end points
 348
 349   $feat->overlaps($other)  # do $feat and $other overlap?
 350   $feat->contains($other)  # is $other completely within $feat?
 351   $feat->equals($other)    # do $feat and $other completely agree?
 352
 353   # one can also add features
 354
 355   $seqobj->add_SeqFeature($feat)     # returns 1 if successful
 356   $seqobj->add_SeqFeature(@features) # returns 1 if successful
 357
 358   # sub features. For complex join() statements, the feature
 359   # is one sequence feature with many sub SeqFeatures
 360
 361   $feat->sub_SeqFeature  # returns array of sub seq features
 362
 363 Please see L<Bio::SeqFeatureI> and L<Bio::SeqFeature::Generic>,
 364 for more information on sequence features.
 365
 366 It is worth mentioning that one can also retrieve the start and end
 367 positions of a feature using a Bio::LocationI object:
 368
 369   $location = $feat->location # $location is a Bio::LocationI object
 370   $location->start;           # start position
 371   $location->end;             # end position
 372
 373 This is useful because one needs a Bio::Location::SplitLocationI object
 374 in order to retrieve the coordinates inside the Genbank or EMBL join()
 375 statements (e.g. "CDS    join(51..142,273..495,1346..1474)"):
 376
 377   if ( $feat->location->isa('Bio::Location::SplitLocationI') &&
 378                $feat->primary_tag eq 'CDS' )  {
 379     foreach $loc ( $feat->location->sub_Location ) {
 380       print $loc->start . ".." . $loc->end . "\n";
 381     }
 382   }
 383
 384 See L<Bio::LocationI> and L<Bio::Location::SplitLocationI> for more
 385 information.
 386
 387 =head1 Implemented Interfaces
 388
 389 This class implements the following interfaces.
 390
 391 =over 4
 392
 393 =item Bio::SeqI
 394
 395 Note that this includes implementing Bio::PrimarySeqI.
 396
 397 =item Bio::IdentifiableI
 398
 399 =item Bio::DescribableI
 400
 401 =item Bio::AnnotatableI
 402
 403 =item Bio::FeatureHolderI
 404
 405 =back
 406
 407 =head1 FEEDBACK
 408
 409
 410 =head2 Mailing Lists
 411
 412 User feedback is an integral part of the evolution of this and other
 413 Bioperl modules. Send your comments and suggestions preferably to one
 414 of the Bioperl mailing lists. Your participation is much appreciated.
 415
 416   bioperl-l@bioperl.org                  - General discussion
 417   http://bioperl.org/wiki/Mailing_lists  - About the mailing lists
 418
 419 =head2 Support
 420
 421 Please direct usage questions or support issues to the mailing list:
 422
 423 I<bioperl-l@bioperl.org>
 424
 425 rather than to the module maintainer directly. Many experienced and
 426 reponsive experts will be able look at the problem and quickly
 427 address it. Please include a thorough description of the problem
 428 with code and data examples if at all possible.
 429
 430 =head2 Reporting Bugs
 431
 432 Report bugs to the Bioperl bug tracking system to help us keep track
 433 the bugs and their resolution.  Bug reports can be submitted via the
 434 web:
 435
 436   https://redmine.open-bio.org/projects/bioperl/
 437
 438 =head1 AUTHOR - Ewan Birney, inspired by Ian Korf objects
 439
 440 Email birney@ebi.ac.uk
 441
 442 =head1 CONTRIBUTORS
 443
 444  Jason Stajich E<lt>jason@bioperl.orgE<gt>
 445  Mark A. Jensen maj -at- fortinbras -dot- us
 446
 447 =head1 APPENDIX
 448
 449
 450 The rest of the documentation details each of the object
 451 methods. Internal methods are usually preceded with a "_".
 452
 453 =cut
 454
 455 #'
 456 # Let the code begin...
 457
 458
 459 package Bio::Seq;
 460 use strict;
 461
 462 use Bio::Annotation::Collection;
 463 use Bio::PrimarySeq;
 464
 465 use base qw(Bio::Root::Root Bio::SeqI Bio::IdentifiableI Bio::DescribableI Bio::AnnotatableI Bio::FeatureHolderI Bio::AnnotationCollectionI);
 466
 467 =head2 new
 468
 469  Title   : new
 470  Usage   : $seq = Bio::Seq->new( -seq => 'ATGGGGGTGGTGGTACCCT',
 471                                  -id  => 'human_id',
 472                                  -accession_number => 'AL000012',
 473                                );
 474
 475  Function: Returns a new Seq object from
 476            basic constructors, being a string for the sequence
 477            and strings for id and accession_number
 478  Returns : a new Bio::Seq object
 479
 480 =cut
 481
 482 sub new {
 483     my($caller,@args) = @_;
 484
 485     if( $caller ne 'Bio::Seq') {
 486         $caller = ref($caller) if ref($caller);
 487     }
 488
 489     # we know our inherietance hierarchy
 490     my $self = Bio::Root::Root->new(@args);
 491     bless $self,$caller;
 492
 493     # this is way too sneaky probably. We delegate the construction of
 494     # the Seq object onto PrimarySeq and then pop primary_seq into
 495     # our primary_seq slot
 496
 497     my $pseq = Bio::PrimarySeq->new(@args);
 498
 499     # as we have just made this, we know it is ok to set hash directly
 500     # rather than going through the method
 501
 502     $self->{'primary_seq'} = $pseq;
 503
 504     # setting this array is now delayed until the final
 505     # moment, again speed ups for non feature containing things
 506     # $self->{'_as_feat'} = [];
 507
 508
 509     my ($ann, $pid,$feat,$species) = &Bio::Root::RootI::_rearrange($self,[qw(ANNOTATION PRIMARY_ID FEATURES SPECIES)], @args);
 510
 511     # for a number of cases - reading fasta files - these are never set. This
 512     # gives a quick optimisation around testing things later on
 513
 514     if( defined $ann || defined $pid || defined $feat || defined $species ) {
 515         $pid && $self->primary_id($pid);
 516         $species && $self->species($species);
 517         $ann && $self->annotation($ann);
 518
 519         if( defined $feat ) {
 520             if( ref($feat) !~ /ARRAY/i ) {
 521                 if( ref($feat) && $feat->isa('Bio::SeqFeatureI') ) {
 522                     $self->add_SeqFeature($feat);
 523                 } else {
 524                     $self->warn("Must specify a valid Bio::SeqFeatureI or ArrayRef of Bio::SeqFeatureI's with the -features init parameter for ".ref($self));
 525                 }
 526             } else {
 527                 foreach my $feature ( @$feat ) {
 528                     $self->add_SeqFeature($feature);
 529                 }
 530             }
 531         }
 532     }
 533
 534     return $self;
 535 }
 536
 537 =head1 PrimarySeq interface
 538
 539
 540 The PrimarySeq interface provides the basic sequence getting
 541 and setting methods for on all sequences.
 542
 543 These methods implement the Bio::PrimarySeq interface by delegating
 544 to the primary_seq inside the object. This means that you
 545 can use a Seq object wherever there is a PrimarySeq, and
 546 of course, you are free to use these functions anyway.
 547
 548 =cut
 549
 550 =head2 seq
 551
 552  Title   : seq
 553  Usage   : $string = $obj->seq()
 554  Function: Get/Set the sequence as a string of letters. The
 555            case of the letters is left up to the implementer.
 556            Suggested cases are upper case for proteins and lower case for
 557            DNA sequence (IUPAC standard),
 558            but implementations are suggested to keep an open mind about
 559            case (some users... want mixed case!)
 560  Returns : A scalar
 561  Args    : Optionally on set the new value (a string). An optional second
 562            argument presets the alphabet (otherwise it will be guessed).
 563            Both parameters may also be given in named parameter style
 564            with -seq and -alphabet being the names.
 565
 566 =cut
 567
 568 sub seq {
 569     return shift->primary_seq()->seq(@_);
 570 }
 571
 572
 573 =head2 validate_seq
 574
 575  Title   : validate_seq
 576  Usage   : if(! $seq->validate_seq($seq_str) ) {
 577                 print "sequence $seq_str is not valid for an object of type ",
 578                       ref($seq), "\n";
 579            }
 580  Function: Validates a given sequence string. A validating sequence string
 581            must be accepted by seq(). A string that does not validate will
 582            lead to an exception if passed to seq().
 583
 584            The implementation provided here does not take alphabet() into
 585            account. Allowed are all letters (A-Z), '-','.','*','=', and '~'.
 586
 587  Example :
 588  Returns : 1 if the supplied sequence string is valid for the object, and
 589            0 otherwise.
 590  Args    : The sequence string to be validated.
 591
 592
 593 =cut
 594
 595 sub validate_seq {
 596     return shift->primary_seq()->validate_seq(@_);
 597 }
 598
 599 =head2 length
 600
 601  Title   : length
 602  Usage   : $len = $seq->length()
 603  Function:
 604  Example :
 605  Returns : Integer representing the length of the sequence.
 606  Args    : None
 607
 608 =cut
 609
 610 sub length {
 611     return shift->primary_seq()->length(@_);
 612 }
 613
 614 =head1 Methods from the Bio::PrimarySeqI interface
 615
 616 =cut
 617
 618 =head2 subseq
 619
 620  Title   : subseq
 621  Usage   : $substring = $obj->subseq(10,40);
 622  Function: Returns the subseq from start to end, where the first base
 623            is 1 and the number is inclusive, ie 1-2 are the first two
 624            bases of the sequence
 625
 626            Start cannot be larger than end but can be equal
 627
 628  Returns : A string
 629  Args    : 2 integers
 630
 631
 632 =cut
 633
 634 sub subseq {
 635     return shift->primary_seq()->subseq(@_);
 636 }
 637
 638 =head2 display_id
 639
 640  Title   : display_id
 641  Usage   : $id = $obj->display_id or $obj->display_id($newid);
 642  Function: Gets or sets the display id, also known as the common name of
 643            the Seq object.
 644
 645            The semantics of this is that it is the most likely string
 646            to be used as an identifier of the sequence, and likely to
 647            have "human" readability.  The id is equivalent to the LOCUS
 648            field of the GenBank/EMBL databanks and the ID field of the
 649            Swissprot/sptrembl database. In fasta format, the >(\S+) is
 650            presumed to be the id, though some people overload the id
 651            to embed other information. Bioperl does not use any
 652            embedded information in the ID field, and people are
 653            encouraged to use other mechanisms (accession field for
 654            example, or extending the sequence object) to solve this.
 655
 656            Notice that $seq->id() maps to this function, mainly for
 657            legacy/convenience issues.
 658  Returns : A string
 659  Args    : None or a new id
 660
 661
 662 =cut
 663
 664 sub display_id {
 665    return shift->primary_seq->display_id(@_);
 666 }
 667
 668
 669
 670 =head2 accession_number
 671
 672  Title   : accession_number
 673  Usage   : $unique_biological_key = $obj->accession_number;
 674  Function: Returns the unique biological id for a sequence, commonly
 675            called the accession_number. For sequences from established
 676            databases, the implementors should try to use the correct
 677            accession number. Notice that primary_id() provides the
 678            unique id for the implemetation, allowing multiple objects
 679            to have the same accession number in a particular implementation.
 680
 681            For sequences with no accession number, this method should return
 682            "unknown".
 683
 684            Can also be used to set the accession number.
 685  Example : $key = $seq->accession_number or $seq->accession_number($key)
 686  Returns : A string
 687  Args    : None or an accession number
 688
 689
 690 =cut
 691
 692 sub accession_number {
 693    return shift->primary_seq->accession_number(@_);
 694 }
 695
 696 =head2 desc
 697
 698  Title   : desc
 699  Usage   : $seqobj->desc($string) or $seqobj->desc()
 700  Function: Sets or gets the description of the sequence
 701  Example :
 702  Returns : The description
 703  Args    : The description or none
 704
 705
 706 =cut
 707
 708 sub desc {
 709    return shift->primary_seq->desc(@_);
 710 }
 711
 712 =head2 primary_id
 713
 714  Title   : primary_id
 715  Usage   : $unique_implementation_key = $obj->primary_id;
 716  Function: Returns the unique id for this object in this
 717            implementation. This allows implementations to manage
 718            their own object ids in a way the implementation can control
 719            clients can expect one id to map to one object.
 720
 721            For sequences with no natural id, this method should return
 722            a stringified memory location.
 723
 724            Can also be used to set the primary_id (or unset to undef).
 725
 726            [Note this method name is likely to change in 1.3]
 727
 728  Example : $id = $seq->primary_id or $seq->primary_id($id)
 729  Returns : A string
 730  Args    : None or an id, or undef to unset the primary id.
 731
 732
 733 =cut
 734
 735 sub primary_id {
 736     # Note: this used to not delegate to the primary seq. This is
 737     # really bad in very subtle ways. E.g., if you created the object
 738     # with a primary id given to the constructor and then later you
 739     # change the primary id, if this method wouldn't delegate you'd
 740     # have different values for primary id in the PrimarySeq object
 741     # compared to this instance. Not good.
 742
 743     # I can't remember why not delegating was ever deemed
 744     # advantageous, but I hereby claim that its problems far outweigh
 745     # its advantages, if there are any. Convince me otherwise if you
 746     # disagree. HL 2004/08/05
 747
 748     return shift->primary_seq->primary_id(@_);
 749 }
 750
 751 =head2 can_call_new
 752
 753  Title   : can_call_new
 754  Usage   : if ( $obj->can_call_new ) {
 755              $newobj = $obj->new( %param );
 756            }
 757  Function: can_call_new returns 1 or 0 depending
 758            on whether an implementation allows new
 759            constructor to be called. If a new constructor
 760            is allowed, then it should take the followed hashed
 761            constructor list.
 762
 763            $myobject->new( -seq => $sequence_as_string,
 764                            -display_id  => $id
 765                            -accession_number => $accession
 766                            -alphabet => 'dna',
 767                          );
 768  Example :
 769  Returns : 1 or 0
 770  Args    : None
 771
 772
 773 =cut
 774
 775 sub can_call_new {
 776     return 1;
 777 }
 778
 779 =head2 alphabet
 780
 781  Title   : alphabet
 782  Usage   : if ( $obj->alphabet eq 'dna' ) { /Do Something/ }
 783  Function: Get/Set the type of sequence being one of
 784            'dna', 'rna' or 'protein'. This is case sensitive.
 785
 786            This is not called <type> because this would cause
 787            upgrade problems from the 0.5 and earlier Seq objects.
 788
 789  Returns : A string either 'dna','rna','protein'. NB - the object must
 790            make a call of the type - if there is no type specified it
 791            has to guess.
 792  Args    : optional string to set : 'dna' | 'rna' | 'protein'
 793
 794
 795 =cut
 796
 797 sub alphabet {
 798    my $self = shift;
 799    return $self->primary_seq->alphabet(@_) if @_ && defined $_[0];
 800    return $self->primary_seq->alphabet();
 801 }
 802
 803 =head2 is_circular
 804
 805  Title   : is_circular
 806  Usage   : if( $obj->is_circular) { /Do Something/ }
 807  Function: Returns true if the molecule is circular
 808  Returns : Boolean value
 809  Args    : none
 810
 811 =cut
 812
 813 sub is_circular {
 814     return shift->primary_seq()->is_circular(@_);
 815 }
 816
 817
 818 =head1 Methods for Bio::IdentifiableI compliance
 819
 820 =cut
 821
 822 =head2 object_id
 823
 824  Title   : object_id
 825  Usage   : $string    = $obj->object_id()
 826  Function: a string which represents the stable primary identifier
 827            in this namespace of this object. For DNA sequences this
 828            is its accession_number, similarly for protein sequences
 829
 830            This is aliased to accession_number().
 831  Returns : A scalar
 832
 833
 834 =cut
 835
 836 sub object_id {
 837     return shift->accession_number(@_);
 838 }
 839
 840 =head2 version
 841
 842  Title   : version
 843  Usage   : $version    = $obj->version()
 844  Function: a number which differentiates between versions of
 845            the same object. Higher numbers are considered to be
 846            later and more relevant, but a single object described
 847            the same identifier should represent the same concept
 848
 849  Returns : A number
 850
 851 =cut
 852
 853 sub version{
 854     return shift->primary_seq->version(@_);
 855 }
 856
 857
 858 =head2 authority
 859
 860  Title   : authority
 861  Usage   : $authority    = $obj->authority()
 862  Function: a string which represents the organisation which
 863            granted the namespace, written as the DNS name for
 864            organisation (eg, wormbase.org)
 865
 866  Returns : A scalar
 867
 868 =cut
 869
 870 sub authority {
 871     return shift->primary_seq()->authority(@_);
 872 }
 873
 874 =head2 namespace
 875
 876  Title   : namespace
 877  Usage   : $string    = $obj->namespace()
 878  Function: A string representing the name space this identifier
 879            is valid in, often the database name or the name
 880            describing the collection
 881
 882  Returns : A scalar
 883
 884
 885 =cut
 886
 887 sub namespace{
 888     return shift->primary_seq()->namespace(@_);
 889 }
 890
 891 =head1 Methods for Bio::DescribableI compliance
 892
 893 =cut
 894
 895 =head2 display_name
 896
 897  Title   : display_name
 898  Usage   : $string    = $obj->display_name()
 899  Function: A string which is what should be displayed to the user
 900            the string should have no spaces (ideally, though a cautious
 901            user of this interface would not assumme this) and should be
 902            less than thirty characters (though again, double checking
 903            this is a good idea)
 904
 905            This is aliased to display_id().
 906  Returns : A scalar
 907
 908 =cut
 909
 910 sub display_name {
 911     return shift->display_id(@_);
 912 }
 913
 914 =head2 description
 915
 916  Title   : description
 917  Usage   : $string    = $obj->description()
 918  Function: A text string suitable for displaying to the user a
 919            description. This string is likely to have spaces, but
 920            should not have any newlines or formatting - just plain
 921            text. The string should not be greater than 255 characters
 922            and clients can feel justified at truncating strings at 255
 923            characters for the purposes of display
 924
 925            This is aliased to desc().
 926  Returns : A scalar
 927
 928 =cut
 929
 930 sub description {
 931     return shift->desc(@_);
 932 }
 933
 934 =head1 Methods for implementing Bio::AnnotatableI
 935
 936 =cut
 937
 938 =head2 annotation
 939
 940  Title   : annotation
 941  Usage   : $ann = $seq->annotation or
 942            $seq->annotation($ann)
 943  Function: Gets or sets the annotation
 944  Returns : Bio::AnnotationCollectionI object
 945  Args    : None or Bio::AnnotationCollectionI object
 946
 947 See L<Bio::AnnotationCollectionI> and L<Bio::Annotation::Collection>
 948 for more information
 949
 950 =cut
 951
 952 sub annotation {
 953     my ($obj,$value) = @_;
 954     if( defined $value ) {
 955         $obj->throw("object of class ".ref($value)." does not implement ".
 956                     "Bio::AnnotationCollectionI. Too bad.")
 957             unless $value->isa("Bio::AnnotationCollectionI");
 958         $obj->{'_annotation'} = $value;
 959     } elsif( ! defined $obj->{'_annotation'}) {
 960         $obj->{'_annotation'} = Bio::Annotation::Collection->new();
 961     }
 962     return $obj->{'_annotation'};
 963 }
 964
 965 =head1 Methods for delegating Bio::AnnotationCollectionI
 966
 967 =head2 get_Annotations()
 968
 969  Usage   : my @annotations = $seq->get_Annotations('key')
 970  Function: Retrieves all the Bio::AnnotationI objects for a specific key
 971            for this object
 972  Returns : list of Bio::AnnotationI - empty if no objects stored for a key
 973  Args    : string which is key for annotations
 974
 975 =cut
 976
 977 sub get_Annotations { shift->annotation->get_Annotations(@_); }
 978
 979 =head2 add_Annotation()
 980
 981  Usage   : $seq->add_Annotation('reference',$object);
 982            $seq->add_Annotation($object,'Bio::MyInterface::DiseaseI');
 983            $seq->add_Annotation($object);
 984            $seq->add_Annotation('disease',$object,'Bio::MyInterface::DiseaseI');
 985  Function: Adds an annotation for a specific key for this sequence object.
 986
 987            If the key is omitted, the object to be added must provide a value
 988            via its tagname().
 989
 990            If the archetype is provided, this and future objects added under
 991            that tag have to comply with the archetype and will be rejected
 992            otherwise.
 993
 994  Returns : none
 995  Args    : annotation key ('disease', 'dblink', ...)
 996            object to store (must be Bio::AnnotationI compliant)
 997            [optional] object archetype to map future storage of object
 998            of these types to
 999
1000 =cut
1001
1002 sub add_Annotation { shift->annotation->add_Annotation(@_) }
1003
1004 =head2 remove_Annotations()
1005
1006  Usage   : $seq->remove_Annotations()
1007  Function: Remove the annotations for the specified key from this sequence
1008            object
1009  Returns : an list of Bio::AnnotationI compliant objects which were stored
1010            under the given key(s) for this sequence object
1011  Args    : the key(s) (tag name(s), one or more strings) for which to
1012            remove annotations (optional; if none given, flushes all
1013            annotations)
1014
1015 =cut
1016
1017 sub remove_Annotations { shift->annotation->remove_Annotations(@_) }
1018
1019 =head2 get_num_of_annotations()
1020
1021  Usage   : my $count = $seq->get_num_of_annotations()
1022  Alias   : num_Annotations
1023  Function: Returns the count of all annotations stored for this sequence
1024            object
1025  Returns : integer
1026  Args    : none
1027
1028 =cut
1029
1030 sub get_num_of_annotations { shift->annotation->get_num_of_annotations(@_) }
1031 sub num_Annotations { shift->get_num_of_annotations }; #DWYM
1032
1033 =head1 Methods to implement Bio::FeatureHolderI
1034
1035 This includes methods for retrieving, adding, and removing features.
1036
1037 =cut
1038
1039 =head2 get_SeqFeatures
1040
1041  Title   : get_SeqFeatures
1042  Usage   :
1043  Function: Get the feature objects held by this feature holder.
1044
1045            Features which are not top-level are subfeatures of one or
1046            more of the returned feature objects, which means that you
1047            must traverse the subfeature arrays of each top-level
1048            feature object in order to traverse all features associated
1049            with this sequence.
1050
1051            Top-level features can be obtained by tag, specified in
1052            the argument.
1053
1054            Use get_all_SeqFeatures() if you want the feature tree
1055            flattened into one single array.
1056
1057  Example :
1058  Returns : an array of Bio::SeqFeatureI implementing objects
1059  Args    : [optional] scalar string (feature tag)
1060
1061
1062 =cut
1063
1064 sub get_SeqFeatures{
1065    my $self = shift;
1066    my $tag = shift;
1067
1068    if( !defined $self->{'_as_feat'} ) {
1069        $self->{'_as_feat'} = [];
1070    }
1071    if ($tag) {
1072        return map { $_->primary_tag eq $tag ? $_ : () } @{$self->{'_as_feat'}};
1073    }
1074    else {
1075        return @{$self->{'_as_feat'}};
1076    }
1077 }
1078
1079 =head2 get_all_SeqFeatures
1080
1081  Title   : get_all_SeqFeatures
1082  Usage   : @feat_ary = $seq->get_all_SeqFeatures();
1083  Function: Returns the tree of feature objects attached to this
1084            sequence object flattened into one single array. Top-level
1085            features will still contain their subfeature-arrays, which
1086            means that you will encounter subfeatures twice if you
1087            traverse the subfeature tree of the returned objects.
1088
1089            Use get_SeqFeatures() if you want the array to contain only
1090            the top-level features.
1091
1092  Returns : An array of Bio::SeqFeatureI implementing objects.
1093  Args    : None
1094
1095
1096 =cut
1097
1098 # this implementation is inherited from FeatureHolderI
1099
1100 =head2 feature_count
1101
1102  Title   : feature_count
1103  Usage   : $seq->feature_count()
1104  Function: Return the number of SeqFeatures attached to a sequence
1105  Returns : integer representing the number of SeqFeatures
1106  Args    : None
1107
1108
1109 =cut
1110
1111 sub feature_count {
1112     my ($self) = @_;
1113
1114     if (defined($self->{'_as_feat'})) {
1115         return ($#{$self->{'_as_feat'}} + 1);
1116     } else {
1117         return 0;
1118     }
1119 }
1120
1121 =head2 add_SeqFeature
1122
1123  Title   : add_SeqFeature
1124  Usage   : $seq->add_SeqFeature($feat);
1125            $seq->add_SeqFeature(@feat);
1126  Function: Adds the given feature object (or each of an array of feature
1127            objects to the feature array of this
1128            sequence. The object passed is required to implement the
1129            Bio::SeqFeatureI interface.
1130  Returns : 1 on success
1131  Args    : A Bio::SeqFeatureI implementing object, or an array of such objects.
1132
1133
1134 =cut
1135
1136 sub add_SeqFeature {
1137    my ($self,@feat) = @_;
1138
1139    $self->{'_as_feat'} = [] unless $self->{'_as_feat'};
1140
1141    foreach my $feat ( @feat ) {
1142        if( !$feat->isa("Bio::SeqFeatureI") ) {
1143            $self->throw("$feat is not a SeqFeatureI and that's what we expect...");
1144        }
1145
1146        # make sure we attach ourselves to the feature if the feature wants it
1147        my $aseq = $self->primary_seq;
1148        $feat->attach_seq($aseq) if $aseq;
1149
1150        push(@{$self->{'_as_feat'}},$feat);
1151    }
1152    return 1;
1153 }
1154
1155 =head2 remove_SeqFeatures
1156
1157  Title   : remove_SeqFeatures
1158  Usage   : $seq->remove_SeqFeatures();
1159  Function: Flushes all attached SeqFeatureI objects.
1160
1161            To remove individual feature objects, delete those from the returned
1162            array and re-add the rest.
1163  Example :
1164  Returns : The array of Bio::SeqFeatureI objects removed from this seq.
1165  Args    : None
1166
1167
1168 =cut
1169
1170 sub remove_SeqFeatures {
1171     my $self = shift;
1172
1173     return () unless $self->{'_as_feat'};
1174     my @feats = @{$self->{'_as_feat'}};
1175     $self->{'_as_feat'} = [];
1176     return @feats;
1177 }
1178
1179 =head1 Methods provided in the Bio::PrimarySeqI interface
1180
1181
1182 These methods are inherited from the PrimarySeq interface
1183 and work as one expects, building new Bio::Seq objects
1184 or other information as expected. See L<Bio::PrimarySeq>
1185 for more information.
1186
1187 Sequence Features are B<not> transfered to the new objects.
1188 This is possibly a mistake. Anyone who feels the urge in
1189 dealing with this is welcome to give it a go.
1190
1191 =head2 revcom
1192
1193  Title   : revcom
1194  Usage   : $rev = $seq->revcom()
1195  Function: Produces a new Bio::Seq object which
1196            is the reversed complement of the sequence. For protein
1197            sequences this throws an exception of "Sequence is a protein.
1198            Cannot revcom"
1199
1200            The id is the same id as the original sequence, and the
1201            accession number is also identical. If someone wants to track
1202            that this sequence has be reversed, it needs to define its own
1203            extensions
1204
1205            To do an in-place edit of an object you can go:
1206
1207            $seq = $seq->revcom();
1208
1209            This of course, causes Perl to handle the garbage collection of
1210            the old object, but it is roughly speaking as efficient as an
1211            in-place edit.
1212
1213  Returns : A new (fresh) Bio::Seq object
1214  Args    : None
1215
1216
1217 =cut
1218
1219 =head2 trunc
1220
1221  Title   : trunc
1222  Usage   : $subseq = $myseq->trunc(10,100);
1223  Function: Provides a truncation of a sequence
1224
1225  Example :
1226  Returns : A fresh Seq object
1227  Args    : A Seq object
1228
1229
1230 =cut
1231
1232 =head2 id
1233
1234  Title   : id
1235  Usage   : $id = $seq->id()
1236  Function: This is mapped on display_id
1237  Returns : value of display_id()
1238  Args    : [optional] value to update display_id
1239
1240
1241 =cut
1242
1243 sub  id {
1244     return shift->display_id(@_);
1245 }
1246
1247
1248 =head1 Seq only methods
1249
1250
1251 These methods are specific to the Bio::Seq object, and not
1252 found on the Bio::PrimarySeq object
1253
1254 =head2 primary_seq
1255
1256  Title   : primary_seq
1257  Usage   : $seq->primary_seq or $seq->primary_seq($newval)
1258  Function: Get or set a PrimarySeq object
1259  Example :
1260  Returns : PrimarySeq object
1261  Args    : None or PrimarySeq object
1262
1263
1264 =cut
1265
1266 sub primary_seq {
1267    my ($obj,$value) = @_;
1268
1269    if( defined $value) {
1270        if( ! ref $value || ! $value->isa('Bio::PrimarySeqI') ) {
1271            $obj->throw("$value is not a Bio::PrimarySeq compliant object");
1272        }
1273
1274        $obj->{'primary_seq'} = $value;
1275        # descend down over all seqfeature objects, seeing whether they
1276        # want an attached seq.
1277
1278        foreach my $sf ( $obj->get_SeqFeatures() ) {
1279            $sf->attach_seq($value);
1280        }
1281
1282    }
1283    return $obj->{'primary_seq'};
1284
1285 }
1286
1287 =head2 species
1288
1289  Title   : species
1290  Usage   : $species = $seq->species() or $seq->species($species)
1291  Function: Gets or sets the species
1292  Returns : L<Bio::Species> object
1293  Args    : None or L<Bio::Species> object
1294
1295 See L<Bio::Species> for more information
1296
1297 =cut
1298
1299 sub species {
1300     my ($self, $species) = @_;
1301     if ($species) {
1302         $self->{'species'} = $species;
1303     } else {
1304         return $self->{'species'};
1305     }
1306 }
1307
1308 =head1 Internal methods
1309
1310 =cut
1311
1312 # keep AUTOLOAD happy
1313 sub DESTROY { }
1314
1315 ############################################################################
1316 # aliases due to name changes or to compensate for our lack of consistency #
1317 ############################################################################
1318
1319 # in all other modules we use the object in the singular --
1320 # lack of consistency sucks
1321 *flush_SeqFeature = \&remove_SeqFeatures;
1322 *flush_SeqFeatures = \&remove_SeqFeatures;
1323
1324 # this is now get_SeqFeatures() (from FeatureHolderI)
1325 *top_SeqFeatures = \&get_SeqFeatures;
1326
1327 # this is now get_all_SeqFeatures() in FeatureHolderI
1328 sub all_SeqFeatures{
1329     return shift->get_all_SeqFeatures(@_);
1330 }
1331
1332 sub accession {
1333     my $self = shift;
1334     $self->warn(ref($self)."::accession is deprecated, ".
1335                 "use accession_number() instead");
1336     return $self->accession_number(@_);
1337 }
1338
1339 1;