#!/usr/local/bin/perl -w
#given names of input qhull results file and pdb file and name of output
#pdb file filters the pdb file to remove the points or amino acids in the
#convex hull producing pdb file with less atom/hetam entries
 
open (QHULLRESULT,"<".$ARGV[0])||die "qhull resultfile could not be opened in rebuild.prl";

open (PDBSOURCE,"<".$ARGV[1])||die "prior pdb could not be opened in rebuild.prl";

open (PDBNEW,">".$ARGV[2])||die "PDB output could not be opened in rebuild.prl";

#assing this to initial pdb file name
$_=$ARGV[3];
/data\/pdb(.*)\.ent/;
#strip off directory bc this is used to create results file names
$resultsfile=$1;

open (SHELLS,">>"."results/".$resultsfile.".shells")||die " shell file open failed in rebuild.prl";
print SHELLS "\n";
print SHELLS $ARGV[1];
print SHELLS " ";


# unremm this line and rem out below to remove only immediate atoms 
# this version hasn't been thoroughly tested
#    prunebyatom();

prunebyacid($ARGV[1]);

close PDBSOURCE;
close PDBNEW;
close QHULLRESULT;
close SHELLS;
# end main program

sub prunebyacid{
	# remove whole amino acids that appear in convex hull, simples way
	# seems two pass, once through to get amino acid # and the chain letter or just het# for each
	# atom # then a second pass to remove all atoms with that amino acid #/chain letter 
	@entriestoremove={};
	$x=0;
	$y=0;
	<QHULLRESULT>; #ditch the number of points in the hull
	
	while (<QHULLRESULT>){
		if (/(\d+)/){
			$targetentry=$1;
			
			while ($x<=$targetentry){
				 #print "targetentry is $targetentry";
				 #print "x is $x";
				$thedataline=<PDBSOURCE>;
				#if the current line is position record
				$_=$thedataline;
				#print "$thedataline \n";
				#if ($_=undef){die;}
				if((/ATOM\s+.+\s+(\w\s+\d+)\s+(-?\d+\.\d+)\s+(-?\d+\.\d+)\s+(-?\d+\.\d+)\s+\d+.+\d+.+\d+.+/) ||
				(/HETATM\s+.+\s+(\d+)\s+(-?\d+\.\d+)\s+(-?\d+\.\d+)\s+(-?\d+\.\d+)\s+\d+.+\d+.+\d+.+/) ||
				   (/ATOM\s+.+\s+(\d+)\s+(-?\d+\.\d+)\s+(-?\d+\.\d+)\s+(-?\d+\.\d+)\s+\d+.+\d+.+\d+.+/)
				    ){
			
					if ($x==$targetentry){
					    if ($entriestoremove[$y-1] ne $1){ 
						#store proper digit in array
						$entriestoremove[$y]=$1;
						#print "stored target ".$1."\n";
						$y++;
					    }
					}
					$x++;
				} # end if position line
			    } #end while <target

		    } #end if digit

	    } #end while of ones to remove

	close PDBSOURCE;
	$removeindex=0;
	$inmatcharea=0;

#now array should have target base #'s and we pass through taking them out
	open (PDBSOURCE,"<".$_[0])||die "prior pdb could not be opened in rebuild.prl";

	while (<PDBSOURCE>){
	   	 $thedata=$_;

		if((/ATOM\s+.+\s+(\w\s+\d+)\s+(-?\d+\.\d+)\s+(-?\d+\.\d+)\s+(-?\d+\.\d+)\s+\d+.+\d+.+\d+.+/)||
		   (/HETATM\s+.+\s+(\d+)\s+(-?\d+\.\d+)\s+(-?\d+\.\d+)\s+(-?\d+\.\d+)\s+\d+.+\d+.+\d+.+/)||
		   (/ATOM\s+.+\s+(\d+)\s+(-?\d+\.\d+)\s+(-?\d+\.\d+)\s+(-?\d+\.\d+)\s+\d+.+\d+.+\d+.+/)
		   ){   
		    $currentlinedata=$1;

		    # if it is current target dont print it and move to next target
		    # note if its the calcium
		    if ($currentlinedata eq $entriestoremove[$removeindex]){
			if ($inmatcharea!=1){
			    $_=$thedata;
			    /\S+\s+\d+\s+\S+\s+(\S+)\s+.*/;
			    #print $1;
			    #print what's gotten removed to file 
			    print SHELLS "$1 ";
			}
			$inmatcharea=1;
			#print "match target";
			#keep it from incrementing past end when done
	       
			$_=$thedata;
			
			if (/\sCA\s+CA\s/){
				print "hit ca\n";
				open (RESULTS,">>"."results/".$resultsfile.".res")||die "result file open failed in rebuild.prl";
				print RESULTS $ARGV[2];
				print RESULTS "\n";
				close RESULTS;
			}
		    } #end if matches target 
		    else{
			if ($inmatcharea==1){
			    #print "remove index incremented";
			    $inmatcharea=0;
			    if ($removeindex+1<$y){
				$removeindex++;
				
			    }
			}
			if ($currentlinedata ne $entriestoremove[$removeindex]){
			    # if it's not one of our target lines just print it
			    # we have to check again because we may have just incremented removeindex
			    print PDBNEW $thedata;
			}
			else {
			#if after incrementing on leaving match area
			# this may now still match ca so we do this check

			$_=$thedata;
			
			if (/\sCA\s+CA\s/){
				print "hit ca";
				open (RESULTS,">>"."results/".$resultsfile.".res")||die "result file open failed in rebuild.prl";
				print RESULTS $ARGV[2]."\n";
				close RESULTS;
			}


			    if ($removeindex+1<$y){
				$removeindex++;
			    }
			}
		    }

		} #end if point entry 
		 #if its not a atom or het data line just print it
		else {
		    print PDBNEW $thedata;
		}
	}# end while pdbsource


} #end sub


sub prunebyatom{
# rebuild pdb file taking out only atoms
$x=0;
<QHULLRESULT>; #ditch the number of points in the hull

while (<QHULLRESULT>){
	if (/(\d+)/){
		$targetentry=$1;
		#print "$targetentry";
		while ($x<=$targetentry){
			$thedataline=<PDBSOURCE>;
			#if the current line is position record
			$_=$thedataline;
			#print "$thedataline";
			if((
				/ATOM\s+.+\d+\s+(-?\d+\.\d+)\s+(-?\d+\.\d+)\s+(-?\d+\.\d+)\s+\d+.+\d+.+\d+.+/)||
				(/HETATM\s+.+\d+\s+(-?\d+\.\d+)\s+(-?\d+\.\d+)\s+(-?\d+\.\d+)\s+\d+.+\d+.+\d+.+/)){
			
				if ($x!=$targetentry){
				#print "print of atom entry";
				print PDBNEW $thedataline;
				}
	
			$x++;	
			}
			#if its not a position record print it to output
			else { 
			#print "print of not data line";
			print PDBNEW $thedataline;
			}
		} #end while <target

	} #end if digit

} #end while of ones to remove

#put everything else into the file
while(<PDBSOURCE>){
    $data=$_;
    print PDBNEW $data;
}

} #end  sub  rebuild on atomic level 

