Bio/PrimarySeq.pm

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