#!/usr/bin/perl

#
# Artem G. Evdokimov
# artem@xtals.org
# please feel free to modify/include this script as long as you
# acknowledge the source code back to me.
# Please note that all directories are hard-coded into this script and 
# you MUST change the directories and file locations according to
# what you have on your system.
#

use CGI::Pretty;
$query = new CGI;
print $query->header;
print<<EOT;
<script language="JavaScript">
function framePrint() {
 window.focus();
 window.print();
}
</script>
EOT

print $query->start_html;
$OK_CHARS='a-zA-Z';
my $sequence = $query->param('**Sequence');

$sequence =~ s/[^$OK_CHARS]//go;
$sequence =~ tr/[a-z]/[A-Z]/;

if ($sequence eq '') {
print "<hr>";
print "No useful sequence entered.";
print "<hr>";
exit();
}

&ReadResidues;

$avemass = &CalcAveMass($sequence);
$monmass = &CalcMonMass($sequence);
my $name = $query->param('name');
print "<center>$name</center><BR><HR> \n";
$pI = &CalcpI($sequence);
$EC = &CalcEC($sequence);
$OD = $EC/$avemass;
print "<TABLE BORDER=0 CELLPADDING=3 CELLSPAING=3><TR> \n";
print "<TD ALIGN=RIGHT>pI = </TD><TD>";
printf ("%.2f \n", $pI);
print "</TD></TR>\n";
print "<TD ALIGN=RIGHT>Average mass = </TD><TD>";
printf ("%.4f \n", $avemass);
print "</TD></TR>\n";
print "<TD ALIGN=RIGHT>Monoisotopic mass = </TD><TD>";
printf ("%.4f \n", $monmass);
print "</TD></TR>\n";
print "<TD ALIGN=RIGHT>OD<sub>280</sub> =  </TD><TD>";
printf ("%.2f \n", $OD);
print "</TD></TR>\n";
print "<TD ALIGN=RIGHT>Extinction coefficient = </TD><TD>";
printf ("%.1f \n", $EC);
print "</TD></TR></TABLE>\n";
print "<P>\n";
print "<font  face=courier size=3><pre>\n";
print "           |......... |......... |......... |......... |......... |......... \n \n";
 
my $schar;
my @schars = split (//, $sequence);
my $scount = 1;
 printf(" %8d ",$scount);
 foreach $schar (@schars) {
  if (($scount%60 == 1)&& !($scount == 1)) {
   my $zcount= $scount-1;
	 print " $zcount \n";
   printf(" %8d ",$scount);
  }
  if ($scount%10 == 1)  {
   print " ";
  }
  print "$schar";
  $scount++;
 }

print "</pre></font>";
print "<P> \n";
 
print <<EOT;
<center><img src="$s_grafile"></center>
<input type="button" value="Print This" onClick="framePrint()">
EOT

print $query->end_html;
exit();

# - - - SUBROUTINES - - - 

sub CalcEC
{
 my $seq = @_[0];
 my $char;
 my @chars = split (//, $seq);
 my %count;
 my $EC = 0;
 foreach $char (@chars) {
  $count{$char}++;
 }
 foreach $char (keys (%count)) {
   if ($char eq 'C') {
    $EC = $EC + $count{$char} * 48;
   }
   if ($char eq 'Y') {
    $EC = $EC + $count{$char} * 1280;
   }
   if ($char eq 'W') {
    $EC = $EC + $count{$char} * 5690;
   }
  }
 return ($EC); 
}

# - - -

sub CalcpI {

 my $filebase = time() . '.junk';

 # substitute MY HTML PATH for the relevant directory in your system
 my $filename = '../MY HTML PATH/ermine_tmp/' . $filebase . '.txt';
 $grafile = $filename . '.png';

# usually one directory up from cgi-bin will give you the Web Root, however
# you might want to hard-code the directory, if it's not default -
# something like this 
# $s_grafile = 'www.mysite.org/' . 'ermine_tmp/' . $filebase . '.txt.png';
  $s_grafile = '../ermine_tmp/' . $filebase . '.txt.png';

 open (FILE, ">$filename");
 
 my $seq = @_[0];
 my $char;
 my @chars = split (//, $seq);
 my %count;
 my $pka_side, $cr, $c, $cn, $charge;
 my $first = @chars[0];
 my $last = @chars[$#chars];
 my $ph, $dph;
 my $reached = 0;
 my $pI = 0;
 
 $pk_c = 0 - $pKa_C{$last};
 $pk_n = $pKa_N{$first};
 
 foreach $char (@chars) {
  $count{$char}++;
 }
 
# print "$first, $pk_n, $last, $pk_c \n";

# begin pH loop from 0 to 14

 for ($ph = 0; $ph<14.5; $ph = $ph+0.1) {
  $charge = 0;
  
  $c = 10**($pk_n-$ph);
  $cr = $c/(1+$c);
  $cn = (+1) * $cr;
  # printf ("N -> %.2f ",$cn);
  $charge = $charge + $cn;
   
  $c = 10**($pk_c+$ph);
  $cr = $c/(1+$c);
  $cn = (-1) * $cr;
  # printf ("C -> %.2f \n",$cn);
  $charge = $charge + $cn;
  
  foreach $char (keys (%count)) {
     
   if ($char =~ m/K|R|H/) {
    $pka_side = $pKa_side{$char};
   } else {
    $pka_side = - $pKa_side{$char};
   }
  
   if ($pka_side > 0) {
    $c = 10**($pka_side-$ph);
    $cr = $c/(1+$c);
    $cn = (+1) * $count{$char} * $cr;
   } elsif ($pka_side < 0) {
    $c = 10**($pka_side+$ph);
    $cr = $c/(1+$c);
    $cn = (-1) * $count{$char} * $cr;
   } else {
    $cn = 0;
   }
   $charge = $charge + $cn;
  }
  
  $val = sprintf ("%5.2f  %5.5f \n",$ph,$charge);
  print (FILE $val);
  if (($charge <= 0) && ($reached == 0))  {
   $pI = $ph - 0.05;
   $reached = 1;
  }     
 } 
 close (FILE);
 
# please note that it is adviseable to set up .crontab to clean up
# junk graphics file directory once daily, to avoid storing all the 
# little .png files
 eval {
 open GNUPLOT, "| ../gnuplot";
 print (GNUPLOT <<EOT);
 set term png small color
 set size 0.4
 set out "$grafile"
 set nokey
 set xlabel "pH"
 set ylabel "Charge"
 plot [0:14] "$filename" with linespoints
 quit
EOT
 close (GNUPLOT);
 } ;

 print `rm -rf $filename`; 
 return ($pI);
}

#-------------------------------------------------


sub CalcAveMass {
 my $seq = @_[0];
 my $mass = 18.0152;
 my $char;
 my @chars = split (//, $seq);
 foreach $char (@chars) {
  $mass = $mass + $mass_p_av{$char};
 }
 return ($mass);
}

#-------------------------------------------------

sub CalcMonMass {
 my $seq = @_[0];
 my $mass = 18.0105;
 my $char;
 my @chars = split (//, $seq);
 foreach $char (@chars) {
  $mass = $mass + $mass_p_mo{$char};
 }
 return ($mass);
}

#-------------------------------------------------

sub ReadResidues {
 my $filename = 'data/residue.dat';
 open (FI, $filename);
 my $name1, $m1;
 
 while (!eof(FI)) {
  $tq = new CGI(\*FI);
  $name1 = $tq->param('**Name_1');
#  print "$name1 ";
  push (@names, $name1);
  
  $m1 = $tq->param('**Mass_pep_average');
  $mass_p_av{$name1} = $m1;
#  print "$m1 ";
  $m1 = $tq->param('**Mass_pep_mono');
  $mass_p_mo{$name1} = $m1;
#  print "$m1 ";
  $m1 = $tq->param('**Mass_average');
  $mass_av{$name1} = $m1;
#  print "$m1 ";
  $m1 = $tq->param('**Mass_mono');
  $mass_mo{$name1} = $m1;
#  print "$m1 \n";
  $m1 = $mass_av{$name1} - $mass_p_av{$name1};
#  print "$m1 ";
  $m1 = $mass_mo{$name1} - $mass_p_mo{$name1};
#  print "$m1 \n";
  $m1 = $tq->param('**pKa_C');
  $pKa_C{$name1} = $m1;
  $m1 = $tq->param('**pKa_N');
  $pKa_N{$name1} = $m1;
  $m1 = $tq->param('**pKa_side');
  $pKa_side{$name1} = $m1;
 }
 close (FI);
}