# 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";

  # get from database
  $db = Bio::DB::GenBank->new();
  $seqobj = $db->get_Seq_by_acc('X78121');

  # 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);

    %valid_type = map {$_, 1} qw( dna rna protein );

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

    my($seq,$id,$acc,$pid,$desc,$alphabet,$given_id,$is_circular) =
	$self->_rearrange([qw(SEQ
			      DISPLAY_ID
			      ACCESSION_NUMBER
			      PRIMARY_ID
			      DESC
			      ALPHABET
			      ID
			      IS_CIRCULAR
			      )],
			  @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; }

    # if alphabet is provided we set it first, so that it won't be guessed
    # when the sequence is set
    $alphabet && $self->alphabet($alphabet);
    # note: the sequence string may be empty
    $self->seq($seq) if defined($seq);
    $id      && $self->display_id($id);
    $acc     && $self->accession_number($acc);
    $pid     && $self->primary_id($pid);
    $desc    && $self->desc($desc);
    $is_circular && $self->is_circular($is_circular);

    return $self;

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

   if( defined $value) {
       if(! $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
       my $is_changed_seq = exists($obj->{'seq'});
       $obj->{'seq'} = $value;
       if(($is_changed_seq && (CORE::length($value) > 0)) ||
	  (! defined($obj->alphabet()))) {
	   $obj->_guess_alphabet();
       }
    }
   return $obj->{'seq'};

   my ($self,$seqstr) = @_;

   if((CORE::length($seqstr) > 0) && ($seqstr !~ /^[A-Za-z\-\.\*\?]+$/)) {
       return 0;
   }
   return 1;

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

   if( ref($start) && $start->isa('Bio::LocationI') ) {
       my $loc = $start;
       if( $loc->length == 0 ) { 
	   $self->warn("Expect location lengths to be > 0");
	   return '';
       } elsif( $loc->end < $loc->start ) { 
	   # what about circular seqs
	   $self->warn("Expect location start to come before location end");
       }
       my $seq = '';
       if( $loc->isa('Bio::Location::SplitLocationI') ) {
	   foreach my $subloc ( $loc->sub_Location ) {
	       my $piece = substr( $self->seq(), $subloc->start - 1, 
				   $subloc->length);
	       if( $subloc->strand < 0 ) { 
		   $piece = Bio::PrimarySeq->new('-seq' => 
						 $piece)->revcom()->seq();
	       }
	       $seq .= $piece;
	   }
       } else { 
	   $seq = substr( $self->seq(), $loc->start - 1, $loc->length);
       }
       if( $loc->strand < 0 ) { 
	   $seq = Bio::PrimarySeq->new('-seq' => $seq)->revcom()->seq();
       }
       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 and length less ".
			"than the total length of sequence [$start:$end] ".
			"Total ".$self->length."");
       }

       # remove one from start, and then length is end-start
       $start--;

       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)= @_;
   my $seq = $self->seq();
   return 0 if ( !defined $seq );
   return CORE::length($seq);

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

    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) {
      $obj->{'primary_id'} = $value;
    }
   if( ! exists $obj->{'primary_id'} ) {
       return "$obj";
   }
   return $obj->{'primary_id'};

    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->{'desc'} = $value;
    }
    return $obj->{'desc'};

   my ($self) = @_;

   return 1;

   my ($self,$value) = @_;

   if( defined $value ) {
	return $self->display_id($value);
   }
   return $self->display_id();

   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]");
   }

   $str2 = $str;

   $str2 =~ s/[ATGCNatgcn]//g;
   $atgc = $total - CORE::length $str2;
   $str = $str2;
   $str2 =~ s/[Uu]//g;

   $u = CORE::length($str) - CORE::length($str2);


   if( ($atgc / $total) > 0.85 ) {
       $type = 'dna';
   } elsif( (($atgc + $u) / $total) > 0.85 ) {

    # old style reading from files
    $seq = Bio::Seq->new( -file => "myfile");

    # new style reading from files.
    $inputstream = Bio::SeqIO->new( -file => "myfile", -format => 'Fasta');
    $seqobj = $inputstream->next_seq();

     # old style writing to files
     print OUTPUT $seq->layout_fasta;

     # new style writing to files
     $outputstream = Bio::SeqIO->new( -fh => \*OUTPUT, -format => 'Fasta');
     $outputstream->write_seq($seqobj);

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

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

 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    :