3 # bioperl module for Bio::PrimarySeq
5 # Cared for by Ewan Birney <birney@ebi.ac.uk>
7 # Copyright Ewan Birney
9 # You may distribute this module under the same terms as perl itself
11 # POD documentation - main docs before the code
15 Bio::PrimarySeq - Bioperl lightweight Sequence Object
19 # Bio::SeqIO for file reading, Bio::DB::GenBank for
28 $seqobj = Bio::PrimarySeq->new ( -seq => 'ATGGGGTGGGCGGTGGGTGGTTTG',
29 -id => 'GeneFragment-12',
30 -accession_number => 'X78121',
33 print "Sequence ", $seqobj->id(), " with accession ",
34 $seqobj->accession_number, "\n";
38 $inputstream = Bio::SeqIO->new(-file => "myseq.fa",
40 $seqobj = $inputstream->next_seq();
41 print "Sequence ", $seqobj->id(), " and desc ", $seqobj->desc, "\n";
43 # to get out parts of the sequence.
45 print "Sequence ", $seqobj->id(), " with accession ",
46 $seqobj->accession_number, " and desc ", $seqobj->desc, "\n";
48 $string = $seqobj->seq();
49 $string2 = $seqobj->subseq(1,40);
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
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).
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
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.
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.
86 bioperl-l@bioperl.org - General discussion
87 http://bioperl.org/wiki/Mailing_lists - About the mailing lists
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
95 http://bugzilla.open-bio.org/
97 =head1 AUTHOR - Ewan Birney
99 Email birney@ebi.ac.uk
103 The rest of the documentation details each of the object
104 methods. Internal methods are usually preceded with a _
109 # Let the code begin...
112 package Bio
::PrimarySeq
;
113 use vars
qw($MATCHPATTERN $GAP_SYMBOLS);
116 $MATCHPATTERN = 'A-Za-z\-\.\*\?=~';
119 use base qw(Bio::Root::Root Bio::PrimarySeqI
120 Bio::IdentifiableI Bio::DescribableI);
123 # setup the allowed values for alphabet()
126 my %valid_type = map {$_, 1} qw( dna rna protein );
131 Usage : $seq = Bio::PrimarySeq->new( -seq => 'ATGGGGGTGGTGGTACCCT',
133 -accession_number => 'AL000012',
136 Function: Returns a new primary seq object from
137 basic constructors, being a string for the sequence
138 and strings for id and accession_number.
140 Note that you can provide an empty sequence string. However, in
141 this case you MUST specify the type of sequence you wish to
142 initialize by the parameter -alphabet. See alphabet() for possible
144 Returns : a new Bio::PrimarySeq object
145 Args : -seq => sequence string
146 -display_id => display id of the sequence (locus name)
147 -accession_number => accession number
148 -primary_id => primary id (Genbank id)
149 -version => version number
150 -namespace => the namespace for the accession
151 -authority => the authority for the namespace
152 -description => description text
153 -desc => alias for description
154 -alphabet => sequence type (alphabet) (dna|rna|protein)
155 -id => alias for display id
156 -is_circular => boolean field for whether or not sequence is circular
157 -direct => boolean field for directly setting sequence (requires alphabet also set)
158 -ref_to_seq => boolean field indicating the sequence is a reference (?!?)
159 -nowarnonempty => boolean field for whether or not to warn when sequence is empty
165 my ($class, @args) = @_;
166 my $self = $class->SUPER::new
(@args);
168 my($seq,$id,$acc,$pid,$ns,$auth,$v,$oid,
170 $alphabet,$given_id,$is_circular,$direct,$ref_to_seq,$len,$nowarnonempty) =
171 $self->_rearrange([qw(SEQ
191 # private var _nowarnonempty, need to be set before calling _guess_alphabet
192 $self->{'_nowarnonempty'} = $nowarnonempty;
194 if( defined $id && defined $given_id ) {
195 if( $id ne $given_id ) {
196 $self->throw("Provided both id and display_id constructor ".
197 "functions. [$id] [$given_id]");
200 if( defined $given_id ) { $id = $given_id; }
202 # let's set the length before the seq -- if there is one, this length is
203 # going to be invalidated
204 defined $len && $self->length($len);
206 # if alphabet is provided we set it first, so that it won't be guessed
207 # when the sequence is set
208 $alphabet && $self->alphabet($alphabet);
210 # if there is an alphabet, and direct is passed in, assume the alphabet
213 if( $direct && $ref_to_seq) {
214 $self->{'seq'} = $$ref_to_seq;
216 $self->_guess_alphabet();
217 } # else it has been set already above
219 # print STDERR "DEBUG: setting sequence to [$seq]\n";
220 # note: the sequence string may be empty
221 $self->seq($seq) if defined($seq);
224 $id && $self->display_id($id);
225 $acc && $self->accession_number($acc);
226 defined $pid && $self->primary_id($pid);
227 $desc && $self->desc($desc);
228 $description && $self->description($description);
229 $is_circular && $self->is_circular($is_circular);
230 $ns && $self->namespace($ns);
231 $auth && $self->authority($auth);
232 defined($v) && $self->version($v);
233 defined($oid) && $self->object_id($oid);
241 return $obj->{'seq'} = shift if @_;
249 Usage : $string = $obj->seq()
250 Function: Returns the sequence as a string of letters. The
251 case of the letters is left up to the implementer.
252 Suggested cases are upper case for proteins and lower case for
253 DNA sequence (IUPAC standard), but you should not rely on this.
255 Args : Optionally on set the new value (a string). An optional second
256 argument presets the alphabet (otherwise it will be guessed).
261 my ($obj,@args) = @_;
263 if( scalar(@args) == 0 ) {
264 return $obj->{'seq'};
267 my ($value,$alphabet) = @args;
270 if(defined($value) && (! $obj->validate_seq($value))) {
271 $obj->throw("Attempting to set the sequence to [$value] ".
272 "which does not look healthy");
274 # if a sequence was already set we make sure that we re-adjust the
275 # alphabet, otherwise we skip guessing if alphabet is already set
276 # note: if the new seq is empty or undef, we don't consider that a
277 # change (we wouldn't have anything to guess on anyway)
279 exists($obj->{'seq'}) && (CORE
::length($value || '') > 0);
280 $obj->{'seq'} = $value;
281 # new alphabet overridden by arguments?
283 # yes, set it no matter what
284 $obj->alphabet($alphabet);
285 } elsif( # if we changed a previous sequence to a new one
287 # or if there is no alphabet yet at all
288 (! defined($obj->alphabet()))) {
289 # we need to guess the (possibly new) alphabet
290 $obj->_guess_alphabet();
291 } # else (seq not changed and alphabet was defined) do nothing
292 # if the seq is changed, make sure we unset a possibly set length
293 $obj->length(undef) if $is_changed_seq || $obj->{'seq'};
295 return $obj->{'seq'};
301 Usage : if(! $seq->validate_seq($seq_str) ) {
302 print "sequence $seq_str is not valid for an object of
303 alphabet ",$seq->alphabet, "\n";
305 Function: Validates a given sequence string. A validating sequence string
306 must be accepted by seq(). A string that does not validate will
307 lead to an exception if passed to seq().
309 The implementation provided here does not take alphabet() into
310 account. Allowed are all letters (A-Z) and '-','.','*','?','=',
314 Returns : 1 if the supplied sequence string is valid for the object, and
316 Args : The sequence string to be validated.
322 my ($self,$seqstr) = @_;
323 if( ! defined $seqstr ){ $seqstr = $self->seq(); }
324 return 0 unless( defined $seqstr);
325 if((CORE
::length($seqstr) > 0) &&
326 ($seqstr !~ /^([$MATCHPATTERN]+)$/)) {
327 $self->warn("seq doesn't validate, mismatch is " .
328 join(",",($seqstr =~ /([^$MATCHPATTERN]+)/g)));
337 Usage : $substring = $obj->subseq(10,40);
338 $substring = $obj->subseq(10,40,NOGAP)
339 $substring = $obj->subseq(-START=>10,-END=>40,-REPLACE_WITH=>'tga')
340 Function: returns the subseq from start to end, where the first sequence
341 character has coordinate 1 number is inclusive, ie 1-2 are the
342 first two characters of the sequence
344 Args : integer for start position
345 integer for end position
347 Bio::LocationI location for subseq (strand honored)
348 Specify -NOGAP=>1 to return subseq with gap characters removed
349 Specify -REPLACE_WITH=>$new_subseq to replace the subseq returned
350 with $new_subseq in the sequence object
357 my ($start,$end,$nogap,$replace) = $self->_rearrange([qw(START
360 REPLACE_WITH)],@args);
362 # if $replace is specified, have the constructor validate it as seq
363 my $dummy = new Bio
::PrimarySeq
(-seq
=>$replace, -alphabet
=>$self->alphabet) if defined($replace);
365 if( ref($start) && $start->isa('Bio::LocationI') ) {
368 foreach my $subloc ($loc->each_Location()) {
369 my $piece = $self->subseq(-START
=>$subloc->start(),
370 '-END'=>$subloc->end(),
371 -REPLACE_WITH
=>$replace,
373 $piece =~ s/[$GAP_SYMBOLS]//g if $nogap;
374 if($subloc->strand() < 0) {
375 $piece = Bio
::PrimarySeq
->new('-seq' => $piece)->revcom()->seq();
380 } elsif( defined $start && defined $end ) {
382 $self->throw("Bad start,end parameters. Start [$start] has to be ".
383 "less than end [$end]");
386 $self->throw("Bad start parameter ($start). Start must be positive.");
388 if( $end > $self->length ) {
389 $self->throw("Bad end parameter ($end). End must be less than the total length of sequence (total=".$self->length.")");
392 # remove one from start, and then length is end-start
394 my @ss_args = map { eval "defined $_" ?
$_ : () } qw( $self->{seq} $start $end-$start $replace);
395 my $seqstr = eval join( '', "substr(", join(',',@ss_args), ")");
396 $seqstr =~ s/[$GAP_SYMBOLS]//g if ($nogap);
400 $self->warn("Incorrect parameters to subseq - must be two integers or a Bio::LocationI object. Got:", $self,$start,$end,$replace,$nogap);
408 Usage : $len = $seq->length();
409 Function: Get the length of the sequence in number of symbols (bases
412 You can also set this attribute, even to a number that does
413 not match the length of the sequence string. This is useful
414 if you don''t want to set the sequence too, or if you want
415 to free up memory by unsetting the sequence. In the latter
416 case you could do e.g.
418 $seq->length($seq->length);
421 Note that if you set the sequence to a value other than
422 undef at any time, the length attribute will be
423 invalidated, and the length of the sequence string will be
424 reported again. Also, we won''t let you lie about the length.
427 Returns : integer representing the length of the sequence.
428 Args : Optionally, the value on set
434 my $len = CORE::length($self->seq() || '');
438 if(defined($val) && $len && ($len != $val)) {
439 $self->throw("You're trying to lie about the length: ".
440 "is $len but you say ".$val);
442 $self->{'_seq_length'} = $val;
443 } elsif(defined($self->{'_seq_length'})) {
444 return $self->{'_seq_length'};
451 Title : display_id or display_name
452 Usage : $id_string = $obj->display_id();
453 Function: returns the display id, aka the common name of the Sequence object.
455 The semantics of this is that it is the most likely string to
456 be used as an identifier of the sequence, and likely to have
457 "human" readability. The id is equivalent to the ID field of
458 the GenBank/EMBL databanks and the id field of the
459 Swissprot/sptrembl database. In fasta format, the >(\S+) is
460 presumed to be the id, though some people overload the id to
461 embed other information. Bioperl does not use any embedded
462 information in the ID field, and people are encouraged to use
463 other mechanisms (accession field for example, or extending
464 the sequence object) to solve this.
466 With the new Bio::DescribeableI interface, display_name aliases
476 my ($obj,$value) = @_;
477 if( defined $value) {
478 $obj->{'display_id'} = $value;
480 return $obj->{'display_id'};
483 =head2 accession_number
485 Title : accession_number or object_id
486 Usage : $unique_key = $obj->accession_number;
487 Function: Returns the unique biological id for a sequence, commonly
488 called the accession_number. For sequences from established
489 databases, the implementors should try to use the correct
490 accession number. Notice that primary_id() provides the
491 unique id for the implemetation, allowing multiple objects
492 to have the same accession number in a particular implementation.
494 For sequences with no accession number, this method should
497 [Note this method name is likely to change in 1.3]
499 With the new Bio::IdentifiableI interface, this is aliased
503 Args : A string (optional) for setting
507 sub accession_number {
508 my( $obj, $acc ) = @_;
511 $obj->{'accession_number'} = $acc;
513 $acc = $obj->{'accession_number'};
514 $acc = 'unknown' unless defined $acc;
522 Usage : $unique_key = $obj->primary_id;
523 Function: Returns the unique id for this object in this
524 implementation. This allows implementations to manage their
525 own object ids in a way the implementaiton can control
526 clients can expect one id to map to one object.
528 For sequences with no natural primary id, this method
529 should return a stringified memory location.
532 Args : A string (optional, for setting)
540 $obj->{'primary_id'} = shift;
542 if( ! defined($obj->{'primary_id'}) ) {
545 return $obj->{'primary_id'};
552 Usage : if( $obj->alphabet eq 'dna' ) { /Do Something/ }
553 Function: Get/Set the alphabet of sequence, one of
554 'dna', 'rna' or 'protein'. This is case sensitive.
556 This is not called <type> because this would cause
557 upgrade problems from the 0.5 and earlier Seq objects.
559 Returns : a string either 'dna','rna','protein'. NB - the object must
560 make a call of the type - if there is no alphabet specified it
562 Args : optional string to set : 'dna' | 'rna' | 'protein'
568 my ($obj,$value) = @_;
569 if (defined $value) {
571 unless ( $valid_type{$value} ) {
572 $obj->throw("Alphabet '$value' is not a valid alphabet (".
573 join(',', map "'$_'", sort keys %valid_type) .
576 $obj->{'alphabet'} = $value;
578 return $obj->{'alphabet'};
583 Title : desc or description
584 Usage : $obj->desc($newval)
585 Function: Get/set description of the sequence.
587 'description' is an alias for this for compliance with the
588 Bio::DescribeableI interface.
591 Returns : value of desc (a string)
592 Args : newvalue (a string or undef, optional)
600 return $self->{'desc'} = shift if @_;
601 return $self->{'desc'};
625 Usage : $id = $seq->id()
626 Function: This is mapped on display_id
635 return shift->display_id(@_);
641 Usage : if( $obj->is_circular) { /Do Something/ }
642 Function: Returns true if the molecule is circular
643 Returns : Boolean value
651 return $self->{'is_circular'} = shift if @_;
652 return $self->{'is_circular'};
656 =head1 Methods for Bio::IdentifiableI compliance
663 Usage : $string = $obj->object_id()
664 Function: A string which represents the stable primary identifier
665 in this namespace of this object. For DNA sequences this
666 is its accession_number, similarly for protein sequences.
668 This is aliased to accession_number().
675 return shift->accession_number(@_);
681 Usage : $version = $obj->version()
682 Function: A number which differentiates between versions of
683 the same object. Higher numbers are considered to be
684 later and more relevant, but a single object described
685 the same identifier should represent the same concept.
692 my ($self,$value) = @_;
693 if( defined $value) {
694 $self->{'_version'} = $value;
696 return $self->{'_version'};
703 Usage : $authority = $obj->authority()
704 Function: A string which represents the organisation which
705 granted the namespace, written as the DNS name for
706 organisation (eg, wormbase.org).
713 my ($obj,$value) = @_;
714 if( defined $value) {
715 $obj->{'authority'} = $value;
717 return $obj->{'authority'};
723 Usage : $string = $obj->namespace()
724 Function: A string representing the name space this identifier
725 is valid in, often the database name or the name
726 describing the collection.
734 my ($self,$value) = @_;
735 if( defined $value) {
736 $self->{'namespace'} = $value;
738 return $self->{'namespace'} || "";
741 =head1 Methods for Bio::DescribableI compliance
743 This comprises of display_name and description.
750 Usage : $string = $obj->display_name()
751 Function: A string which is what should be displayed to the user.
752 The string should have no spaces (ideally, though a cautious
753 user of this interface would not assumme this) and should be
754 less than thirty characters (though again, double checking
755 this is a good idea).
757 This is aliased to display_id().
763 return shift->display_id(@_);
769 Usage : $string = $obj->description()
770 Function: A text string suitable for displaying to the user a
771 description. This string is likely to have spaces, but
772 should not have any newlines or formatting - just plain
773 text. The string should not be greater than 255 characters
774 and clients can feel justified at truncating strings at 255
775 characters for the purposes of display.
777 This is aliased to desc().
783 return shift->desc(@_);
786 =head1 Methods Inherited from Bio::PrimarySeqI
788 These methods are available on Bio::PrimarySeq, although they are
789 actually implemented on Bio::PrimarySeqI
794 Usage : $rev = $seq->revcom()
795 Function: Produces a new Bio::SeqI implementing object which
796 is the reversed complement of the sequence. For protein
797 sequences this throws an exception of
798 "Sequence is a protein. Cannot revcom".
800 The id is the same id as the orginal sequence, and the
801 accession number is also indentical. If someone wants to
802 track that this sequence has be reversed, it needs to
803 define its own extensions.
805 To do an inplace edit of an object you can go:
807 $seqobj = $seqobj->revcom();
809 This of course, causes Perl to handle the garbage
810 collection of the old object, but it is roughly speaking as
811 efficient as an inplace edit.
813 Returns : A new (fresh) Bio::SeqI object
821 Usage : $subseq = $myseq->trunc(10,100);
822 Function: Provides a truncation of a sequence,
825 Returns : A fresh Bio::SeqI implementing object.
831 =head1 Internal methods
833 These are internal methods to PrimarySeq
837 =head2 _guess_alphabet
839 Title : _guess_alphabet
841 Function: Determines (and sets) the type of sequence: dna, rna, protein
843 Returns : one of strings 'dna', 'rna' or 'protein'.
849 sub _guess_alphabet {
853 #return if $self->alphabet;
855 my $str = $self->seq();
856 # Remove char's that clearly denote ambiguity
859 my $total = CORE::length($str);
861 if (!$self->{'_nowarnonempty'}) {
862 $self->warn("Got a sequence with no letters in it ".
863 "cannot guess alphabet");
868 my $u = ($str =~ tr/Uu//);
869 # The assumption here is that most of sequences comprised of mainly
870 # ATGC, with some N, will be 'dna' despite the fact that N could
872 my $atgc = ($str =~ tr/ATGCNatgcn//);
874 if( ($atgc / $total) > 0.85 ) {
876 } elsif( (($atgc + $u) / $total) > 0.85 ) {
882 $self->alphabet($type);
886 ############################################################################
887 # aliases due to name changes or to compensate for our lack of consistency #
888 ############################################################################
893 $self->warn(ref($self)."::accession is deprecated, ".
894 "use accession_number() instead");
895 return $self->accession_number(@_);