#!/usr/bin/perl

# Randal R. Ketchem, Ph.D.                619.784.9879 (voice1)
# Department of Molecular Biology, MB-9   619.784.8754 (voice2)
# The Scripps Research Institute          619.784.9985 (FAX)
# 10550 N. Torrey Pines Road              ketchemr@scripps.edu (email)
# La Jolla, CA 92037-1027                 http://www.scripps.edu/~ketchemr

# modified 04/05/99 by Patty Fagan Jones to filter for multiply assigned peaks only.

# Create 7 column file for DIANA of multiple assignments from Felix97 using the
# volume and peak files generated from felix.  This script as written will only
# read the first two multiple assignments on one peak.  If there are more than
# two assignments they will be ignored.  Input peak file must be exported from
# Felix97 in "ASCII Printout All Peaks" format to retain the multiple assignment
# information.  The first line of the volume and peak files must be deleted before
# running mkmult_assig.  The output of this script is a file which should be called
# mult_assig, and is input for makerst.new 
# 

# Bin file
# 0.103E+09 3.4
# 0.763E+07 4.2
# 0.0       5.5

# Volume file
#    14     0.622E+08

# Peak file
# 609    8.587    3.905  118.536 31.270 60.805 37.568  2 1:ALA_49:HN      1:THR_82:HB      1:ALA_49:N       1:ALA_49:HN      1:ALA_49:HA      1:ALA_49:N

require "ctime.pl";
require "getopt.pl";

if($#ARGV < 7)
{
   $prgName = `basename $0`;
   chop($prgName);
   die <<EOHELP;

This program will create a 7 column restraint file for DIANA from peaks
which were multiply assigned in Felix97.  It uses the same input files
as makerst.new, except that the format of the peak file from Felix97 must
be changed to "ASCII Printout All Peaks" to retain multiple assignment
information. All distances for restraints from multiply assigned peaks
are set to 6.0 angstroms.  It reads a stereo file to get the stereospecific
assignments.

Volume file:
    14     0.622E+08     0.753E+06     0.753E+06

Peak File:
609    8.587    3.905  118.536 31.270 60.805 37.568  2 1:ALA_49:HN      1:THR_82:HB      1:ALA_49:N       1:ALA_49:HN      1:ALA_49:HA      1:ALA_49:N

Bin file:
#0.103E+09 3.4
#0.763E+07 4.2
#0.0       5.5

Stereo file:
29 LEU HD11 MD1
29 LEU HD12 MD2
54 VAL HG11 MG1
54 VAL HG12 MG2

Usage: $prgName -b [bin file] -c [volume column] -p [peak file] -s [stereo file] -v [volume file]
Volume column must be 1..[number of columns]. The peak number column is 0.
EOHELP
}

&Getopt('bcpsv');
$binFile = $opt_b;
$column = $opt_c;
$peakFile = $opt_p;
$commentPeakFile = $peakFile;
$commentPeakFile =~ s/.*\///;
$stereoFile = $opt_s;
$volumeFile = $opt_v;

# Store the bins and their distances in a hash.
open(BINFILE, $binFile) || die("Cannot open $binFile.\n");
while(<BINFILE>)
{
   next if(/^\s*#/);
   next if(/^\s*$/);
   chop;
   @words = split(' ');
   $binNum{$words[0]} = $words[1];
}
close(BINFILE);
die("No volume bins in bin file $binFile.\n") unless(%binNum);

# Store the stereos in a hash.
if($stereoFile)
{
   open(STEREOFILE, $stereoFile) || die("Cannot open $stereoFile.\n");
   while(<STEREOFILE>)
   {
      next if(/^\s*#/);
      next if(/^\s*$/);
      chop;
      ($resNum, $resName, $atomName, $stereoName) = split(' ');
      $stereo{"$resNum $resName $atomName"} = $stereoName;
   }
   close(STEREOFILE);
}

# Store all the peak numbers and their volumes in a hash.
open(VOLUMEFILE, $volumeFile) || die("Cannot open $volumeFile.\n");
while(<VOLUMEFILE>)
{
   next if(/^\s*#/);
   next if(/^\s*$/);
   next if(/^[a-z]|^[A-Z]/);
   chop;
   @words = split(' ');
   $numCols = scalar(@words);
   die("The column number must be [1..$numCols]")
         if($column < 1 || $column > $numCols);
   $volumeNum{$words[0]} = $words[$column];
}
close(VOLUMEFILE);
die("No volumes in volume file $volumeFile.\n") unless(%volumeNum);

# Loop through all of the peaks.
# This is the section which was edited by Patty Fagan Jones
# 7/14/99 to filter for multiply assigned peaks
open(PEAKFILE, $peakFile) || die("Cannot open $peakFile.\n");
while(<PEAKFILE>)

{
   next if(/^\s*#/);		#skip commented lines
   next if(/^\s*$/);		#skip blank lines
   chop;
   @words = split(' ');
   $peakNum = $words[0];
   $multiTag = $words[7];

   next if($multiTag == 1);	#skip all singly assigned peaks

   $resInfo1 = $words[8];
   $resInfo2 = $words[9];
   $resInfo3 = $words[11];
   $resInfo4 = $words[12];

   @res = split(':', $resInfo1);	      # 1:ILE_13:HN -> 1 ILE_13 HN
   $chainID1 = $res[0];
   $atomName = $res[2];
   @splitRes = split('_', $res[1]);	      # ILE_13 -> ILE 13
   $resName1 = substr($splitRes[0], 0, 3);    # Grab first three characters of the residue name
   $resNum1 = $splitRes[1];
   $atomName1 = &ConvertAtomName($resNum1, $resName1, $atomName);

   @res = split(':', $resInfo2);
   $chainID2 = $res[0];
   $atomName = $res[2];
   @splitRes = split('_', $res[1]);
   $resName2 = substr($splitRes[0], 0, 3);
   $resNum2 = $splitRes[1];
   $atomName2 = &ConvertAtomName($resNum2, $resName2, $atomName);

   @res = split(':', $resInfo3);
   $chainID3 = $res[0];
   $atomName = $res[2];
   @splitRes = split('_', $res[1]);
   $resName3 = substr($splitRes[0], 0, 3);
   $resNum3 = $splitRes[1];
   $atomName3 = &ConvertAtomName($resNum3, $resName3, $atomName);

   @res = split(':', $resInfo4);
   $chainID4 = $res[0];
   $atomName = $res[2];
   @splitRes = split('_', $res[1]);
   $resName4 = substr($splitRes[0], 0, 3);
   $resNum4 = $splitRes[1];
   $atomName4 = &ConvertAtomName($resNum4, $resName4, $atomName);

# Make sure the volume exists.
   if($volumeNum{$peakNum} == 0)
   {
      die <<EONOVOL;

Argh. I'm melting! I'm melting!
No volume for peak number $peakNum.
$resNum1 $resName1 $atomName1 $resNum2 $resName2 $atomName2


EONOVOL
   }

# Divide methyl volumes by 3, etc.
   if($atomName1 =~ /MB|MG|QQG|MD|QQD|ME/) { $volumeNum{$peakNum} /= 3.0; }
   if($atomName1 =~ /CG|CZ/) { $volumeNum{$peakNum} /= 2.0; }
   if($atomName2 =~ /MB|MG|QQG|MD|QQD|ME/) { $volumeNum{$peakNum} /= 3.0; }
   if($atomName2 =~ /CG|CZ/) { $volumeNum{$peakNum} /= 2.0; }

# Assign a distance of 6.0 to each multiple assignment since we cannot know
# the intensity of each of the two or more contributing peaks.

#   foreach $bin (sort bynumber keys(%binNum)) { 
    if($volumeNum{$peakNum} >= $bin) { $distance = 6.0; }
# } 

# Only use intra assignments.
   if($chainID1 eq $chainID2)
   {

# Build an array of peak numbers for printing later.
      $peaknumber = sprintf("%4d", $peakNum);
      push(@peaknumbers, $peaknumber);
# Build a 7col array for printing later. Save the peak number as a comment.
      $line = sprintf("%3d %-5s%-5s%4d %-5s%-5s%7.2f \n%3d %-5s%-5s%4d %-5s%-5s%7.2f \n",
            $resNum1, $resName1, $atomName1, $resNum2, $resName2, $atomName2, $distance,
	    $resNum3, $resName3, $atomName3, $resNum4, $resName4, $atomName4, $distance
            );
      push(@lines, $line);
   }
}
close(PEAKFILE);

# Print some header info.
$date = &ctime(time);
print "# Date: $date";
#print "# Bin File: $binFile\n";
#print "# Stereo File: $stereoFile\n" if($stereoFile);
#print "# Peak File: $peakFile\n";
#print "# Volume File: $volumeFile\n";
#print "# resNum1 resName1 atomName1 ";
#print "resNum2 resName2 atomName2 distance # peakFile:peakNum\n";

# Print the peak number and NOESY name for each multiply assigned peak
#  foreach $peaknumber (@peaknumbers)
#  { print "# Peak number $peaknumber, taken from peak file $peakFile\n"; }
# Print the sorted peaks.
#  foreach $line (@lines) 
#  { print $line; }

# Print the peak number and NOESY name for each multiply assigned peak, followed
# by the multiple peak assignments.
for ($i = 0; $i < @peaknumbers; $i++) 
{
   print "# Peak number $peaknumbers[$i], taken from peak file $peakFile\n";
   print $lines[$i];
}

# Subroutine to sort by number as opposed to character.
 sub bynumber { $a <=> $b; }

# Subroutine to convert the atom names into names disgeo will understand.
sub ConvertAtomName
{
   local($resNum, $resName, $oldName) = @_;
# Convert to upper case.
   $oldName = "\U${oldName}\E";
# If no conversion is done, the atom will be surrounded by 'xxx'. Convert
# everything!
   local($newName) = "xxx${oldName}xxx";


   if($resName1 eq "ARG") { $resName1 = "ARG+" };
   if($resName1 eq "ASPN") { $resName1 = "ASP-" };
   if($resName1 eq "GLUN") { $resName1 = "GLU-" };
   if($resName1 eq "LYS") { $resName1 = "LYS+" };

   if($resName2 eq "ARG") { $resName2 = "ARG+" };
   if($resName2 eq "ASPN") { $resName2 = "ASP-" };
   if($resName2 eq "GLUN") { $resName2 = "GLU-" };
   if($resName2 eq "LYS") { $resName2 = "LYS+" };

   if($resName eq "ARG") { $resName = "ARG+" };
   if($resName eq "ASPN") { $resName = "ASP-" };
   if($resName eq "GLUN") { $resName = "GLU-" };
   if($resName eq "LYS") { $resName = "LYS+" };



   if($resName eq "ALA")
   {
      if($oldName eq "HN") { $newName = "HN" };
      if($oldName eq "HA") { $newName = "HA" };
      if($oldName eq "HB1") { $newName = "QB" };
      if($oldName eq "HB2") { $newName = "QB" };
      if($oldName eq "HB3") { $newName = "QB" };
      if($oldName eq "HB*") { $newName = "QB" };
   }
   if($resName eq "ARG+")
   {
      if($oldName eq "HN") { $newName = "HN" };
      if($oldName eq "HA") { $newName = "HA" };
      if($oldName eq "HB1") { $newName = "QB" };
      if($oldName eq "HB2") { $newName = "QB" };
      if($oldName eq "HB*") { $newName = "QB" };
      if($oldName eq "HG1") { $newName = "QG" };
      if($oldName eq "HG2") { $newName = "QG" };
      if($oldName eq "HG*") { $newName = "QG" };
      if($oldName eq "HD1") { $newName = "QD" };
      if($oldName eq "HD2") { $newName = "QD" };
      if($oldName eq "HD*") { $newName = "QD" };
      if($oldName eq "HE") { $newName = "HE" };
# Missing all past delta
   }
   if($resName eq "ASN")
   {
      if($oldName eq "HN") { $newName = "HN" };
      if($oldName eq "HA") { $newName = "HA" };
      if($oldName eq "HB1") { $newName = "QB" };
      if($oldName eq "HB2") { $newName = "QB" };
      if($oldName eq "HD21") { $newName = "QD2" };
      if($oldName eq "HD22") { $newName = "QD2" };
   }
   if($resName eq "ASP-")
   {
      if($oldName eq "HN") { $newName = "HN" };
      if($oldName eq "HA") { $newName = "HA" };
      if($oldName eq "HB1") { $newName = "QB" };
      if($oldName eq "HB2") { $newName = "QB" };
      if($oldName eq "HB*") { $newName = "QB" };
   }
   if($resName eq "CYS")
   {
   }
   if($resName eq "GLN")
   {
      if($oldName eq "HN") { $newName = "HN" };
      if($oldName eq "HA") { $newName = "HA" };
      if($oldName eq "HB1") { $newName = "QB" };
      if($oldName eq "HB2") { $newName = "QB" };
      if($oldName eq "HB*") { $newName = "QB" };
      if($oldName eq "HG1") { $newName = "QG" };
      if($oldName eq "HG2") { $newName = "QG" };
      if($oldName eq "HG*") { $newName = "QG" };
      if($oldName eq "HE21") { $newName = "QE2" };
      if($oldName eq "HE22") { $newName = "QE2" };
   }
   if($resName eq "GLU-")
   {
      if($oldName eq "HN") { $newName = "HN" };
      if($oldName eq "HA") { $newName = "HA" };
      if($oldName eq "HB1") { $newName = "QB" };
      if($oldName eq "HB2") { $newName = "QB" };
      if($oldName eq "HB*") { $newName = "QB" };
      if($oldName eq "HG1") { $newName = "QG" };
      if($oldName eq "HG2") { $newName = "QG" };
      if($oldName eq "HG*") { $newName = "QG" };
   }
   if($resName eq "GLY")
   {
      if($oldName eq "HN") { $newName = "HN" };
      if($oldName eq "HA1") { $newName = "QA" };
      if($oldName eq "HA2") { $newName = "QA" };
      if($oldName eq "HA*") { $newName = "QA" };
   }
   if($resName eq "HIS")
   {
      if($oldName eq "HN") { $newName = "HN" };
      if($oldName eq "HA") { $newName = "HA" };
      if($oldName eq "HB1") { $newName = "QB" };
      if($oldName eq "HB2") { $newName = "QB" };
      if($oldName eq "HD2") { $newName = "HD2" };
      if($oldName eq "HE1") { $newName = "HE1" };
# Missing all past delta
   }
   if($resName eq "ILE")
   {
      if($oldName eq "HN") { $newName = "HN" };
      if($oldName eq "HA") { $newName = "HA" };
      if($oldName eq "HB") { $newName = "HB" };
      if($oldName eq "HG11") { $newName = "QG1" };
      if($oldName eq "HG12") { $newName = "QG1" };
      if($oldName eq "HG1*") { $newName = "QG1" };
      if($oldName eq "HG21") { $newName = "QG2" };
      if($oldName eq "HG22") { $newName = "QG2" };
      if($oldName eq "HG23") { $newName = "QG2" };
      if($oldName eq "HD1") { $newName = "QD" };
      if($oldName eq "HD2") { $newName = "QD" };
      if($oldName eq "HD3") { $newName = "QD" };
      if($oldName eq "HD*") { $newName = "QD1" };   
      if($oldName eq "HG2*") { $newName = "QG2" };
   }
   if($resName eq "LEU")
   {
      if($oldName eq "HN") { $newName = "HN" };
      if($oldName eq "HA") { $newName = "HA" };
      if($oldName eq "HB1") { $newName = "QB" };
      if($oldName eq "HB2") { $newName = "QB" };
      if($oldName eq "HB*") { $newName = "QB" };
      if($oldName eq "HG") { $newName = "HG" };
      if($oldName eq "HD11") { $newName = "QD1" };
      if($oldName eq "HD12") { $newName = "QD1" };
      if($oldName eq "HD13") { $newName = "QD1" };
      if($oldName eq "HD21") { $newName = "QD2" };
      if($oldName eq "HD22") { $newName = "QD2" };
      if($oldName eq "HD23") { $newName = "QD2" };
      if($oldName eq "HD2*") { $newName = "QD2" };
      if($oldName eq "HD1*") { $newName = "QD1" };
      if($oldName eq "HD*") { $newName = "QQD" };
   }
   if($resName eq "LYS+")
   {
      if($oldName eq "HN") { $newName = "HN" };
      if($oldName eq "HA") { $newName = "HA" };
      if($oldName eq "HB1") { $newName = "QB" };
      if($oldName eq "HB2") { $newName = "QB" };
      if($oldName eq "HB*") { $newName = "QB" };
      if($oldName eq "HG1") { $newName = "QG" };
      if($oldName eq "HG2") { $newName = "QG" };
      if($oldName eq "HG*") { $newName = "QG" };
      if($oldName eq "HD1") { $newName = "QD" };
      if($oldName eq "HD2") { $newName = "QD" };
      if($oldName eq "HD*") { $newName = "QD" };
      if($oldName eq "HE1") { $newName = "QE" };
      if($oldName eq "HE2") { $newName = "QE" };
      if($oldName eq "HE*") { $newName = "QE" };
# Missing NH3
   }
   if($resName eq "MET")
   {
      if($oldName eq "HN") { $newName = "HN" };
      if($oldName eq "HA") { $newName = "HA" };
      if($oldName eq "HB1") { $newName = "QB" };
      if($oldName eq "HB2") { $newName = "QB" };
      if($oldName eq "HB*") { $newName = "QB" };
      if($oldName eq "HG1") { $newName = "QG" };
      if($oldName eq "HG2") { $newName = "QG" };
      if($oldName eq "HG*") { $newName = "QG" };
      if($oldName eq "HE1") { $newName = "QE" };
      if($oldName eq "HE2") { $newName = "QE" };
      if($oldName eq "HE3") { $newName = "QE" };
      if($oldName eq "HE*") { $newName = "QE" }; 
   }
   if($resName eq "PHE")
   {
      if($oldName eq "HN") { $newName = "HN" };
      if($oldName eq "HA") { $newName = "HA" };
      if($oldName eq "HB1") { $newName = "QB" };
      if($oldName eq "HB2") { $newName = "QB" };
      if($oldName eq "HB*") { $newName = "QB" };
      if($oldName eq "HD1") { $newName = "CG" };
      if($oldName eq "HD2") { $newName = "CG" };
      if($oldName eq "HE1") { $newName = "CZ" };
      if($oldName eq "HE2") { $newName = "CZ" };
      if($oldName eq "HZ") { $newName = "HZ" };
      if($oldName eq "HD*") { $newName = "CG" };
      if($oldName eq "HE*") { $newName = "CZ" };   
   }
   if($resName eq "PRO")
   {
      if($oldName eq "HN") { $newName = "HN" };
      if($oldName eq "HA") { $newName = "HA" };
      if($oldName eq "HB1") { $newName = "QB" };
      if($oldName eq "HB2") { $newName = "QB" };
      if($oldName eq "HB*") { $newName = "QB" };
      if($oldName eq "HG1") { $newName = "QG" };
      if($oldName eq "HG2") { $newName = "QG" };
      if($oldName eq "HG*") { $newName = "QG" };
      if($oldName eq "HD1") { $newName = "QD" };
      if($oldName eq "HD2") { $newName = "QD" };
      if($oldName eq "HD*") { $newName = "QD" };
   }
   if($resName eq "SER")
   {
      if($oldName eq "HN") { $newName = "HN" };
      if($oldName eq "HA") { $newName = "HA" };
      if($oldName eq "HB1") { $newName = "QB" };
      if($oldName eq "HB2") { $newName = "QB" };
      if($oldName eq "HB*") { $newName = "QB" };
# Missing HO
   }
   if($resName eq "THR")
   {
      if($oldName eq "HN") { $newName = "HN" };
      if($oldName eq "HA") { $newName = "HA" };
      if($oldName eq "HB") { $newName = "HB" };
      if($oldName eq "HG1") { $newName = "HG1" };
      if($oldName eq "HG21") { $newName = "QG2" };
      if($oldName eq "HG22") { $newName = "QG2" };
      if($oldName eq "HG23") { $newName = "QG2" };
      if($oldName eq "HG2*") { $newName = "QG2" };
# Missing HO
   }
   if($resName eq "TRP")
   {
   }
   if($resName eq "TYR")
   {
      if($oldName eq "HN") { $newName = "HN" };
      if($oldName eq "HA") { $newName = "HA" };
      if($oldName eq "HB1") { $newName = "QB" };
      if($oldName eq "HB2") { $newName = "QB" };
      if($oldName eq "HB*") { $newName = "QB" };
      if($oldName eq "HD1") { $newName = "CG" };
      if($oldName eq "HD2") { $newName = "CG" };
      if($oldName eq "HE1") { $newName = "CZ" };
      if($oldName eq "HE2") { $newName = "CZ" };
      if($oldName eq "HD*") { $newName = "CG" };
      if($oldName eq "HE*") { $newName = "CZ" };
      if($oldName eq "HH") { $newName = "HH" };
# Missing HO
   }
   if($resName eq "VAL")
   {
      if($oldName eq "HN") { $newName = "HN" };
      if($oldName eq "HA") { $newName = "HA" };
      if($oldName eq "HB") { $newName = "HB" };
      if($oldName eq "HB*") { $newName = "HB" };
      if($oldName eq "HG11") { $newName = "QG1" };
      if($oldName eq "HG12") { $newName = "QG1" };
      if($oldName eq "HG13") { $newName = "QG1" };
      if($oldName eq "HG21") { $newName = "QG2" };
      if($oldName eq "HG22") { $newName = "QG2" };
      if($oldName eq "HG23") { $newName = "QG2" };
      if($oldName eq "HG1*") { $newName = "QG1" };
      if($oldName eq "HG2*") { $newName = "QG2" };
      if($oldName eq "HG*") { $newName = "QQG" };
   }

# Stereo's
   if($stereo{"$resNum $resName $oldName"})
   { $newName = $stereo{"$resNum $resName $oldName"} };

# Return the new atom name.
   $newName;


}