Bio/Seq.pm

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