3 # bioperl module for Bio::PrimarySeq
5 # Please direct questions and support issues to <bioperl-l@bioperl.org>
7 # Cared for by Ewan Birney <birney@ebi.ac.uk>
9 # Copyright Ewan Birney
11 # You may distribute this module under the same terms as perl itself
13 # POD documentation - main docs before the code
17 Bio::PrimarySeq - Bioperl lightweight Sequence Object
21 # Bio::SeqIO for file reading, Bio::DB::GenBank for
30 $seqobj = Bio::PrimarySeq->new ( -seq => 'ATGGGGTGGGCGGTGGGTGGTTTG',
31 -id => 'GeneFragment-12',
32 -accession_number => 'X78121',
35 print "Sequence ", $seqobj->id(), " with accession ",
36 $seqobj->accession_number, "\n";
40 $inputstream = Bio::SeqIO->new(-file => "myseq.fa",
42 $seqobj = $inputstream->next_seq();
43 print "Sequence ", $seqobj->id(), " and desc ", $seqobj->desc, "\n";
45 # to get out parts of the sequence.
47 print "Sequence ", $seqobj->id(), " with accession ",
48 $seqobj->accession_number, " and desc ", $seqobj->desc, "\n";
50 $string = $seqobj->seq();
51 $string2 = $seqobj->subseq(1,40);
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
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).
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
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.
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.
88 bioperl-l@bioperl.org - General discussion
89 http://bioperl.org/wiki/Mailing_lists - About the mailing lists
93 Please direct usage questions or support issues to the mailing list:
95 L<bioperl-l@bioperl.org>
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.
102 =head2 Reporting Bugs
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
108 http://bugzilla.open-bio.org/
110 =head1 AUTHOR - Ewan Birney
112 Email birney@ebi.ac.uk
116 The rest of the documentation details each of the object
117 methods. Internal methods are usually preceded with a _
122 # Let the code begin...
125 package Bio
::PrimarySeq
;
126 use vars
qw($MATCHPATTERN $GAP_SYMBOLS);
129 $MATCHPATTERN = 'A-Za-z\-\.\*\?=~';
132 use base qw(Bio::Root::Root Bio::PrimarySeqI
133 Bio::IdentifiableI Bio::DescribableI);
136 # setup the allowed values for alphabet()
139 my %valid_type = map {$_, 1} qw( dna rna protein );
144 Usage : $seq = Bio::PrimarySeq->new( -seq => 'ATGGGGGTGGTGGTACCCT',
146 -accession_number => 'AL000012',
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.
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
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
178 my ($class, @args) = @_;
179 my $self = $class->SUPER::new
(@args);
181 my($seq,$id,$acc,$pid,$ns,$auth,$v,$oid,
183 $alphabet,$given_id,$is_circular,$direct,$ref_to_seq,$len,$nowarnonempty) =
184 $self->_rearrange([qw(SEQ
204 # private var _nowarnonempty, need to be set before calling _guess_alphabet
205 $self->{'_nowarnonempty'} = $nowarnonempty;
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]");
213 if( defined $given_id ) { $id = $given_id; }
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);
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);
223 # if there is an alphabet, and direct is passed in, assume the alphabet
226 if( $direct && $ref_to_seq) {
227 $self->{'seq'} = $$ref_to_seq;
229 $self->_guess_alphabet();
230 } # else it has been set already above
232 # print STDERR "DEBUG: setting sequence to [$seq]\n";
233 # note: the sequence string may be empty
234 $self->seq($seq) if defined($seq);
237 $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);
254 return $obj->{'seq'} = shift if @_;
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.
268 Args : Optionally on set the new value (a string). An optional second
269 argument presets the alphabet (otherwise it will be guessed).
274 my ($obj,@args) = @_;
276 if( scalar(@args) == 0 ) {
277 return $obj->{'seq'};
280 my ($value,$alphabet) = @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");
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)
292 exists($obj->{'seq'}) && (CORE
::length($value || '') > 0);
293 $obj->{'seq'} = $value;
294 # new alphabet overridden by arguments?
296 # yes, set it no matter what
297 $obj->alphabet($alphabet);
298 } elsif( # if we changed a previous sequence to a new one
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'};
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("seq doesn't validate, mismatch is " .
341 join(",",($seqstr =~ /([^$MATCHPATTERN]+)/g)));
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
357 Args : integer for start position
358 integer for end position
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
370 my ($start,$end,$nogap,$replace) = $self->_rearrange([qw(START
373 REPLACE_WITH)],@args);
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);
378 if( ref($start) && $start->isa('Bio::LocationI') ) {
381 foreach my $subloc ($loc->each_Location()) {
382 my $piece = $self->subseq(-START
=>$subloc->start(),
383 '-END'=>$subloc->end(),
384 -REPLACE_WITH
=>$replace,
386 $piece =~ s/[$GAP_SYMBOLS]//g if $nogap;
387 if($subloc->strand() < 0) {
388 $piece = Bio
::PrimarySeq
->new('-seq' => $piece)->revcom()->seq();
393 } elsif( defined $start && defined $end ) {
395 $self->throw("Bad start,end parameters. Start [$start] has to be ".
396 "less than end [$end]");
399 $self->throw("Bad start parameter ($start). Start must be positive.");
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.")");
405 # remove one from start, and then length is end-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);
413 $self->warn("Incorrect parameters to subseq - must be two integers or a Bio::LocationI object. Got:", $self,$start,$end,$replace,$nogap);
421 Usage : $len = $seq->length();
422 Function: Get the length of the sequence in number of symbols (bases
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.
431 $seq->length($seq->length);
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.
440 Returns : integer representing the length of the sequence.
441 Args : Optionally, the value on set
447 my $len = CORE::length($self->seq() || '');
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);
455 $self->{'_seq_length'} = $val;
456 } elsif(defined($self->{'_seq_length'})) {
457 return $self->{'_seq_length'};
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.
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.
479 With the new Bio::DescribeableI interface, display_name aliases
489 my ($obj,$value) = @_;
490 if( defined $value) {
491 $obj->{'display_id'} = $value;
493 return $obj->{'display_id'};
496 =head2 accession_number
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.
507 For sequences with no accession number, this method should
510 [Note this method name is likely to change in 1.3]
512 With the new Bio::IdentifiableI interface, this is aliased
516 Args : A string (optional) for setting
520 sub accession_number {
521 my( $obj, $acc ) = @_;
524 $obj->{'accession_number'} = $acc;
526 $acc = $obj->{'accession_number'};
527 $acc = 'unknown' unless defined $acc;
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.
541 For sequences with no natural primary id, this method
542 should return a stringified memory location.
545 Args : A string (optional, for setting)
553 $obj->{'primary_id'} = shift;
555 if( ! defined($obj->{'primary_id'}) ) {
558 return $obj->{'primary_id'};
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.
569 This is not called <type> because this would cause
570 upgrade problems from the 0.5 and earlier Seq objects.
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
575 Args : optional string to set : 'dna' | 'rna' | 'protein'
581 my ($obj,$value) = @_;
582 if (defined $value) {
584 unless ( $valid_type{$value} ) {
585 $obj->throw("Alphabet '$value' is not a valid alphabet (".
586 join(',', map "'$_'", sort keys %valid_type) .
589 $obj->{'alphabet'} = $value;
591 return $obj->{'alphabet'};
596 Title : desc or description
597 Usage : $obj->desc($newval)
598 Function: Get/set description of the sequence.
600 'description' is an alias for this for compliance with the
601 Bio::DescribeableI interface.
604 Returns : value of desc (a string)
605 Args : newvalue (a string or undef, optional)
613 return $self->{'desc'} = shift if @_;
614 return $self->{'desc'};
638 Usage : $id = $seq->id()
639 Function: This is mapped on display_id
648 return shift->display_id(@_);
654 Usage : if( $obj->is_circular) { /Do Something/ }
655 Function: Returns true if the molecule is circular
656 Returns : Boolean value
664 return $self->{'is_circular'} = shift if @_;
665 return $self->{'is_circular'};
669 =head1 Methods for Bio::IdentifiableI compliance
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.
681 This is aliased to accession_number().
688 return shift->accession_number(@_);
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.
705 my ($self,$value) = @_;
706 if( defined $value) {
707 $self->{'_version'} = $value;
709 return $self->{'_version'};
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).
726 my ($obj,$value) = @_;
727 if( defined $value) {
728 $obj->{'authority'} = $value;
730 return $obj->{'authority'};
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.
747 my ($self,$value) = @_;
748 if( defined $value) {
749 $self->{'namespace'} = $value;
751 return $self->{'namespace'} || "";
754 =head1 Methods for Bio::DescribableI compliance
756 This comprises of display_name and description.
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).
770 This is aliased to display_id().
776 return shift->display_id(@_);
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.
790 This is aliased to desc().
796 return shift->desc(@_);
799 =head1 Methods Inherited from Bio::PrimarySeqI
801 These methods are available on Bio::PrimarySeq, although they are
802 actually implemented on Bio::PrimarySeqI
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".
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.
818 To do an inplace edit of an object you can go:
820 $seqobj = $seqobj->revcom();
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.
826 Returns : A new (fresh) Bio::SeqI object
834 Usage : $subseq = $myseq->trunc(10,100);
835 Function: Provides a truncation of a sequence,
838 Returns : A fresh Bio::SeqI implementing object.
844 =head1 Internal methods
846 These are internal methods to PrimarySeq
850 =head2 _guess_alphabet
852 Title : _guess_alphabet
854 Function: Determines (and sets) the type of sequence: dna, rna, protein
856 Returns : one of strings 'dna', 'rna' or 'protein'.
862 sub _guess_alphabet {
866 #return if $self->alphabet;
868 my $str = $self->seq();
869 # Remove char's that clearly denote ambiguity
872 my $total = CORE::length($str);
874 if (!$self->{'_nowarnonempty'}) {
875 $self->warn("Got a sequence with no letters in it ".
876 "cannot guess alphabet");
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
885 my $atgc = ($str =~ tr/ATGCNatgcn//);
887 if( ($atgc / $total) > 0.85 ) {
889 } elsif( (($atgc + $u) / $total) > 0.85 ) {
895 $self->alphabet($type);
899 ############################################################################
900 # aliases due to name changes or to compensate for our lack of consistency #
901 ############################################################################
906 $self->warn(ref($self)."::accession is deprecated, ".
907 "use accession_number() instead");
908 return $self->accession_number(@_);