#!/usr/bin/perl

#Perl script to count number of DDM entries above a user supplied
#cutoff. It also returns the range of DD values. User can also supply 
#ranges of residues to be included
#Input format is MSI rstrnt file (such as the output from Randy
#Ketchem's ddmrstrnt perl program)

#M. Nelson 7/8/98

if($#ARGV != 7)
{
	$prgName = $0;
	$prgName = `basename $prgName`;
	chomp($prgName);
	die <<EODIE;

This program reads an MSI rstrnt file and counts the number of DDM
entries with an absolute value larger than cutoff. It will also
return the range of DDs. The user must supply residue
ranges as well. These ranges (flags a and b) specify which residues
should be included in the count, etc. Only DDs involving one residue
from range a and one residue from range b will be included. 

dd.rstrnt:
	!BIOSYM restraint 1
	!Converted from the reco_comp.ddm.mtx discom matrix.
	!Tue Apr  7 11:39:30 PDT 1998
	#distance
	1:GLY_1:CA      1:SER_5:CA        -4.070   -4.070 1.00 1.00 1000.000
	1:GLY_1:CA      1:GLY_6:CA        -4.450   -4.450 1.00 1.00 1000.000

Usage: $prgName -c cutoff -f dd.rstrnt -a start-end -b start-end [>outfile]
EODIE
}

require "getopt.pl";

&Getopt('cfab');
$cutoff = $opt_c;
$infile = $opt_f;
$range1 = $opt_a;
$range2 = $opt_b;

($start1, $end1) = split('-',$range1);
($start2, $end2) = split('-',$range2);

#go through input file counting entries above cutoff with $sigcounter
#the two range variables are set to unreasonable values at start to
#ensure that get true answer from file, not from starting setting.

$sigcounter = 0;
$highrange = -2000;
$lowrange = 2000;
open(RSTFILE, $infile) || die("Can't open $infile: $!\n");
while(<RSTFILE>){
	@line = split;
	($junk1,$res1,$junk2) = split(":",$line[0]);
	($resname1, $resnum1) = split("_", $res1);
	($junk1,$res2,$junk2) = split(":",$line[1]);
	($resname2, $resnum2) = split("_", $res2);

	if ($resnum1 >= $start1 && $resnum1 <= $end1 && $resnum2 >= $start2 &&
	$resnum2 <= $end2){
		if ($line[2] >= $cutoff || $line[2] <= -$cutoff){
			$sigcounter++;
		}
		if ($line[2] > $highrange){
			$highrange = $line[2];
			$highres1 = $res1;
			$highres2 = $res2;
		} elsif ($line[2] < $lowrange){
			$lowrange = $line[2];
			$lowres1 = $res1;
			$lowres2 = $res2;
		}
	}

}

print "DDs involving one residue in $range1 and one residue in $range2:\n";
print "---------------------------------------------------------------------\n";
print "DDs with absolute value greater than $cutoff: $sigcounter\n";
print "Most positive DD: $highres1, $highres2: $highrange\n";
print "Most negative DD: $lowres1, $lowres2: $lowrange\n";