# The Bio::SeqIO for file reading, Bio::DB::GenBank for
  # database reading

  use Bio::Seq;
  use Bio::SeqIO;
  use Bio::DB::GenBank;

  #make from memory
  $seqobj = Bio::PrimarySeq->new ( -seq => 'ATGGGGTGGGCGGTGGGTGGTTTG',
				   -id  => 'GeneFragment-12',
				   -accession_number => 'X78121',
				   -alphabet => 'dna',
				   -is_circular => 1
				   );
  print "Sequence ", $seqobj->id(), " with accession ", 
    $seqobj->accession_number, "\n";

  # read from file
  $inputstream = Bio::SeqIO->new(-file => "myseq.fa",-format => 'Fasta');
  $seqobj = $inputstream->next_seq();
  print "Sequence ", $seqobj->id(), " and desc ", $seqobj->desc, "\n";

  # to get out parts of the sequence.

  print "Sequence ", $seqobj->id(), " with accession ", 
    $seqobj->accession_number, " and desc ", $seqobj->desc, "\n";

  $string  = $seqobj->seq();
  $string2 = $seqobj->subseq(1,40);

   my ($self) = @_;
   my ($str,$str2,$total,$atgc,$u,$type);

   $str = $self->seq();
   $str =~ s/\-\.\?//g;

   $total = CORE::length($str);
   if( $total == 0 ) {
       $self->throw("Got a sequence with no letters in - ".
		    "cannot guess alphabet [$str]");
   }
   
   $u = ($str =~ tr/Uu//);
   $atgc = ($str =~ tr/ATGCNatgcn//);
   
   if( ($atgc / $total) > 0.85 ) {
       $type = 'dna';
   } elsif( (($atgc + $u) / $total) > 0.85 ) {
       $type = 'rna';
   } else {
       $type = 'protein';
   }

   $self->alphabet($type);
   return $type;
}

############################################################################
# aliases due to name changes or to compensate for our lack of consistency #
############################################################################

    my $self = shift;

    $self->warn(ref($self)."::accession is deprecated, ".
		"use accession_number() instead");
    return $self->accession_number(@_);
}

1;

    my( $obj, $acc ) = @_;

    if (defined $acc) {
        $obj->{'accession_number'} = $acc;
    } else {
        $acc = $obj->{'accession_number'};
        $acc = 'unknown' unless defined $acc;
    }
    return $acc;

    my ($obj,$value) = @_;
    if (defined $value) {
	$value = lc $value;
	unless ( $valid_type{$value} ) {
	    $obj->throw("Molecular type '$value' is not a valid type (".
			join(',', map "'$_'", sort keys %valid_type) .
			") lowercase");
	}
	$obj->{'alphabet'} = $value;
    }
    return $obj->{'alphabet'};

    my ($obj,$value) = @_;
    if( defined $value) {
	$obj->{'authority'} = $value;
    }
    return $obj->{'authority'};

   my ($self) = @_;

   return 1;

    my $self = shift;

    return $self->{'desc'} = shift if @_;
    return $self->{'desc'};

    return shift->desc(@_);

    my $obj = shift;
    return $obj->{'seq'} = shift if @_;
    return undef;

   my ($obj,$value) = @_;
   if( defined $value) {
      $obj->{'display_id'} = $value;
    }
    return $obj->{'display_id'};

    return shift->display_id(@_);

   return shift->display_id(@_);

    my $self = shift;
    return $self->{'is_circular'} = shift if @_;
    return $self->{'is_circular'};

    my $self = shift;
    my $len = CORE::length($self->seq() || '');
    
    if(@_) {
	my $val = shift;
	if(defined($val) && $len && ($len != $val)) {
	    $self->throw("You're trying to lie about the length: ".
			 "is $len but you say ".$val);
	}
	$self->{'_seq_length'} = $val;
    } elsif(defined($self->{'_seq_length'})) {
	return $self->{'_seq_length'};
    }
    return $len;

    my ($self,$value) = @_;
    if( defined $value) {
	$self->{'namespace'} = $value;
    }
    return $self->{'namespace'} || "";

    my ($class, @args) = @_;
    my $self = $class->SUPER::new(@args);

    my($seq,$id,$acc,$pid,$ns,$auth,$v,$oid,
       $desc,$alphabet,$given_id,$is_circular,$direct,$ref_to_seq,$len) =
	$self->_rearrange([qw(SEQ
			      DISPLAY_ID
			      ACCESSION_NUMBER
			      PRIMARY_ID
			      NAMESPACE
			      AUTHORITY
			      VERSION
			      OBJECT_ID
			      DESC
			      ALPHABET
			      ID
			      IS_CIRCULAR
			      DIRECT
			      REF_TO_SEQ
			      LENGTH
			      )],
			  @args);
    if( defined $id && defined $given_id ) {
	if( $id ne $given_id ) {
	    $self->throw("Provided both id and display_id constructor ".
			 "functions. [$id] [$given_id]");	
	}
    }
    if( defined $given_id ) { $id = $given_id; }

    # let's set the length before the seq -- if there is one, this length is
    # going to be invalidated
    defined $len && $self->length($len);

    # if alphabet is provided we set it first, so that it won't be guessed
    # when the sequence is set
    $alphabet && $self->alphabet($alphabet);
    
    # if there is an alphabet, and direct is passed in, assumme the alphabet
    # and sequence is ok 

    if( $direct && $ref_to_seq) {
	$self->{'seq'} = $$ref_to_seq;
	if( ! $alphabet ) {
	    $self->_guess_alphabet();
	} # else it has been set already above
    } else {
#	print STDERR "DEBUG: setting sequence to [$seq]\n";
	# note: the sequence string may be empty
	$self->seq($seq) if defined($seq);
    }

    $id          && $self->display_id($id);
    $acc         && $self->accession_number($acc);
    defined $pid && $self->primary_id($pid);
    $desc        && $self->desc($desc);
    $is_circular && $self->is_circular($is_circular);
    $ns          && $self->namespace($ns);
    $auth        && $self->authority($auth);
    defined($v)  && $self->version($v);
    defined($oid) && $self->object_id($oid);

    return $self;

    return shift->accession_number(@_);

   my ($obj,$value) = @_;
   if( defined $value) {
      $obj->{'primary_id'} = $value;
    }
   if( ! exists $obj->{'primary_id'} ) {
       return "$obj";
   }
   return $obj->{'primary_id'};

   my ($obj,@args) = @_;

   if( scalar(@args) == 0 ) {
       return $obj->{'seq'};
   }

   my ($value,$alphabet) = @args;


   if(@args) {
       if(defined($value) && (! $obj->validate_seq($value))) {
	   $obj->throw("Attempting to set the sequence to [$value] ".
		       "which does not look healthy");
       }
       # if a sequence was already set we make sure that we re-adjust the
       # mol.type, otherwise we skip guessing if mol.type is already set
       # note: if the new seq is empty or undef, we don't consider that a
       # change (we wouldn't have anything to guess on anyway)
       my $is_changed_seq =
	   exists($obj->{'seq'}) && (CORE::length($value || '') > 0);
       $obj->{'seq'} = $value;
       # new alphabet overridden by arguments?
       if($alphabet) {
	   # yes, set it no matter what
	   $obj->alphabet($alphabet);
       } elsif( # if we changed a previous sequence to a new one
		$is_changed_seq ||
		# or if there is no alphabet yet at all
		(! defined($obj->alphabet()))) {
	   # we need to guess the (possibly new) alphabet
	   $obj->_guess_alphabet();
       } # else (seq not changed and alphabet was defined) do nothing
       # if the seq is changed, make sure we unset a possibly set length
       $obj->length(undef) if $is_changed_seq;
   }
   return $obj->{'seq'};

   my ($self,$start,$end,$replace) = @_;

   if( ref($start) && $start->isa('Bio::LocationI') ) {
       my $loc = $start;
       $replace = $end; # do we really use this anywhere? scary. HL
       my $seq = "";
       foreach my $subloc ($loc->each_Location()) {
	   my $piece = $self->subseq($subloc->start(),
				     $subloc->end(), $replace);
	   if($subloc->strand() < 0) {
	       $piece = Bio::PrimarySeq->new('-seq' => $piece)->revcom()->seq();
	   }
	   $seq .= $piece;
       }
       return $seq;
   } elsif(  defined  $start && defined $end ) {
       if( $start > $end ){
	   $self->throw("in subseq, start [$start] has to be ".
			"greater than end [$end]");
       }
       if( $start <= 0 || $end > $self->length ) {
	   $self->throw("You have to have start positive\n\tand length less ".
			"than the total length of sequence [$start:$end] ".
			"Total ".$self->length."");
       }

       # remove one from start, and then length is end-start
       $start--;
       if( defined $replace ) {
	   return substr( $self->seq(), $start, ($end-$start), $replace);
       } else {
	   return substr( $self->seq(), $start, ($end-$start));
       }
   } else {
       $self->warn("Incorrect parameters to subseq - must be two integers ".
		   "or a Bio::LocationI object not ($start,$end)");
   }

    my ($self,$seqstr) = @_;
    if( ! defined $seqstr ){ $seqstr = $self->seq(); }
    return 0 unless( defined $seqstr); 
    if((CORE::length($seqstr) > 0) && ($seqstr !~ /^([A-Za-z\-\.\*\?]+)$/)) {
	$self->warn("seq doesn't validate, mismatch is " .
		   ($seqstr =~ /([^A-Za-z\-\.\*\?]+)/g));
	return 0;
    }
    return 1;

    my ($self,$value) = @_;
    if( defined $value) {
	$self->{'_version'} = $value;
    }
    return $self->{'_version'};

  bioperl-l@bioperl.org             - General discussion
  http://bio.perl.org/MailList.html - About the mailing lists

  bioperl-bugs@bio.perl.org
  http://bugzilla.bioperl.org/

 Title   : revcom
 Usage   : $rev = $seq->revcom()
 Function: Produces a new Bio::SeqI implementing object which
           is the reversed complement of the sequence. For protein
           sequences this throws an exception of
           "Sequence is a protein. Cannot revcom"

           The id is the same id as the orginal sequence, and the
           accession number is also indentical. If someone wants to
           track that this sequence has be reversed, it needs to
           define its own extensions

           To do an inplace edit of an object you can go:

           $seqobj = $seqobj->revcom();

           This of course, causes Perl to handle the garbage
           collection of the old object, but it is roughly speaking as
           efficient as an inplace edit.

 Returns : A new (fresh) Bio::SeqI object
 Args    : none

 Title   : trunc
 Usage   : $subseq = $myseq->trunc(10,100);
 Function: Provides a truncation of a sequence,

 Example :
 Returns : a fresh Bio::SeqI implementing object
 Args    :