2 # bioperl module for Bio::PrimarySeq
4 # Please direct questions and support issues to <bioperl-l@bioperl.org>
6 # Cared for by Ewan Birney <birney@ebi.ac.uk>
8 # Copyright Ewan Birney
10 # You may distribute this module under the same terms as perl itself
12 # POD documentation - main docs before the code
16 Bio::PrimarySeq - Bioperl lightweight Sequence Object
20 # Bio::SeqIO for file reading, Bio::DB::GenBank for
29 $seqobj = Bio::PrimarySeq->new ( -seq => 'ATGGGGTGGGCGGTGGGTGGTTTG',
30 -id => 'GeneFragment-12',
31 -accession_number => 'X78121',
34 print "Sequence ", $seqobj->id(), " with accession ",
35 $seqobj->accession_number, "\n";
39 $inputstream = Bio::SeqIO->new(-file => "myseq.fa",
41 $seqobj = $inputstream->next_seq();
42 print "Sequence ", $seqobj->id(), " and desc ", $seqobj->desc, "\n";
44 # to get out parts of the sequence.
46 print "Sequence ", $seqobj->id(), " with accession ",
47 $seqobj->accession_number, " and desc ", $seqobj->desc, "\n";
49 $string = $seqobj->seq();
50 $string2 = $seqobj->subseq(1,40);
54 PrimarySeq is a lightweight Sequence object, storing the sequence, its
55 name, a computer-useful unique name, and other fundamental attributes.
56 It does not contain sequence features or other information. To have a
57 sequence with sequence features you should use the Seq object which uses
60 Although new users will use Bio::PrimarySeq a lot, in general you will
61 be using it from the Bio::Seq object. For more information on Bio::Seq
62 see L<Bio::Seq>. For interest you might like to know that
63 Bio::Seq has-a Bio::PrimarySeq and forwards most of the function calls
64 to do with sequence to it (the has-a relationship lets us get out of a
65 otherwise nasty cyclical reference in Perl which would leak memory).
67 Sequence objects are defined by the Bio::PrimarySeqI interface, and this
68 object is a pure Perl implementation of the interface. If that's
69 gibberish to you, don't worry. The take home message is that this
70 object is the bioperl default sequence object, but other people can
71 use their own objects as sequences if they so wish. If you are
72 interested in wrapping your own objects as compliant Bioperl sequence
73 objects, then you should read the Bio::PrimarySeqI documentation
75 The documentation of this object is a merge of the Bio::PrimarySeq and
76 Bio::PrimarySeqI documentation. This allows all the methods which you can
77 call on sequence objects here.
83 User feedback is an integral part of the evolution of this and other
84 Bioperl modules. Send your comments and suggestions preferably to one
85 of the Bioperl mailing lists. Your participation is much appreciated.
87 bioperl-l@bioperl.org - General discussion
88 http://bioperl.org/wiki/Mailing_lists - About the mailing lists
92 Please direct usage questions or support issues to the mailing list:
94 I<bioperl-l@bioperl.org>
96 rather than to the module maintainer directly. Many experienced and
97 reponsive experts will be able look at the problem and quickly
98 address it. Please include a thorough description of the problem
99 with code and data examples if at all possible.
101 =head2 Reporting Bugs
103 Report bugs to the Bioperl bug tracking system to help us keep track
104 the bugs and their resolution. Bug reports can be submitted via the
107 https://redmine.open-bio.org/projects/bioperl/
109 =head1 AUTHOR - Ewan Birney
111 Email birney@ebi.ac.uk
115 The rest of the documentation details each of the object
116 methods. Internal methods are usually preceded with a _
121 # Let the code begin...
124 package Bio
::PrimarySeq
;
125 use vars
qw($MATCHPATTERN $GAP_SYMBOLS);
128 $MATCHPATTERN = 'A-Za-z\-\.\*\?=~';
131 use base qw(Bio::Root::Root Bio::PrimarySeqI
132 Bio::IdentifiableI Bio::DescribableI);
135 # setup the allowed values for alphabet()
138 my %valid_type = map {$_, 1} qw( dna rna protein );
143 Usage : $seq = Bio::PrimarySeq->new( -seq => 'ATGGGGGTGGTGGTACCCT',
145 -accession_number => 'AL000012',
148 Function: Returns a new primary seq object from
149 basic constructors, being a string for the sequence
150 and strings for id and accession_number.
152 Note that you can provide an empty sequence string. However, in
153 this case you MUST specify the type of sequence you wish to
154 initialize by the parameter -alphabet. See alphabet() for possible
156 Returns : a new Bio::PrimarySeq object
157 Args : -seq => sequence string
158 -display_id => display id of the sequence (locus name)
159 -accession_number => accession number
160 -primary_id => primary id (Genbank id)
161 -version => version number
162 -namespace => the namespace for the accession
163 -authority => the authority for the namespace
164 -description => description text
165 -desc => alias for description
166 -alphabet => sequence type (alphabet) (dna|rna|protein)
167 -id => alias for display id
168 -is_circular => boolean field for whether or not sequence is circular
169 -direct => boolean field for directly setting sequence (requires alphabet also set)
170 -ref_to_seq => boolean field indicating the sequence is a reference (?!?)
171 -nowarnonempty => boolean field for whether or not to warn when sequence is empty
177 my ($class, @args) = @_;
178 my $self = $class->SUPER::new
(@args);
180 my($seq,$id,$acc,$pid,$ns,$auth,$v,$oid,
182 $alphabet,$given_id,$is_circular,$direct,$ref_to_seq,$len,$nowarnonempty) =
183 $self->_rearrange([qw(SEQ
203 # private var _nowarnonempty, need to be set before calling _guess_alphabet
204 $self->{'_nowarnonempty'} = $nowarnonempty;
206 if( defined $id && defined $given_id ) {
207 if( $id ne $given_id ) {
208 $self->throw("Provided both id and display_id constructor ".
209 "functions. [$id] [$given_id]");
212 if( defined $given_id ) { $id = $given_id; }
214 # let's set the length before the seq -- if there is one, this length is
215 # going to be invalidated
216 defined $len && $self->length($len);
218 # if alphabet is provided we set it first, so that it won't be guessed
219 # when the sequence is set
220 $alphabet && $self->alphabet($alphabet);
222 # bernd's idea: define ids so that invalid sequence messages
223 # can be more informative...
224 defined $id && $self->display_id($id);
225 $acc && $self->accession_number($acc);
226 defined $pid && $self->primary_id($pid);
228 # if there is an alphabet, and direct is passed in, assume the alphabet
231 if( $direct && $ref_to_seq) {
232 $self->{'seq'} = $$ref_to_seq;
234 $self->_guess_alphabet();
235 } # else it has been set already above
237 # print STDERR "DEBUG: setting sequence to [$seq]\n";
238 # note: the sequence string may be empty
239 $self->seq($seq) if defined($seq);
242 $desc && $self->desc($desc);
243 $description && $self->description($description);
244 $is_circular && $self->is_circular($is_circular);
245 $ns && $self->namespace($ns);
246 $auth && $self->authority($auth);
247 defined($v) && $self->version($v);
248 defined($oid) && $self->object_id($oid);
256 return $obj->{'seq'} = shift if @_;
264 Usage : $string = $obj->seq()
265 Function: Returns the sequence as a string of letters. The
266 case of the letters is left up to the implementer.
267 Suggested cases are upper case for proteins and lower case for
268 DNA sequence (IUPAC standard), but you should not rely on this.
270 Args : Optionally on set the new value (a string). An optional second
271 argument presets the alphabet (otherwise it will be guessed).
276 my ($obj,@args) = @_;
278 if( scalar(@args) == 0 ) {
279 return $obj->{'seq'};
282 my ($value,$alphabet) = @args;
285 if(defined($value) && (! $obj->validate_seq($value))) {
286 $obj->throw("Attempting to set the sequence '".(defined($obj->id) ||
287 "[unidentified sequence]")."' to [$value] which does not look healthy");
289 # if a sequence was already set we make sure that we re-adjust the
290 # alphabet, otherwise we skip guessing if alphabet is already set
291 # note: if the new seq is empty or undef, we don't consider that a
292 # change (we wouldn't have anything to guess on anyway)
294 exists($obj->{'seq'}) && (CORE
::length($value || '') > 0);
295 $obj->{'seq'} = $value;
296 # new alphabet overridden by arguments?
298 # yes, set it no matter what
299 $obj->alphabet($alphabet);
300 } elsif ($is_changed_seq || (! defined($obj->alphabet()))) {
301 # if we changed a previous sequence to a new one or if there is no
302 # alphabet yet at all, 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'};
308 return $obj->{'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";
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().
322 The implementation provided here does not take alphabet() into
323 account. Allowed are all letters (A-Z) and '-','.','*','?','=',
327 Returns : 1 if the supplied sequence string is valid for the object, and
329 Args : The sequence string to be validated.
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("sequence '".(defined($self->id) || "[unidentified sequence]").
341 "' doesn't validate, mismatch is " .
342 join(",",($seqstr =~ /([^$MATCHPATTERN]+)/g)));
351 Usage : $substring = $obj->subseq(10,40);
352 $substring = $obj->subseq(10,40,NOGAP)
353 $substring = $obj->subseq(-START=>10,-END=>40,-REPLACE_WITH=>'tga')
354 Function: returns the subseq from start to end, where the first sequence
355 character has coordinate 1 number is inclusive, ie 1-2 are the
356 first two characters of the sequence
358 Args : integer for start position
359 integer for end position
361 Bio::LocationI location for subseq (strand honored)
362 Specify -NOGAP=>1 to return subseq with gap characters removed
363 Specify -REPLACE_WITH=>$new_subseq to replace the subseq returned
364 with $new_subseq in the sequence object
371 my ($start,$end,$nogap,$replace) = $self->_rearrange([qw(START
374 REPLACE_WITH)],@args);
376 # if $replace is specified, have the constructor validate it as seq
377 my $dummy = new Bio
::PrimarySeq
(-seq
=>$replace, -alphabet
=>$self->alphabet) if defined($replace);
379 if( ref($start) && $start->isa('Bio::LocationI') ) {
382 foreach my $subloc ($loc->each_Location()) {
383 my $piece = $self->subseq(-START
=>$subloc->start(),
384 '-END'=>$subloc->end(),
385 -REPLACE_WITH
=>$replace,
387 $piece =~ s/[$GAP_SYMBOLS]//g if $nogap;
388 if($subloc->strand() < 0) {
389 $piece = Bio
::PrimarySeq
->new('-seq' => $piece)->revcom()->seq();
394 } elsif( defined $start && defined $end ) {
396 $self->throw("Bad start,end parameters. Start [$start] has to be ".
397 "less than end [$end]");
400 $self->throw("Bad start parameter ($start). Start must be positive.");
403 # remove one from start, and then length is end-start
405 my @ss_args = map { eval "defined $_" ?
$_ : () } qw( $self->{seq} $start $end-$start $replace);
406 my $seqstr = eval join( '', "substr(", join(',',@ss_args), ")");
408 if( $end > $self->length) {
409 if ($self->is_circular) {
411 my $end = $end - $self->length;
412 my @ss_args = map { eval "defined $_" ? $_ : () } qw( $self->{seq} $start $end-$start $replace);
413 my $appendstr = eval join( '', "substr(", join(',',@ss_args), ")");
414 $seqstr .= $appendstr;
416 $self->throw("Bad end parameter ($end). End must be less than the total length of sequence (total=".$self->length.")")
420 $seqstr =~ s/[$GAP_SYMBOLS]//g if ($nogap);
424 $self->warn("Incorrect parameters to subseq - must be two integers or a Bio::LocationI object. Got:", $self,$start,$end,$replace,$nogap);
432 Usage : $len = $seq->length();
433 Function: Get the length of the sequence in number of symbols (bases
436 You can also set this attribute, even to a number that does
437 not match the length of the sequence string. This is useful
438 if you don''t want to set the sequence too, or if you want
439 to free up memory by unsetting the sequence. In the latter
440 case you could do e.g.
442 $seq->length($seq->length);
445 Note that if you set the sequence to a value other than
446 undef at any time, the length attribute will be
447 invalidated, and the length of the sequence string will be
448 reported again. Also, we won''t let you lie about the length.
451 Returns : integer representing the length of the sequence.
452 Args : Optionally, the value on set
458 my $len = CORE::length($self->seq() || '');
462 if(defined($val) && $len && ($len != $val)) {
463 $self->throw("You're trying to lie about the length: ".
464 "is $len but you say ".$val);
466 $self->{'_seq_length'} = $val;
467 } elsif(defined($self->{'_seq_length'})) {
468 return $self->{'_seq_length'};
475 Title : display_id or display_name
476 Usage : $id_string = $obj->display_id();
477 Function: returns the display id, aka the common name of the Sequence object.
479 The semantics of this is that it is the most likely string to
480 be used as an identifier of the sequence, and likely to have
481 "human" readability. The id is equivalent to the ID field of
482 the GenBank/EMBL databanks and the id field of the
483 Swissprot/sptrembl database. In fasta format, the >(\S+) is
484 presumed to be the id, though some people overload the id to
485 embed other information. Bioperl does not use any embedded
486 information in the ID field, and people are encouraged to use
487 other mechanisms (accession field for example, or extending
488 the sequence object) to solve this.
490 With the new Bio::DescribeableI interface, display_name aliases
500 my ($obj,$value) = @_;
501 if( defined $value) {
502 $obj->{'display_id'} = $value;
504 return $obj->{'display_id'};
507 =head2 accession_number
509 Title : accession_number or object_id
510 Usage : $unique_key = $obj->accession_number;
511 Function: Returns the unique biological id for a sequence, commonly
512 called the accession_number. For sequences from established
513 databases, the implementors should try to use the correct
514 accession number. Notice that primary_id() provides the
515 unique id for the implemetation, allowing multiple objects
516 to have the same accession number in a particular implementation.
518 For sequences with no accession number, this method should
521 [Note this method name is likely to change in 1.3]
523 With the new Bio::IdentifiableI interface, this is aliased
527 Args : A string (optional) for setting
531 sub accession_number {
532 my( $obj, $acc ) = @_;
535 $obj->{'accession_number'} = $acc;
537 $acc = $obj->{'accession_number'};
538 $acc = 'unknown' unless defined $acc;
546 Usage : $unique_key = $obj->primary_id;
547 Function: Returns the unique id for this object in this
548 implementation. This allows implementations to manage their
549 own object ids in a way the implementaiton can control
550 clients can expect one id to map to one object.
552 For sequences with no natural primary id, this method
553 should return a stringified memory location.
556 Args : A string (optional, for setting)
564 $obj->{'primary_id'} = shift;
566 if( ! defined($obj->{'primary_id'}) ) {
569 return $obj->{'primary_id'};
576 Usage : if( $obj->alphabet eq 'dna' ) { /Do Something/ }
577 Function: Get/Set the alphabet of sequence, one of
578 'dna', 'rna' or 'protein'. This is case sensitive.
580 This is not called <type> because this would cause
581 upgrade problems from the 0.5 and earlier Seq objects.
583 Returns : a string either 'dna','rna','protein'. NB - the object must
584 make a call of the type - if there is no alphabet specified it
586 Args : optional string to set : 'dna' | 'rna' | 'protein'
592 my ($obj,$value) = @_;
593 if (defined $value) {
595 unless ( $valid_type{$value} ) {
596 $obj->throw("Alphabet '$value' is not a valid alphabet (".
597 join(',', map "'$_'", sort keys %valid_type) .") lowercase");
599 $obj->{'alphabet'} = $value;
601 return $obj->{'alphabet'};
606 Title : desc or description
607 Usage : $obj->desc($newval)
608 Function: Get/set description of the sequence.
610 'description' is an alias for this for compliance with the
611 Bio::DescribeableI interface.
614 Returns : value of desc (a string)
615 Args : newvalue (a string or undef, optional)
623 return $self->{'desc'} = shift if @_;
624 return $self->{'desc'};
648 Usage : $id = $seq->id()
649 Function: This is mapped on display_id
658 return shift->display_id(@_);
664 Usage : if( $obj->is_circular) { /Do Something/ }
665 Function: Returns true if the molecule is circular
666 Returns : Boolean value
674 return $self->{'is_circular'} = shift if @_;
675 return $self->{'is_circular'};
679 =head1 Methods for Bio::IdentifiableI compliance
686 Usage : $string = $obj->object_id()
687 Function: A string which represents the stable primary identifier
688 in this namespace of this object. For DNA sequences this
689 is its accession_number, similarly for protein sequences.
691 This is aliased to accession_number().
698 return shift->accession_number(@_);
704 Usage : $version = $obj->version()
705 Function: A number which differentiates between versions of
706 the same object. Higher numbers are considered to be
707 later and more relevant, but a single object described
708 the same identifier should represent the same concept.
715 my ($self,$value) = @_;
716 if( defined $value) {
717 $self->{'_version'} = $value;
719 return $self->{'_version'};
726 Usage : $authority = $obj->authority()
727 Function: A string which represents the organisation which
728 granted the namespace, written as the DNS name for
729 organisation (eg, wormbase.org).
736 my ($obj,$value) = @_;
737 if( defined $value) {
738 $obj->{'authority'} = $value;
740 return $obj->{'authority'};
746 Usage : $string = $obj->namespace()
747 Function: A string representing the name space this identifier
748 is valid in, often the database name or the name
749 describing the collection.
757 my ($self,$value) = @_;
758 if( defined $value) {
759 $self->{'namespace'} = $value;
761 return $self->{'namespace'} || "";
764 =head1 Methods for Bio::DescribableI compliance
766 This comprises of display_name and description.
773 Usage : $string = $obj->display_name()
774 Function: A string which is what should be displayed to the user.
775 The string should have no spaces (ideally, though a cautious
776 user of this interface would not assumme this) and should be
777 less than thirty characters (though again, double checking
778 this is a good idea).
780 This is aliased to display_id().
786 return shift->display_id(@_);
792 Usage : $string = $obj->description()
793 Function: A text string suitable for displaying to the user a
794 description. This string is likely to have spaces, but
795 should not have any newlines or formatting - just plain
796 text. The string should not be greater than 255 characters
797 and clients can feel justified at truncating strings at 255
798 characters for the purposes of display.
800 This is aliased to desc().
806 return shift->desc(@_);
809 =head1 Methods Inherited from Bio::PrimarySeqI
811 These methods are available on Bio::PrimarySeq, although they are
812 actually implemented on Bio::PrimarySeqI
817 Usage : $rev = $seq->revcom()
818 Function: Produces a new Bio::SeqI implementing object which
819 is the reversed complement of the sequence. For protein
820 sequences this throws an exception of
821 "Sequence is a protein. Cannot revcom".
823 The id is the same id as the orginal sequence, and the
824 accession number is also indentical. If someone wants to
825 track that this sequence has be reversed, it needs to
826 define its own extensions.
828 To do an inplace edit of an object you can go:
830 $seqobj = $seqobj->revcom();
832 This of course, causes Perl to handle the garbage
833 collection of the old object, but it is roughly speaking as
834 efficient as an inplace edit.
836 Returns : A new (fresh) Bio::SeqI object
844 Usage : $subseq = $myseq->trunc(10,100);
845 Function: Provides a truncation of a sequence,
848 Returns : A fresh Bio::SeqI implementing object.
854 =head1 Internal methods
856 These are internal methods to PrimarySeq
860 =head2 _guess_alphabet
862 Title : _guess_alphabet
864 Function: Automatically guess and set the type of sequence: dna, rna or protein
866 Returns : one of strings 'dna', 'rna' or 'protein'.
872 sub _guess_alphabet {
876 # Remove char's that clearly don't denote nucleic or amino acids
877 my $str = $self->seq();
880 # Check for sequences without valid letters
881 my $total = CORE::length($str);
883 if (!$self->{'_nowarnonempty'}) {
884 $self->warn("Got a sequence with no letters in it ".
885 "cannot guess alphabet");
890 if ($str =~ m/[EFIJLOPQXZ]/i) {
891 # Start with a safe method to find proteins.
892 # Unambiguous IUPAC letters for proteins are: E,F,I,J,L,O,P,Q,X,Z
895 # Alphabet is unsure, could still be DNA, RNA or protein.
896 # DNA and RNA contain mostly A, T, U, G, C and N, but the other letters
897 # they use are also among the 15 valid letters that a protein sequence
898 # can contain at this stage. Make our best guess based on sequence
899 # composition. If it contains over 70% of ACGTUN, it is likely nucleic.
900 if( ($str =~ tr/ATUGCNatugcn//) / $total > 0.7 ) {
901 if ( $str =~ m/U/i ) {
911 $self->alphabet($type);
915 ############################################################################
916 # aliases due to name changes or to compensate for our lack of consistency #
917 ############################################################################
922 $self->warn(ref($self)."::accession is deprecated, ".
923 "use accession_number() instead");
924 return $self->accession_number(@_);