#!/usr/bin/perl -w

#ContSum.pl sumarize contact and output histogram

# Calculate the contact hitogram for a list of files in pdb format
# the files have twenty types of atoms CA and 19 sidechains.
if(!@ARGV){
	printf "ContSum.pl: calculate contacts.\n";
	printf "example: ContSum.pl -if 1a32.sc.pdb\n";
	exit();
}

$Num_Atom_Type = 20;
$CaOnly=0;   # whether only Ca included
$Length=2000;
$isAB=0;
$isSC=1;
$isEst=0;

foreach $i (0..$#ARGV) {
    if($ARGV[$i] eq "-if") {
	$InFile=$ARGV[++$i];
    }
    if($ARGV[$i] eq "-of") {
        open OUTFILE, ">>$ARGV[$i+1]";
	$i++;
    }
    if($ARGV[$i] eq "-na") {
        $Num_Atom_Type = $ARGV[++$i];
    }
    if($ARGV[$i] eq "-ca") {
        $CaOnly=1;
    }
    if($ARGV[$i] eq "-ab") {
        $isAB=1;
    }
    if($ARGV[$i] eq "-sc"){
	$isSC=1;
    }
    if($ARGV[$i] eq "-e"){
	$isEst=1;
    }
    if($ARGV[$i] eq "-l") {
        $Length = $ARGV[++$i];
    }
}

#@R_sc={1.5,2.0,2.5,3.1,3.4,0.0,3.1,2.3,3.5,2.6,3.0,2.5,1.9,3.1,4.1,1.9,1.9,2.0,3.9,3.8};

$ResLab=0;
$Res{"ALA"}=$ResLab;
$Res{$ResLab++}="ALA";
$Res{"CYS"}=$ResLab;
$Res{$ResLab++}="CYS";
$Res{"ASP"}=$ResLab;
$Res{$ResLab++}="ASP";
$Res{"GLU"}=$ResLab;
$Res{$ResLab++}="GLU";
$Res{"PHE"}=$ResLab;
$Res{$ResLab++}="PHE";
$Res{"GLY"}=$ResLab;
$Res{$ResLab++}="GLY";
$Res{"HIS"}=$ResLab;
$Res{$ResLab++}="HIS";
$Res{"ILE"}=$ResLab;
$Res{$ResLab++}="ILE";
$Res{"LYS"}=$ResLab;
$Res{$ResLab++}="LYS";
$Res{"LEU"}=$ResLab;
$Res{$ResLab++}="LEU";
$Res{"MET"}=$ResLab;
$Res{$ResLab++}="MET";
$Res{"ASN"}=$ResLab;
$Res{$ResLab++}="ASN";
$Res{"PRO"}=$ResLab;
$Res{$ResLab++}="PRO";
$Res{"GLN"}=$ResLab;
$Res{$ResLab++}="GLN";
$Res{"ARG"}=$ResLab;
$Res{$ResLab++}="ARG";
$Res{"SER"}=$ResLab;
$Res{$ResLab++}="SER";
$Res{"THR"}=$ResLab;
$Res{$ResLab++}="THR";
$Res{"VAL"}=$ResLab;
$Res{$ResLab++}="VAL";
$Res{"TRP"}=$ResLab;
$Res{$ResLab++}="TRP";
$Res{"TYR"}=$ResLab;
$Res{$ResLab++}="TYR";

# get the contact distance from file avgDisMatr.txt
if($Num_Atom_Type==40){
    open CONT, "/home/jinfeng/SharedData/avgDisMatr.txt";
}
elsif($Num_Atom_Type==20) {
    open CONT, "/home/jinfeng/SharedData/ContDis.txt";
}
$j=0;
while(defined($line=<CONT>)) {
    chomp($line);
    $line =~ s/^\s+//g;
    @contList=split(/\s+/,$line);
    for($i=0;$i<$Num_Atom_Type;$i++) {
	$contDis->[$j][$i]=$contList[$i];
    }
    ++$j;
}

if($InFile=~ /(.+)\.pdb/) {
    $name=$1;
}
else{$name=$InFile;}

if($isAB==1){
    $filename=$name . ".ab.pdb";
    $temp = $name . ".ab.abc";
    open OUTFILE, ">$temp";
}
elsif($isSC==1){
    if($isEst==0){
	$filename=$name . ".sc.pdb";
	$temp = $name . ".sc.abc";
	open OUTFILE, ">$temp";
    }
    elsif($isEst==1){
	$filename=$name . ".sce.pdb";
	$temp = $name . ".sce.abc";
	open OUTFILE, ">$temp";
    }
}

if(-e "$filename") {
    open PDBFILE, "$filename";
}
elsif(-e "/home/jinfeng/SharedData/pc30_res2.0_SC/$filename") {
    open PDBFILE, "/home/jinfeng/SharedData/pc30_res2.0_SC/$filename";
}
elsif(-e "/home/jinfeng/SharedData/pc30_res2.0_SC_E/$filename") {
    open PDBFILE, "/home/jinfeng/SharedData/pc30_res2.0_SC_E/$filename";
}
else {
    printf "$filename is not found!\n";
    exit;
}

$lineCnt=0;
while(<PDBFILE>){
    $line[$lineCnt++]=$_;
}
close PDBFILE;
$initLength = substr($line[0], 22,4); $initLength =~ s/^\s+|\s+$//g;
$tLength=$initLength+$Length;
for($i=0;$i<$lineCnt;++$i) {
    if(substr($line[$i],16,4) =~ /[BCDE234][A-Z][A-Z][A-Z]/) {next;}
    $atomName= substr($line[$i], 12,4); $atomName=~ s/^\s+|\s+$//g; 
    $resName = substr($line[$i], 17,3); $resName =~ s/^\s+|\s+$//g;
    $resName =~ tr/a-z/A-Z/;
    $resNum  = substr($line[$i], 22,4); $resNum  =~ s/^\s+|\s+$//g;
    if($resNum==$tLength) {
	printf "Reach target length $Length $tLength\n";
	last;
    }
    $chainID = substr($line[$i], 21,1); $chainID = "*" if($chainID eq " ");
    $x       = substr($line[$i], 30,8); $x       =~ s/^\s+|\s+$//g;
    $y       = substr($line[$i], 38,8); $y       =~ s/^\s+|\s+$//g;
    $z       = substr($line[$i], 46,8); $z       =~ s/^\s+|\s+$//g;
    
    if(($atomName eq "CB" ||$atomName eq "SC")&&$CaOnly==1) {
	next;
    }
    for($j=$i+3;$j<$lineCnt;++$j) {
	if(substr($line[$j],16,4) =~ /[BCDE234][A-Z][A-Z][A-Z]/) {next;}
	$atomName2= substr($line[$j], 12,4); $atomName2=~ s/^\s+|\s+$//g; 
	if(($atomName2 eq "CB" ||$atomName2 eq "SC")&&$CaOnly==1) {
	    next;
	}
	$resName2 = substr($line[$j], 17,3); $resName2 =~ s/^\s+|\s+$//g;
	$resName2 =~ tr/a-z/A-Z/;
	$resNum2  = substr($line[$j], 22,4); $resNum2  =~ s/^\s+|\s+$//g;
	if($resNum2==$tLength) {
	    last;  
	}
	$chainID2 = substr($line[$j], 21,1); $chainID2 = "*" if($chainID eq " ");
	$x2       = substr($line[$j], 30,8); $x2       =~ s/^\s+|\s+$//g;
	$y2       = substr($line[$j], 38,8); $y2       =~ s/^\s+|\s+$//g;
	$z2       = substr($line[$j], 46,8); $z2       =~ s/^\s+|\s+$//g;
	$distance=sqrt(($x-$x2)*($x-$x2)+($y-$y2)*($y-$y2)+($z-$z2)*($z-$z2));
	if(($resNum2-$resNum)<3) {next;}
	if($atomName eq "CA") {
	    $sub1=$Res{"GLY"};
	}
	elsif($atomName eq "SC"||$atomName eq "CB") {
	    $sub1=$Res{$resName};
	}

	if($atomName2 eq "CA") {
	    $sub2=$Res{"GLY"};
	}
	elsif($atomName2 eq "SC"||$atomName2 eq "CB") {
	    $sub2 =$Res{$resName2};
	}
	if($isEst==1){
	    if($distance<=(1.6*$contDis->[$sub1][$sub2])) {
		if($atomName eq "CA" && $atomName2 eq "CA"){
		    if($distance<3.0){$disbin=0;}
		    elsif($distance<1.0*$contDis->[$sub1][$sub2]){$disbin=1;}
		    elsif($distance<1.3*$contDis->[$sub1][$sub2]){$disbin=2;}
		    else {$disbin=3;}
		}
		elsif($atomName eq "SC" && $atomName2 eq "SC"){
		    if($distance < 2.0){$disbin=0;}
		    elsif($distance<1.0*$contDis->[$sub1][$sub2]){$disbin=1;}
		    elsif($distance<1.3*$contDis->[$sub1][$sub2]){$disbin=2;}
		    else {$disbin=3;}
		}
		else {
		    if($distance< 2.5 ){$disbin=0;}
		    elsif($distance<1.0*$contDis->[$sub1][$sub2]){$disbin=1;}
		    elsif($distance<1.3*$contDis->[$sub1][$sub2]){$disbin=2;}
		    else {$disbin=3;}
		}
		printf OUTFILE "$resNum $resName $atomName $resNum2 $resName2 $atomName2 %6.3f %d\n", $distance,$disbin;
	    }
	}
	else {
	    if($distance<=(1.5*$contDis->[$sub1][$sub2])) {
		if($atomName eq "CA" && $atomName2 eq "CA"){
		    if($distance<3.0){$disbin=0;}
		    elsif($distance<1.0*$contDis->[$sub1][$sub2]){$disbin=1;}
		    elsif($distance<1.2*$contDis->[$sub1][$sub2]){$disbin=2;}
		    else {$disbin=3;}
		}
		elsif($atomName eq "SC" && $atomName2 eq "SC"){
		    if($distance < 2.0){$disbin=0;}
		    elsif($distance<1.0*$contDis->[$sub1][$sub2]){$disbin=1;}
		    elsif($distance<1.2*$contDis->[$sub1][$sub2]){$disbin=2;}
		    else {$disbin=3;}
		}
		else {
		    if($distance< 2.5 ){$disbin=0;}
		    elsif($distance<1.0*$contDis->[$sub1][$sub2]){$disbin=1;}
		    elsif($distance<1.2*$contDis->[$sub1][$sub2]){$disbin=2;}
		    else {$disbin=3;}
		}
		printf OUTFILE "$resNum $resName $atomName $resNum2 $resName2 $atomName2 %6.3f %d\n", $distance,$disbin;
	    }
	}
    }
} 

close(OUTFILE);
