use Bio::Tools::pICalculator;
  use Bio::SeqIO;

  my $in = Bio::SeqIO->new( -fh => \*STDIN ,-format => 'Fasta' );

  my $calc = Bio::Tools::pICalculator->new(-places => 2,
                                           -pKset => 'EMBOSS');

  while ( my $seq = $in->next_seq ) {
     $calc->seq($seq);
     my $iep = $calc->iep;
     print sprintf( "%s\t%s\t%.2f\n",
                    $seq->id,
                    $iep,
                    $calc->charge_at_pH($iep) );

     for( my $i = 0; $i <= 14; $i += 0.5 ){
        print sprintf( "pH = %.2f\tCharge = %.2f\n",
                       $i,
                       $calc->charge_at_pH($i) );
     }
  }

   my( $class, %opts ) = @_;
   my $self = $class->SUPER::new(%opts);
   $self = bless {}, ref $self || $self;
   $self->seq( $opts{-seq} ) if exists $opts{-seq};
   $self->pKset( $opts{-pKset} || 'EMBOSS' );
   exists $opts{-places} ? $self->places( $opts{-places} ) :
     $self->places(2);
   return $self;

   my( $this, $seq ) = @_;
   unless( defined $seq && UNIVERSAL::isa($seq,'Bio::Seq') ){
      die $seq . " is not a valid Bio::Seq object\n";
   }
   $this->{-seq} = $seq;
   $this->{-count} = count_charged_residues( $seq );
   return $this->{-seq};

   my ( $this, $pKset ) = @_;
   if( ref $pKset eq 'HASH' ){         # user defined pK values
      $this->{-pKset} = $pKset;
   }elsif( $pKset =~ /^emboss$/i ){    # from EMBOSS's iep program
      $this->{-pKset} = $Emboss_pK;
   }elsif( $pKset =~ /^dtaselect$/i ){ # from DTASelect (scripps)
      $this->{-pKset} = $DTASelect_pK;
   }else{                              # default to EMBOSS
      $this->{-pKset} = $Emboss_pK;
   }
   return $this->{-pKset};

   my $this = shift;
   $this->{-places} = shift if @_;
   return $this->{-places};

   my $this = shift;
   return _calculate_iep($this->{-pKset},
	    		 $this->{-places},
			 $this->{-seq},
			 $this->{-count}
		        );

   my $this = shift;
   return _calculate_charge_at_pH( shift, $this->{-pKset},
				  $this->{-count} );

   my $seq = shift;
   my $sequence = $seq->seq;
   my $count;
   for ( qw( K R H D E C Y ) ){ # charged AA's
      $count->{$_}++ while $sequence =~ /$_/ig;
   }
   return $count;

    my( $pK, $places, $seq, $count ) = @_;
    my $pH = 7.0;
    my $step = 3.5;
    my $last_charge = 0.0;
    my $format = "%.${places}f";

    unless( defined $count ){
       $count = count_charged_residues($seq);
    }
    while(1){
       my $charge = _calculate_charge_at_pH( $pH, $pK, $count );
       last if sprintf($format,$charge) == 
	 sprintf($format,$last_charge);
       $charge > 0 ? ( $pH += $step ) : ( $pH -= $step );
       $step /= 2.0;
       $last_charge = $charge;
    }
    return sprintf( $format, $pH );

   no warnings; # don't complain if a given key doesn't exist
   my( $pH, $pK, $count ) = @_;
   my $charge =    _partial_charge( $pK->{N_term}, $pH )
   + $count->{K} * _partial_charge( $pK->{K},      $pH )
   + $count->{R} * _partial_charge( $pK->{R},      $pH )
   + $count->{H} * _partial_charge( $pK->{H},      $pH )
   - $count->{D} * _partial_charge( $pH,      $pK->{D} )
   - $count->{E} * _partial_charge( $pH,      $pK->{E} )
   - $count->{C} * _partial_charge( $pH,      $pK->{C} )
   - $count->{Y} * _partial_charge( $pH,      $pK->{Y} )
   -               _partial_charge( $pH, $pK->{C_term} );
   return $charge;

   my $cr = 10 ** ( $_[0] - $_[1] );
   return $cr / ( $cr + 1 );

  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   : new
 Usage   : Bio::Tools::pICalculator->new
 Function: Instantiates the Bio::Tools::pICalculator object
 Example : $calc = Bio::Tools::pICalculator->new( -pKset => \%pKvalues,
                                                  # a Bio::Seq object
                                                  -seq => $seq,
                                                  -places => 2 );
           or:

           $calc = Bio::Tools::pICalculator->new( -pKset => 'string',
                                                  # a Bio::Seq object
                                                  -seq => $seq,
                                                  -places => 1 );

           Constructs a new pICalculator. Arguments are a flattened hash.
           Valid, optional keys are:

           pKset - A reference to a hash with key value pairs for the 
                   pK values of the charged amino acids. Required keys
                   are:

                   N_term   C_term   K   R   H   D   E   C   Y

           pKset - A string ( 'DTASelect' or 'EMBOSS' ) that will 
                   specify an internal set of pK values to be used. The 
                   default is 'EMBOSS'

           seq - A Bio::Seq sequence object to analyze

           places - The number of decimal places to use in the
                    isoelectric point calculation. The default is 2.

 Returns : The description
 Args    : The description or none