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