view release on metacpan or  search on metacpan

lib/InSilicoSpectro/InSilico/IsoelPoint.pm  view on Meta::CPAN

		 iterative=>['Lehninger', 'Sillero', 'Rodwell', 'EMBOSS', 'Solomon', 'Patrickios'],
		);

sub new {
  my ($pkg,%h)=@_;
  my $pi={};
  bless $pi, $pkg;

#### Check Lehninger's values
  $pi->{pK}{Sillero}=
    {'Carboxyl'=>3.2,D=>4.0,E=>4.5,'Amino'=>8.2,K=>10.4,R=>12.0,H=>6.4,C=>9.0,Y=>10.0}; # Sillero
  $pi->{pK}{Rodwell}=
    {'Carboxyl'=>3.1,D=>3.86,E=>4.25,'Amino'=>8.0,K=>11.5,R=>11.5,H=>6.0,C=>8.33,Y=>10.07}; # Rodwell
  $pi->{pK}{Lehninger}=
    #    {'Carboxyl'=>3.1,D=>4.4,E=>4.4,'Amino'=>8.0,K=>10.0,R=>12,H=>6.5,C=>8.5,Y=>10.0}; # DTASelect (Lehninger)
    {'Carboxyl'=>2.34,D=>3.86,E=>4.25,'Amino'=>9.69,K=>10.5,R=>12.4,H=>6.0,C=>8.33,Y=>10.0}; # Lehninger
  $pi->{pK}{EMBOSS}=
    {'Carboxyl'=>3.6,D=>3.9,E=>4.1,'Amino'=>8.6,K=>10.8,R=>12.5,H=>6.5,C=>8.5,Y=>10.1}; # EMBOSS
  $pi->{pK}{Solomon}=
    {'Carboxyl'=>2.4,D=>3.9,E=>4.3,'Amino'=>9.6,K=>10.5,R=>12.5,H=>6.0,C=>8.3,Y=>10.1}; # Solomon
  $pi->{pK}{Patrickios}=
    {'Carboxyl'=>'A',D=>'A',E=>'A',K=>'B',R=>'B','Amino'=>'B',pKa=>4.2,pKb=>11.2}; # Patrickios' method
#### Check Lehninger's values

  $pi->set('method','iterative');
  $pi->set('current','Lehninger');
  $pi->set('peptide','');# Peptide
  $pi->set('data',{expseqs=>[],exptimes=>[]});# Experimental data

  $pi->set('calfile','-');# default file

  $pi->set('calibrator',InSilicoSpectro::InSilico::ExpCalibrator->new(fitting=>'bypass',
							 expseqs=>[],expdata=>[],prdata=>[]));# new calibrator

  foreach (keys %h) { $pi->set($_, $h{$_}) }
  return($pi);
}


sub predict {

  my ($this,%h)=@_;
  foreach (keys %h) { $this->set($_, $h{$_}) }

  return $this->{calibrator}->fit($this->predictor);
}

sub predictor {

  my ($this,%h)=@_;
  foreach (keys %h) { $this->set($_, $h{$_}) }

  my $method=$this->{method};

  return $this->$method;
}


sub Patrickios {
# Isoelectric point prediction
# Patrickios' formula

  my ($this,%h)=@_;
  foreach (keys %h) { $this->set($_, $h{$_}) }

  my $pKa=$this->{pK}{Patrickios}{pKa};
  my $pKb=$this->{pK}{Patrickios}{pKb};
  my $A=1; # Carboxyl
  my $B=1; # Amino

  foreach (split '',$this->{peptide}) {
    if (exists $this->{pK}{Patrickios}{$_}) {
      $A++ if $this->{pK}{Patrickios}{$_} eq 'A';
      $B++ if $this->{pK}{Patrickios}{$_} eq 'B';
    }
  }
  my $R=$A/$B;
  return($pKa - log10((1/2)*((1-$R)/$R)+sqrt((1-$R)*(1-$R)/$R/$R+(4/$R)*10**($pKa-$pKb))));
}

sub iterative {
# Isoelectric point prediction
# Iterative

  my ($this,%h)=@_;

  my $charge;
  my $pH;
  my $step=0.5;
  my %count;
  my @aa;
  my ($gamma,$last_charge,$error); 

  foreach (split '',$this->{peptide}) { $count{$_}++ if /^[KRHDECY]$/  }
  $count{'Amino'}=1;
  $count{'Carboxyl'}=1;

  $pH=7;
  $step=3.5;
  $last_charge=0;
  $gamma=0.00001;

  do {
    $charge=0;
    foreach (keys %count) {
      $charge+= $count{$_}*pcharge($this->{pK}{$this->{current}}{$_},$pH) if /^[KRH]$|Amino/;
      $charge-= $count{$_}*pcharge($pH,$this->{pK}{$this->{current}}{$_}) if /^[DECY]$|Carboxyl/;
    }
    ($charge > 0)? ($pH+=$step) : ($pH-=$step);
    $step/=2;
    $error=abs($charge-$last_charge);
    $last_charge=$charge;
 #   print "$pH $charge $error\n";
    } until $error < $gamma;

  return $pH;
}

sub pcharge {

  my $val=10**($_[1] - $_[0]);