Bio/PrimarySeq.pm

   1 # $Id$
   2 #
   3 # bioperl module for Bio::PrimarySeq
   4 #
   5 # Please direct questions and support issues to <bioperl-l@bioperl.org>
   6 #
   7 # Cared for by Ewan Birney <birney@ebi.ac.uk>
   8 #
   9 # Copyright Ewan Birney
  10 #
  11 # You may distribute this module under the same terms as perl itself
  12
  13 # POD documentation - main docs before the code
  14
  15 =head1 NAME
  16
  17 Bio::PrimarySeq - Bioperl lightweight Sequence Object
  18
  19 =head1 SYNOPSIS
  20
  21   # Bio::SeqIO for file reading, Bio::DB::GenBank for
  22   # database reading
  23
  24   use Bio::Seq;
  25   use Bio::SeqIO;
  26   use Bio::DB::GenBank;
  27
  28   # make from memory
  29
  30   $seqobj = Bio::PrimarySeq->new ( -seq => 'ATGGGGTGGGCGGTGGGTGGTTTG',
  31                                    -id  => 'GeneFragment-12',
  32                                    -accession_number => 'X78121',
  33                                    -alphabet => 'dna',
  34                                    -is_circular => 1 );
  35   print "Sequence ", $seqobj->id(), " with accession ",
  36     $seqobj->accession_number, "\n";
  37
  38   # read from file
  39
  40   $inputstream = Bio::SeqIO->new(-file => "myseq.fa",
  41                                  -format => 'Fasta');
  42   $seqobj = $inputstream->next_seq();
  43   print "Sequence ", $seqobj->id(), " and desc ", $seqobj->desc, "\n";
  44
  45   # to get out parts of the sequence.
  46
  47   print "Sequence ", $seqobj->id(), " with accession ",
  48     $seqobj->accession_number, " and desc ", $seqobj->desc, "\n";
  49
  50   $string  = $seqobj->seq();
  51   $string2 = $seqobj->subseq(1,40);
  52
  53 =head1 DESCRIPTION
  54
  55 PrimarySeq is a lightweight Sequence object, storing the sequence, its
  56 name, a computer-useful unique name, and other fundamental attributes.
  57 It does not contain sequence features or other information.  To have a
  58 sequence with sequence features you should use the Seq object which uses
  59 this object.
  60
  61 Although new users will use Bio::PrimarySeq a lot, in general you will
  62 be using it from the Bio::Seq object. For more information on Bio::Seq
  63 see L<Bio::Seq>. For interest you might like to know that
  64 Bio::Seq has-a Bio::PrimarySeq and forwards most of the function calls
  65 to do with sequence to it (the has-a relationship lets us get out of a
  66 otherwise nasty cyclical reference in Perl which would leak memory).
  67
  68 Sequence objects are defined by the Bio::PrimarySeqI interface, and this
  69 object is a pure Perl implementation of the interface. If that's
  70 gibberish to you, don't worry. The take home message is that this
  71 object is the bioperl default sequence object, but other people can
  72 use their own objects as sequences if they so wish. If you are
  73 interested in wrapping your own objects as compliant Bioperl sequence
  74 objects, then you should read the Bio::PrimarySeqI documentation
  75
  76 The documentation of this object is a merge of the Bio::PrimarySeq and
  77 Bio::PrimarySeqI documentation.  This allows all the methods which you can
  78 call on sequence objects here.
  79
  80 =head1 FEEDBACK
  81
  82 =head2 Mailing Lists
  83
  84 User feedback is an integral part of the evolution of this and other
  85 Bioperl modules. Send your comments and suggestions preferably to one
  86 of the Bioperl mailing lists.  Your participation is much appreciated.
  87
  88   bioperl-l@bioperl.org                  - General discussion
  89   http://bioperl.org/wiki/Mailing_lists  - About the mailing lists
  90
  91 =head2 Support
  92
  93 Please direct usage questions or support issues to the mailing list:
  94
  95 I<bioperl-l@bioperl.org>
  96
  97 rather than to the module maintainer directly. Many experienced and
  98 reponsive experts will be able look at the problem and quickly
  99 address it. Please include a thorough description of the problem
 100 with code and data examples if at all possible.
 101
 102 =head2 Reporting Bugs
 103
 104 Report bugs to the Bioperl bug tracking system to help us keep track
 105 the bugs and their resolution.  Bug reports can be submitted via the
 106 web:
 107
 108   http://bugzilla.open-bio.org/
 109
 110 =head1 AUTHOR - Ewan Birney
 111
 112 Email birney@ebi.ac.uk
 113
 114 =head1 APPENDIX
 115
 116 The rest of the documentation details each of the object
 117 methods. Internal methods are usually preceded with a _
 118
 119 =cut
 120
 121
 122 # Let the code begin...
 123
 124
 125 package Bio::PrimarySeq;
 126 use vars qw($MATCHPATTERN $GAP_SYMBOLS);
 127 use strict;
 128
 129 $MATCHPATTERN = 'A-Za-z\-\.\*\?=~';
 130 $GAP_SYMBOLS = '-~';
 131
 132 use base qw(Bio::Root::Root Bio::PrimarySeqI
 133             Bio::IdentifiableI Bio::DescribableI);
 134
 135 #
 136 # setup the allowed values for alphabet()
 137 #
 138
 139 my %valid_type = map {$_, 1} qw( dna rna protein );
 140
 141 =head2 new
 142
 143  Title   : new
 144  Usage   : $seq    = Bio::PrimarySeq->new( -seq => 'ATGGGGGTGGTGGTACCCT',
 145                                            -id  => 'human_id',
 146                                            -accession_number => 'AL000012',
 147                                            );
 148
 149  Function: Returns a new primary seq object from
 150            basic constructors, being a string for the sequence
 151            and strings for id and accession_number.
 152
 153            Note that you can provide an empty sequence string. However, in
 154            this case you MUST specify the type of sequence you wish to
 155            initialize by the parameter -alphabet. See alphabet() for possible
 156            values.
 157  Returns : a new Bio::PrimarySeq object
 158  Args    : -seq         => sequence string
 159            -display_id  => display id of the sequence (locus name)
 160            -accession_number => accession number
 161            -primary_id  => primary id (Genbank id)
 162            -version     => version number
 163            -namespace   => the namespace for the accession
 164            -authority   => the authority for the namespace
 165            -description => description text
 166            -desc        => alias for description
 167            -alphabet    => sequence type (alphabet) (dna|rna|protein)
 168            -id          => alias for display id
 169            -is_circular => boolean field for whether or not sequence is circular
 170            -direct      => boolean field for directly setting sequence (requires alphabet also set)
 171            -ref_to_seq  => boolean field indicating the sequence is a reference (?!?)
 172            -nowarnonempty => boolean field for whether or not to warn when sequence is empty
 173
 174 =cut
 175
 176
 177 sub new {
 178     my ($class, @args) = @_;
 179     my $self = $class->SUPER::new(@args);
 180
 181     my($seq,$id,$acc,$pid,$ns,$auth,$v,$oid,
 182        $desc,$description,
 183        $alphabet,$given_id,$is_circular,$direct,$ref_to_seq,$len,$nowarnonempty) =
 184         $self->_rearrange([qw(SEQ
 185                               DISPLAY_ID
 186                               ACCESSION_NUMBER
 187                               PRIMARY_ID
 188                               NAMESPACE
 189                               AUTHORITY
 190                               VERSION
 191                               OBJECT_ID
 192                               DESC
 193                               DESCRIPTION
 194                               ALPHABET
 195                               ID
 196                               IS_CIRCULAR
 197                               DIRECT
 198                               REF_TO_SEQ
 199                               LENGTH
 200             NOWARNONEMPTY
 201                               )],
 202                           @args);
 203
 204     # private var _nowarnonempty, need to be set before calling _guess_alphabet
 205     $self->{'_nowarnonempty'} = $nowarnonempty;
 206
 207     if( defined $id && defined $given_id ) {
 208       if( $id ne $given_id ) {
 209         $self->throw("Provided both id and display_id constructor ".
 210             "functions. [$id] [$given_id]");
 211       }
 212     }
 213     if( defined $given_id ) { $id = $given_id; }
 214
 215     # let's set the length before the seq -- if there is one, this length is
 216     # going to be invalidated
 217     defined $len && $self->length($len);
 218
 219     # if alphabet is provided we set it first, so that it won't be guessed
 220     # when the sequence is set
 221     $alphabet && $self->alphabet($alphabet);
 222
 223     # if there is an alphabet, and direct is passed in, assume the alphabet
 224     # and sequence is ok
 225
 226     if( $direct && $ref_to_seq) {
 227       $self->{'seq'} = $$ref_to_seq;
 228         if( ! $alphabet ) {
 229           $self->_guess_alphabet();
 230         } # else it has been set already above
 231     } else {
 232                  #      print STDERR "DEBUG: setting sequence to [$seq]\n";
 233                  # note: the sequence string may be empty
 234                  $self->seq($seq) if defined($seq);
 235          }
 236
 237     defined $id  && $self->display_id($id);
 238     $acc         && $self->accession_number($acc);
 239     defined $pid && $self->primary_id($pid);
 240     $desc        && $self->desc($desc);
 241     $description && $self->description($description);
 242     $is_circular && $self->is_circular($is_circular);
 243     $ns          && $self->namespace($ns);
 244     $auth        && $self->authority($auth);
 245     defined($v)  && $self->version($v);
 246     defined($oid) && $self->object_id($oid);
 247
 248
 249     return $self;
 250 }
 251
 252 sub direct_seq_set {
 253     my $obj = shift;
 254     return $obj->{'seq'} = shift if @_;
 255     return;
 256 }
 257
 258
 259 =head2 seq
 260
 261  Title   : seq
 262  Usage   : $string    = $obj->seq()
 263  Function: Returns the sequence as a string of letters. The
 264            case of the letters is left up to the implementer.
 265            Suggested cases are upper case for proteins and lower case for
 266            DNA sequence (IUPAC standard), but you should not rely on this.
 267  Returns : A scalar
 268  Args    : Optionally on set the new value (a string). An optional second
 269            argument presets the alphabet (otherwise it will be guessed).
 270
 271 =cut
 272
 273 sub seq {
 274    my ($obj,@args) = @_;
 275
 276    if( scalar(@args) == 0 ) {
 277        return $obj->{'seq'};
 278    }
 279
 280    my ($value,$alphabet) = @args;
 281
 282    if(@args) {
 283        if(defined($value) && (! $obj->validate_seq($value))) {
 284            $obj->throw("Attempting to set the sequence to [$value] ".
 285                                                         "which does not look healthy");
 286                 }
 287        # if a sequence was already set we make sure that we re-adjust the
 288        # alphabet, otherwise we skip guessing if alphabet is already set
 289        # note: if the new seq is empty or undef, we don't consider that a
 290        # change (we wouldn't have anything to guess on anyway)
 291                 my $is_changed_seq =
 292                   exists($obj->{'seq'}) && (CORE::length($value || '') > 0);
 293                 $obj->{'seq'} = $value;
 294        # new alphabet overridden by arguments?
 295                 if($alphabet) {
 296            # yes, set it no matter what
 297                         $obj->alphabet($alphabet);
 298                 } elsif( # if we changed a previous sequence to a new one
 299                                   $is_changed_seq ||
 300                                   # or if there is no alphabet yet at all
 301                                   (! defined($obj->alphabet()))) {
 302                         # we need to guess the (possibly new) alphabet
 303                         $obj->_guess_alphabet();
 304                 } # else (seq not changed and alphabet was defined) do nothing
 305                 # if the seq is changed, make sure we unset a possibly set length
 306                 $obj->length(undef) if $is_changed_seq || $obj->{'seq'};
 307    }
 308    return $obj->{'seq'};
 309 }
 310
 311 =head2 validate_seq
 312
 313  Title   : validate_seq
 314  Usage   : if(! $seq->validate_seq($seq_str) ) {
 315                 print "sequence $seq_str is not valid for an object of
 316                 alphabet ",$seq->alphabet, "\n";
 317            }
 318  Function: Validates a given sequence string. A validating sequence string
 319            must be accepted by seq(). A string that does not validate will
 320            lead to an exception if passed to seq().
 321
 322            The implementation provided here does not take alphabet() into
 323            account. Allowed are all letters (A-Z) and '-','.','*','?','=',
 324            and '~'.
 325
 326  Example :
 327  Returns : 1 if the supplied sequence string is valid for the object, and
 328            0 otherwise.
 329  Args    : The sequence string to be validated.
 330
 331
 332 =cut
 333
 334 sub validate_seq {
 335         my ($self,$seqstr) = @_;
 336         if( ! defined $seqstr ){ $seqstr = $self->seq(); }
 337         return 0 unless( defined $seqstr);
 338         if((CORE::length($seqstr) > 0) &&
 339            ($seqstr !~ /^([$MATCHPATTERN]+)$/)) {
 340             $self->warn("seq doesn't validate, mismatch is " .
 341                         join(",",($seqstr =~ /([^$MATCHPATTERN]+)/g)));
 342                 return 0;
 343         }
 344         return 1;
 345 }
 346
 347 =head2 subseq
 348
 349  Title   : subseq
 350  Usage   : $substring = $obj->subseq(10,40);
 351            $substring = $obj->subseq(10,40,NOGAP)
 352            $substring = $obj->subseq(-START=>10,-END=>40,-REPLACE_WITH=>'tga')
 353  Function: returns the subseq from start to end, where the first sequence
 354            character has coordinate 1 number is inclusive, ie 1-2 are the
 355            first two characters of the sequence
 356  Returns : a string
 357  Args    : integer for start position
 358            integer for end position
 359                  OR
 360            Bio::LocationI location for subseq (strand honored)
 361            Specify -NOGAP=>1 to return subseq with gap characters removed
 362            Specify -REPLACE_WITH=>$new_subseq to replace the subseq returned
 363            with $new_subseq in the sequence object
 364
 365 =cut
 366
 367 sub subseq {
 368    my $self = shift;
 369    my @args = @_;
 370    my ($start,$end,$nogap,$replace) = $self->_rearrange([qw(START
 371                                                             END
 372                                                             NOGAP
 373                                                             REPLACE_WITH)],@args);
 374
 375    # if $replace is specified, have the constructor validate it as seq
 376    my $dummy = new Bio::PrimarySeq(-seq=>$replace, -alphabet=>$self->alphabet) if defined($replace);
 377
 378    if( ref($start) && $start->isa('Bio::LocationI') ) {
 379        my $loc = $start;
 380        my $seq = "";
 381        foreach my $subloc ($loc->each_Location()) {
 382            my $piece = $self->subseq(-START=>$subloc->start(),
 383                                      '-END'=>$subloc->end(),
 384                                      -REPLACE_WITH=>$replace,
 385                                      -NOGAP=>$nogap);
 386            $piece =~ s/[$GAP_SYMBOLS]//g if $nogap;
 387            if($subloc->strand() < 0) {
 388                $piece = Bio::PrimarySeq->new('-seq' => $piece)->revcom()->seq();
 389            }
 390            $seq .= $piece;
 391        }
 392        return $seq;
 393    } elsif(  defined  $start && defined $end ) {
 394        if( $start > $end ){
 395            $self->throw("Bad start,end parameters. Start [$start] has to be ".
 396                         "less than end [$end]");
 397        }
 398        if( $start <= 0 ) {
 399            $self->throw("Bad start parameter ($start). Start must be positive.");
 400        }
 401        if( $end > $self->length ) {
 402            $self->throw("Bad end parameter ($end). End must be less than the total length of sequence (total=".$self->length.")");
 403        }
 404
 405        # remove one from start, and then length is end-start
 406        $start--;
 407        my @ss_args = map { eval "defined $_"  ? $_ : () } qw( $self->{seq} $start $end-$start $replace);
 408        my $seqstr = eval join( '', "substr(", join(',',@ss_args), ")");
 409        $seqstr =~ s/[$GAP_SYMBOLS]//g if ($nogap);
 410        return $seqstr;
 411
 412    } else {
 413        $self->warn("Incorrect parameters to subseq - must be two integers or a Bio::LocationI object. Got:", $self,$start,$end,$replace,$nogap);
 414        return;
 415    }
 416 }
 417
 418 =head2 length
 419
 420  Title   : length
 421  Usage   : $len = $seq->length();
 422  Function: Get the length of the sequence in number of symbols (bases
 423            or amino acids).
 424
 425            You can also set this attribute, even to a number that does
 426            not match the length of the sequence string. This is useful
 427            if you don''t want to set the sequence too, or if you want
 428            to free up memory by unsetting the sequence. In the latter
 429            case you could do e.g.
 430
 431                $seq->length($seq->length);
 432                $seq->seq(undef);
 433
 434            Note that if you set the sequence to a value other than
 435            undef at any time, the length attribute will be
 436            invalidated, and the length of the sequence string will be
 437            reported again. Also, we won''t let you lie about the length.
 438
 439  Example :
 440  Returns : integer representing the length of the sequence.
 441  Args    : Optionally, the value on set
 442
 443 =cut
 444
 445 sub length {
 446     my $self = shift;
 447     my $len = CORE::length($self->seq() || '');
 448
 449     if(@_) {
 450                  my $val = shift;
 451                  if(defined($val) && $len && ($len != $val)) {
 452                          $self->throw("You're trying to lie about the length: ".
 453                                                           "is $len but you say ".$val);
 454                  }
 455                  $self->{'_seq_length'} = $val;
 456     } elsif(defined($self->{'_seq_length'})) {
 457                  return $self->{'_seq_length'};
 458     }
 459     return $len;
 460 }
 461
 462 =head2 display_id
 463
 464  Title   : display_id or display_name
 465  Usage   : $id_string = $obj->display_id();
 466  Function: returns the display id, aka the common name of the Sequence object.
 467
 468            The semantics of this is that it is the most likely string to
 469            be used as an identifier of the sequence, and likely to have
 470            "human" readability.  The id is equivalent to the ID field of
 471            the GenBank/EMBL databanks and the id field of the
 472            Swissprot/sptrembl database. In fasta format, the >(\S+) is
 473            presumed to be the id, though some people overload the id to
 474            embed other information. Bioperl does not use any embedded
 475            information in the ID field, and people are encouraged to use
 476            other mechanisms (accession field for example, or extending
 477            the sequence object) to solve this.
 478
 479            With the new Bio::DescribeableI interface, display_name aliases
 480            to this method.
 481
 482  Returns : A string
 483  Args    : None
 484
 485
 486 =cut
 487
 488 sub display_id {
 489    my ($obj,$value) = @_;
 490    if( defined $value) {
 491       $obj->{'display_id'} = $value;
 492         }
 493         return $obj->{'display_id'};
 494 }
 495
 496 =head2 accession_number
 497
 498  Title   : accession_number or object_id
 499  Usage   : $unique_key = $obj->accession_number;
 500  Function: Returns the unique biological id for a sequence, commonly
 501            called the accession_number. For sequences from established
 502            databases, the implementors should try to use the correct
 503            accession number. Notice that primary_id() provides the
 504            unique id for the implemetation, allowing multiple objects
 505            to have the same accession number in a particular implementation.
 506
 507            For sequences with no accession number, this method should
 508            return "unknown".
 509
 510            [Note this method name is likely to change in 1.3]
 511
 512            With the new Bio::IdentifiableI interface, this is aliased
 513            to object_id
 514
 515  Returns : A string
 516  Args    : A string (optional) for setting
 517
 518 =cut
 519
 520 sub accession_number {
 521     my( $obj, $acc ) = @_;
 522
 523     if (defined $acc) {
 524                  $obj->{'accession_number'} = $acc;
 525     } else {
 526                  $acc = $obj->{'accession_number'};
 527                  $acc = 'unknown' unless defined $acc;
 528     }
 529     return $acc;
 530 }
 531
 532 =head2 primary_id
 533
 534  Title   : primary_id
 535  Usage   : $unique_key = $obj->primary_id;
 536  Function: Returns the unique id for this object in this
 537            implementation. This allows implementations to manage their
 538            own object ids in a way the implementaiton can control
 539            clients can expect one id to map to one object.
 540
 541            For sequences with no natural primary id, this method
 542            should return a stringified memory location.
 543
 544  Returns : A string
 545  Args    : A string (optional, for setting)
 546
 547 =cut
 548
 549 sub primary_id {
 550     my $obj = shift;
 551
 552     if(@_) {
 553                  $obj->{'primary_id'} = shift;
 554     }
 555     if( ! defined($obj->{'primary_id'}) ) {
 556                  return "$obj";
 557     }
 558     return $obj->{'primary_id'};
 559 }
 560
 561
 562 =head2 alphabet
 563
 564  Title   : alphabet
 565  Usage   : if( $obj->alphabet eq 'dna' ) { /Do Something/ }
 566  Function: Get/Set the alphabet of sequence, one of
 567            'dna', 'rna' or 'protein'. This is case sensitive.
 568
 569            This is not called <type> because this would cause
 570            upgrade problems from the 0.5 and earlier Seq objects.
 571
 572  Returns : a string either 'dna','rna','protein'. NB - the object must
 573            make a call of the type - if there is no alphabet specified it
 574            has to guess.
 575  Args    : optional string to set : 'dna' | 'rna' | 'protein'
 576
 577
 578 =cut
 579
 580 sub alphabet {
 581     my ($obj,$value) = @_;
 582     if (defined $value) {
 583                  $value = lc $value;
 584                  unless ( $valid_type{$value} ) {
 585                          $obj->throw("Alphabet '$value' is not a valid alphabet (".
 586                                                          join(',', map "'$_'", sort keys %valid_type) .
 587                                                          ") lowercase");
 588                  }
 589                  $obj->{'alphabet'} = $value;
 590     }
 591     return $obj->{'alphabet'};
 592 }
 593
 594 =head2 desc
 595
 596  Title   : desc or description
 597  Usage   : $obj->desc($newval)
 598  Function: Get/set description of the sequence.
 599
 600            'description' is an alias for this for compliance with the
 601            Bio::DescribeableI interface.
 602
 603  Example :
 604  Returns : value of desc (a string)
 605  Args    : newvalue (a string or undef, optional)
 606
 607
 608 =cut
 609
 610 sub desc{
 611     my $self = shift;
 612
 613     return $self->{'desc'} = shift if @_;
 614     return $self->{'desc'};
 615 }
 616
 617 =head2 can_call_new
 618
 619  Title   : can_call_new
 620  Usage   :
 621  Function:
 622  Example :
 623  Returns : true
 624  Args    :
 625
 626
 627 =cut
 628
 629 sub can_call_new {
 630    my ($self) = @_;
 631
 632    return 1;
 633 }
 634
 635 =head2 id
 636
 637  Title   : id
 638  Usage   : $id = $seq->id()
 639  Function: This is mapped on display_id
 640  Example :
 641  Returns :
 642  Args    :
 643
 644
 645 =cut
 646
 647 sub  id {
 648    return shift->display_id(@_);
 649 }
 650
 651 =head2 is_circular
 652
 653  Title   : is_circular
 654  Usage   : if( $obj->is_circular) { /Do Something/ }
 655  Function: Returns true if the molecule is circular
 656  Returns : Boolean value
 657  Args    : none
 658
 659 =cut
 660
 661 sub is_circular{
 662     my $self = shift;
 663
 664     return $self->{'is_circular'} = shift if @_;
 665     return $self->{'is_circular'};
 666 }
 667
 668
 669 =head1 Methods for Bio::IdentifiableI compliance
 670
 671 =cut
 672
 673 =head2 object_id
 674
 675  Title   : object_id
 676  Usage   : $string    = $obj->object_id()
 677  Function: A string which represents the stable primary identifier
 678            in this namespace of this object. For DNA sequences this
 679            is its accession_number, similarly for protein sequences.
 680
 681            This is aliased to accession_number().
 682  Returns : A scalar
 683
 684
 685 =cut
 686
 687 sub object_id {
 688     return shift->accession_number(@_);
 689 }
 690
 691 =head2 version
 692
 693  Title   : version
 694  Usage   : $version    = $obj->version()
 695  Function: A number which differentiates between versions of
 696            the same object. Higher numbers are considered to be
 697            later and more relevant, but a single object described
 698            the same identifier should represent the same concept.
 699
 700  Returns : A number
 701
 702 =cut
 703
 704 sub version{
 705     my ($self,$value) = @_;
 706     if( defined $value) {
 707                  $self->{'_version'} = $value;
 708     }
 709     return $self->{'_version'};
 710 }
 711
 712
 713 =head2 authority
 714
 715  Title   : authority
 716  Usage   : $authority    = $obj->authority()
 717  Function: A string which represents the organisation which
 718            granted the namespace, written as the DNS name for
 719            organisation (eg, wormbase.org).
 720
 721  Returns : A scalar
 722
 723 =cut
 724
 725 sub authority {
 726     my ($obj,$value) = @_;
 727     if( defined $value) {
 728                  $obj->{'authority'} = $value;
 729     }
 730     return $obj->{'authority'};
 731 }
 732
 733 =head2 namespace
 734
 735  Title   : namespace
 736  Usage   : $string    = $obj->namespace()
 737  Function: A string representing the name space this identifier
 738            is valid in, often the database name or the name
 739            describing the collection.
 740
 741  Returns : A scalar
 742
 743
 744 =cut
 745
 746 sub namespace{
 747     my ($self,$value) = @_;
 748     if( defined $value) {
 749                  $self->{'namespace'} = $value;
 750     }
 751     return $self->{'namespace'} || "";
 752 }
 753
 754 =head1 Methods for Bio::DescribableI compliance
 755
 756 This comprises of display_name and description.
 757
 758 =cut
 759
 760 =head2 display_name
 761
 762  Title   : display_name
 763  Usage   : $string    = $obj->display_name()
 764  Function: A string which is what should be displayed to the user.
 765            The string should have no spaces (ideally, though a cautious
 766            user of this interface would not assumme this) and should be
 767            less than thirty characters (though again, double checking
 768            this is a good idea).
 769
 770            This is aliased to display_id().
 771  Returns : A scalar
 772
 773 =cut
 774
 775 sub display_name {
 776     return shift->display_id(@_);
 777 }
 778
 779 =head2 description
 780
 781  Title   : description
 782  Usage   : $string    = $obj->description()
 783  Function: A text string suitable for displaying to the user a
 784            description. This string is likely to have spaces, but
 785            should not have any newlines or formatting - just plain
 786            text. The string should not be greater than 255 characters
 787            and clients can feel justified at truncating strings at 255
 788            characters for the purposes of display.
 789
 790            This is aliased to desc().
 791  Returns : A scalar
 792
 793 =cut
 794
 795 sub description {
 796     return shift->desc(@_);
 797 }
 798
 799 =head1 Methods Inherited from Bio::PrimarySeqI
 800
 801 These methods are available on Bio::PrimarySeq, although they are
 802 actually implemented on Bio::PrimarySeqI
 803
 804 =head2 revcom
 805
 806  Title   : revcom
 807  Usage   : $rev = $seq->revcom()
 808  Function: Produces a new Bio::SeqI implementing object which
 809            is the reversed complement of the sequence. For protein
 810            sequences this throws an exception of
 811            "Sequence is a protein. Cannot revcom".
 812
 813            The id is the same id as the orginal sequence, and the
 814            accession number is also indentical. If someone wants to
 815            track that this sequence has be reversed, it needs to
 816            define its own extensions.
 817
 818            To do an inplace edit of an object you can go:
 819
 820            $seqobj = $seqobj->revcom();
 821
 822            This of course, causes Perl to handle the garbage
 823            collection of the old object, but it is roughly speaking as
 824            efficient as an inplace edit.
 825
 826  Returns : A new (fresh) Bio::SeqI object
 827  Args    : none
 828
 829 =cut
 830
 831 =head2 trunc
 832
 833  Title   : trunc
 834  Usage   : $subseq = $myseq->trunc(10,100);
 835  Function: Provides a truncation of a sequence,
 836
 837  Example :
 838  Returns : A fresh Bio::SeqI implementing object.
 839  Args    :
 840
 841
 842 =cut
 843
 844 =head1 Internal methods
 845
 846 These are internal methods to PrimarySeq
 847
 848 =cut
 849
 850 =head2 _guess_alphabet
 851
 852  Title   : _guess_alphabet
 853  Usage   :
 854  Function: Determines (and sets) the type of sequence: dna, rna, protein
 855  Example :
 856  Returns : one of strings 'dna', 'rna' or 'protein'.
 857  Args    : none
 858
 859
 860 =cut
 861
 862 sub _guess_alphabet {
 863    my ($self) = @_;
 864    my $type;
 865
 866         #return if $self->alphabet;
 867
 868    my $str = $self->seq();
 869         # Remove char's that clearly denote ambiguity
 870    $str =~ s/[-.?]//gi;
 871
 872    my $total = CORE::length($str);
 873    if( $total == 0 ) {
 874      if (!$self->{'_nowarnonempty'}) {
 875        $self->warn("Got a sequence with no letters in it ".
 876            "cannot guess alphabet");
 877      }
 878      return '';
 879    }
 880
 881    my $u = ($str =~ tr/Uu//);
 882         # The assumption here is that most of sequences comprised of mainly
 883    # ATGC, with some N, will be 'dna' despite the fact that N could
 884         # also be Asparagine
 885    my $atgc = ($str =~ tr/ATGCNatgcn//);
 886
 887    if( ($atgc / $total) > 0.85 ) {
 888        $type = 'dna';
 889    } elsif( (($atgc + $u) / $total) > 0.85 ) {
 890        $type = 'rna';
 891    } else {
 892        $type = 'protein';
 893    }
 894
 895    $self->alphabet($type);
 896    return $type;
 897 }
 898
 899 ############################################################################
 900 # aliases due to name changes or to compensate for our lack of consistency #
 901 ############################################################################
 902
 903 sub accession {
 904     my $self = shift;
 905
 906     $self->warn(ref($self)."::accession is deprecated, ".
 907                 "use accession_number() instead");
 908     return $self->accession_number(@_);
 909 }
 910
 911 1;
 912