make_rpf_precision_map.pl

use strict;
use warnings;

die "Usage: perl make_rpf_precision_map.pl <path> <pdb_file> <QL_file>\n" if(scalar @ARGV < 3);
my ($path,$pdb_file,$ql_file) = @ARGV;
chdir($path);

my $bg = '#EFF8FF';
my $but_bg = '#265580';
my $but_fg = $bg;
my $cbt_act = '#D1E4F2';
my $but_act = '#2D4455';

&calc_rasmol($pdb_file, $ql_file, 'summary');
sub calc_rasmol
{
    my ($pdb_file, $ql_file, $key) = @_;
  #  my $basepath = Cwd::cwd;
    
    print "calc_rasmol: pdb_file = $pdb_file, ql_file = $ql_file, key = $key\n";
    
    unless (-f $pdb_file and -f $ql_file)
    {
	print ("Invalid file input for pdb and ql files\n");
	return;
    }
    
    my $filter = 'H';

    my @count_array = &get_count_array($ql_file, $filter, $key);
    #print "calc_rasmol: number of residues = $#count_array\n";

    unless ($#count_array > 0)
    {
	print("Too few residues in ql file!\n");
	return;
    }
    
    #my $script_fname = &create_rasmol_script_filename($ql_file);
    #print "calc_rasmol: creating script_fname = $script_fname\n";

#    my $template = "rasmol_script";
#    my $script_fname = getTempFileName($template);
    my $script_fname = "precision_script.pml";
    my $script = &create_pymol_script(\@count_array, $pdb_file );
    my $write_success = &write_script_file($script_fname, $script);
    
}


sub create_pymol_script
{
    my ($count_array_ref, $pdb_file ) = @_;
    my @count_array = @$count_array_ref;
    
    my $script = '';
    my $label_s = '';	
    my $color;
    my ($i, $j);

    #options for rasmol
    my $image_mode = "cartoon";    #trace or cartoons

    my $RED = 370;
    my $ORANGE = 268;
    my $YELLOW = 165;
    my $GREEN = 121;
    my $BLUE_GREEN = 78;
    my $BLUE = 35; 
    
    $script = $script."load $pdb_file\n";
    $script = $script."hide everything, all\ndss\nset stick_transparency,0.6\nset cartoon_transparency,0.5\nshow $image_mode, all\nset bg_rgb, [1,1,1]\n";
    $script = $script."color blue, all\n";
    $script = $script."set label_color,black\nset label_font_id, 4\nset label_size,3\n";
    #loop through count_array
    for ($i = 0; $i<=$#count_array; $i++)
    {
	#print "count_array[$i] = $count_array[$i]\n";
	next unless ($count_array[$i]);

	my $count = $count_array[$i];
	$label_s .= "label ${i}/ca,\"%s\" %(resi)\n" if($count >= $YELLOW);
	#default is blue
	next unless ($count > $BLUE);
	
	if ($count >= $RED)
	{
	    $color = "[255, 0, 0]";
	    $script = $script.&create_color_script($color, $i, 1);
	    next;
	}
	elsif ($count >= $ORANGE)
	{
	    my $green_component = (($RED - $count) * 128)/($RED - $ORANGE);
	    $green_component = sprintf("%.0f", $green_component);
	    $color = "[255, $green_component, 0]";
	    $script = $script.&create_color_script($color, $i, 1);
	}
	elsif ($count >= $YELLOW)
	{
	    my $green_component =  127 + (($ORANGE - $count) * 128)/($ORANGE - $YELLOW);
	    $green_component = sprintf("%.0f", $green_component);
	    $color = "[255, $green_component, 0]";
	    $script = $script.&create_color_script($color, $i, 1);
	}
	elsif ($count >= $GREEN)
	{
	    my $red_component = (($count - $GREEN) * 255)/($YELLOW - $GREEN);
	    $red_component = sprintf("%.0f", $red_component);
	    $color = "[$red_component, 255, 0]";
	    $script = $script.&create_color_script($color, $i, 0);
	}
	elsif ($count >= $BLUE_GREEN)
	{
	    my $blue_component = (($GREEN - $count) * 255)/($GREEN - $BLUE_GREEN);
	    $blue_component = sprintf("%.0f", $blue_component);
	    $color = "[0, 255, $blue_component]";
	    $script = $script.&create_color_script($color, $i, 0);	   
	}
	elsif ($count > $BLUE)
	{
	    my $green_component = (($count - $BLUE) * 255)/($BLUE_GREEN - $BLUE);
	    $green_component = sprintf("%.0f", $green_component);
	    $color = "[0, $green_component, 255]";
	    $script = $script.&create_color_script($color, $i, 0);
	}

    }
    $script = $script . $label_s;	 
  $script .= "orient\n";
    $script = $script."ray 800,800\npng precision.png\nquit\n";
    return($script);
}


sub create_color_script
{
    my ($color, $residue,$extra) = @_;
    
    my $script = "set_color color${residue}, $color\n";
    $script = $script."color color${residue}, resi $residue\n";
    if($extra){
	$script .= "show sticks, res $residue\n";
    }	
    return ($script);
}

sub get_count_array
{
    my ($file, $filter, $key) = @_;
#    my $basepath = Cwd::cwd;

    unless (-e $file)
    {
	
	return;
    }
    
    my ($data_ref, $heading_ref, $residue_ref) = &parse_ql($file);
    my @heading = @$heading_ref;
    my @residue = @$residue_ref;

    my @data = @{$data_ref};
    
    #array of counts for each residue
    my @count;
    my $dummy = 0;
    my $section = 1;
    my ($i, $j, $first_atom, $second_atom);

    #print "\n\ncalc_count: number of residues = $#residue\n";

    for ($i=0; $i<=$#residue; $i++)
    {
	next unless $residue[$i];
	#print "Residue[$i] = $residue[$i]\n";
	
	$count[$i] = &getCount(\@data, $i, $residue_ref, $heading_ref, $filter, $key);
	#print "get_count_array: Residue $i - $residue[$i], count is $count[$i]\n";
    }

    return (@count);
}


#given i, returns the count of residue i, exluding intra residue contacts
sub getCount
{
    my ($data_ref, $i, $residue_ref, $heading_ref, $filter, $key) = @_;
    
    my @data = @$data_ref;
    my @residue = @$residue_ref;
    my @heading = @$heading_ref;
    
    my %dist_hash;

    my $dist = 0;

    my $count = 0;
    #print "\n******getCount: getting count of residue $i*******\n";

    my ($j, $i_atom, $j_atom, $k, $dist_hash_ref);
    $dist_hash_ref = \%dist_hash;
    
    #print "**Residue $i on left\n";
    #Left residue = $i, Right residue = $j - loop

    if ($key =~ /summary/i) #for the overal summary
    {
	for ($k = 0; $k<=$#heading; $k++)
	{
	    #next if ($k == 1);
	    #next if ($k == 3);
	    
	    &getLeftCountofi($data_ref, $i, $k, $filter, $dist_hash_ref, $residue_ref, 0); # 0 - no intra 
	}
	for ($k = 0; $k<=$#heading; $k++)
	{
	    #next if ($k == 1);
	    #next if ($k == 3);
	    
	    &getRightCountofi($data_ref, $i, $k, $filter, $dist_hash_ref, $residue_ref, 0);
	}
	
    }
    else
    {
	for (my $m = 0; $m<=$#heading; $m++) # for each peak list 
	{
	    if ($heading[$m] =~ /$key/)
	    {
		$k = $m;
		last;
	    }
	}

	&getLeftCountofi($data_ref, $i, $k, $filter, $dist_hash_ref, $residue_ref, 0);
	&getRightCountofi($data_ref, $i, $k, $filter, $dist_hash_ref, $residue_ref, 0);
    }


    #sum up all distances from dist_hash
    my $key;
    my $count = 0;
    
    #clean-up by removing the non-symmetric ones 
    my %nonsym; 
    foreach $key (keys %$dist_hash_ref)
    {
	my ($a1, $a2) = split '\s', $key; 

	if ($dist_hash_ref->{$a2 . ' ' . $a1} != $dist_hash_ref->{$key})
	{
	    $nonsym{$key} = 1; 
	} 
    }

    foreach $key (keys %$dist_hash_ref)
    {
	next if ($nonsym{$key} == 1);  

	my $dist = $dist_hash_ref->{$key};
	#$dist = 30000 * ($dist ** (-6))/2;
	$dist = 15000 * ($dist ** (-6))/2; #devided by 2 - remove the symmetry one 
	$count = $count + $dist;
    }

    #print "total distance for residue $i is $dist\n";
    #print "getCount: returning count of $i as $count\n";
    return $count;

}

#filter is used in calc_rasmol analysis
#getRightCountofi (report x-i (x<i)
#getLeftCountofi (report i-x (x>=i)) 
#i-x (i<=x)
sub getLeftCountofi
{
    my ($data_ref, $i, $k, $filter, $dist_hash_ref, $residue_ref, $i_atom_filter, $j_filter, $j_atom_filter, $distance_comp_filter, $distance_filter, $intra) = @_;
    
    #print "getLeftCountofi: i = $i, k = $k, i_atom_filter = $i_atom_filter, j_filter = $j_filter, j_atom_filter = $j_atom_fitler, distance_comp_filter = $distance_comp_filter, distance = $distance_filter, filter = $filter\n";

    
    my @residue = @$residue_ref;
    my @data = @$data_ref;

    my ($i_atom, $j, $j_atom);

    foreach $i_atom (keys %{$data[$k][$i]})
    {
	if ($i_atom =~ /^$filter$/i)
	{
	    #print "filtered i: $i-$i_atom\n";
	    next;
	}
	
	for ($j = 0; $j<scalar(@{$data[$k][$i]{$i_atom}}); $j++)
	{
	    next unless $residue[$j];
	    next if (($i == $j) && !$intra);  
	    #next if (($i == $j)); 
		
	    if ($j_filter)
	    {
		next unless ($j == $j_filter);
	    }

	    foreach $j_atom (keys %{$data[$k][$i]{$i_atom}[$j]})
	    {
				
		if ($j_atom =~ /^$filter$/i)
		{
		    #print "filtered j: $j-$j_atom\n";
		    next;
		}

		#atom filter		    
		if ($j == $i)  #intra -  
		{
		    if ($j_atom_filter eq '')
		    {
			$j_atom_filter = '\w'; 
		    }
		    if ($i_atom_filter eq '')
		    {
			$i_atom_filter = '\w'; 
		    }
		    
		    next unless ((($i_atom =~ /$i_atom_filter/i) && ($j_atom=~ /$j_atom_filter/i)) ||(($i_atom =~ /$j_atom_filter/i) && ($j_atom=~ /$i_atom_filter/i)));  		    
		} else
		{
		    if (($j_atom_filter)) # j atom only 
		    {
			next unless (($j_atom =~ /$j_atom_filter/i));
		    } 
		    if (($i_atom_filter)) # i atom only 
		    {
			next unless (($i_atom =~ /$i_atom_filter/i));
		    }
		}
		

		if ($distance_filter)
		{
		    my $distance = $data[$k][$i]{$i_atom}[$j]{$j_atom};
		    #print "getLeftCountofi: distance = $distance, only distances $distance_comp_filter than $distance_filter.  ";

		    if ($distance_comp_filter eq '=')
		    {
			$distance_comp_filter = '==';
			#print "distance_comp_filter is now $distance_comp_filter\n";
		    }
		    unless (eval("$distance $distance_comp_filter $distance_filter"))
		    {
			#print "distance $distance filtered\n";
			next;
		    }
		}
		
		#print "residue $i on left - distance from $i-$i_atom to $j-$j_atom is $data[$i]{$i_atom}[$j]{$j_atom}\n";
		$dist_hash_ref->{"$i-$i_atom $j-$j_atom"} = $data[$k][$i]{$i_atom}[$j]{$j_atom};
		#print "getLeftCountofi: dist_hash_ref->{$i-$i_atom $j-$j_atom} = ".$dist_hash_ref->{"$i-$i_atom $j-$j_atom"}."\n";
	    }
	}
    }
}


#x-i, x < i 
sub getRightCountofi
{

    my ($data_ref, $i, $k, $filter, $dist_hash_ref, $residue_ref, $i_atom_filter, $j_filter, $j_atom_filter, $distance_comp_filter, $distance_filter, $intra) = @_;

    #print "getRightCountofi: i = $i, k = $k, i_atom_filter = $i_atom_filter, j_filter = $j_filter, j_atom_filter = $j_atom_fitler, distance_comp_filter = $distance_comp_filter, distance = $distance_filter, filter = $filter\n";

    my @residue = @$residue_ref;
    my @data = @$data_ref;

    my ($i_atom, $j_atom);

    #print "**Residue $i on right\n";
    
    #right residue = i, left residue = j
    for (my $j = 0; $j<scalar(@{$data[$k]}); $j++)
    {
	next unless $residue[$j];
	#next if ($i == $j);
	next if (($i == $j) && !$intra);#remove repeat entry from getRightCountofI

	if ($j_filter)
	{
	    next unless ($j == $j_filter);
	}
	
	foreach $j_atom (keys %{$data[$k][$j]})
	{
	    if ($j_atom_filter)
	    {
		next unless ($j_atom =~ /$j_atom_filter/i);
	    }
	    
	    if ($j_atom =~ /^$filter$/i)
	    {
		#print "filtered j: $j-$j_atom\n";
		next;
	    }
	    
	    foreach $i_atom (keys %{$data[$k][$j]{$j_atom}[$i]})
	    {
		if ($i_atom_filter)
		{
		    next unless ($i_atom =~ /$i_atom_filter/i);
		}
		
		if ($i_atom =~ /^$filter$/i)
		{
		    #print "filtered i: $i-$i_atom\n";
		    next;
		}
		
		if ($distance_filter)
		{
		    my $distance = $data[$k][$j]{$j_atom}[$i]{$i_atom};
		    #print "getRightCountofi: distance $distance, only distances $distance_comp_filter than $distance_filter.  ";
		    if ($distance_comp_filter eq '=')
		    {
			$distance_comp_filter = '==';
			#print "distance_comp_filter is now $distance_comp_filter\n";
		    }

		    unless (eval("$distance $distance_comp_filter $distance_filter"))
		    {
			#print "distance $distance is filtered\n";
			next;
		    }

		    #next unless (eval("$distance $distance_comp_filter $distance_filter"));
		}

		#print "residue $i on right - distance from $j-$j_atom to $i-$i_atom is $data[$j]{$j_atom}[$i]{$i_atom}\n";
		unless ($dist_hash_ref->{"$j-$j_atom $i-$i_atom"})
		{
		    #$dist_hash_ref->{"$i-$i_atom $j-$j_atom"} = $data[$k][$j]{$j_atom}[$i]{$i_atom};
		    $dist_hash_ref->{"$j-$j_atom $i-$i_atom"} = $data[$k][$j]{$j_atom}[$i]{$i_atom};
		    #print "getRightCountofi: dist_hash_ref->{$i-$i_atom $j-$j_atom} = ".$dist_hash_ref->{"$i-$i_atom $j-$j_atom"}."\n";
		}
	    }
	}
    }
}

sub write_script_file
{
    my ($scr_fname, $script) = @_;
    
    unless (open(FILE, ">$scr_fname"))
    {
	print("Cannot open Rasmol script file \"$scr_fname\" for writing\n");
	return 0;
    }
    
    print FILE $script;
    close(FILE);
    return 1;
}

sub parse_ql
{
    my $file = shift;
 
    #two dimensional array with each element as a hash having three keys - distance, residue number, and atom name
    
    my @data;
    
    #my $data_count = 0;
    my $first_line_seen;
    my $pair_seen;
    my $section = -1;
    my @section_heading;
    my @atom;
    my $count = 0;

    my $other_first_line;

    my ($first, $first_name, $second, $second_name, $distance);
    my ($i, $j, $option);

    	    
    unless (-e $file)
    {
	print "No QL file specified\n";
	return;
    }

    my $opened = open(QL, $file);

    unless ($opened)
    {
	print "Cannot open QL file \"$file\"";
	return;
    }
    
    
    #read until first line
    my $line;	
    while (chomp($line = <QL>))
    {
	if ($line =~ /missing/i)
	{
	    $line =~ /missing .+ for (.+)$/i;
	    
	    $section++;
	    $section_heading[$section] = $1;

	    #$data_count = 0;
	    $first_line_seen = 0;
	    $other_first_line = 1;
	    $pair_seen = 0;
	}
	elsif ($line =~ /^atom/i)
	{
	    #only look at first pair
	    if (!$pair_seen)
	    {
		if (!$first_line_seen)
		{
		    #first line
		    
		    $line =~ /^atom:\s*\d+\s*(\d+)\w-(\w+)\s*.*/;
		    $first = $1;
		    $first_name = $2;
		    $first_line_seen = 1;
		    
		}
		else
		{
		    #second line
		    
		    $line =~ /^atom:\s*\d+\s*(\d+)\w-(\w+)\s*.*/;
		    $second = $1;
		    $second_name = $2;
		    $pair_seen = 1;
		}
		
	    }
	    else
	    {
		#other lines
		if ($other_first_line)
		{
		    $line =~ /^atom:\s*\d+\s*(\d+)\w-(\w+)\s*.*/;
		    my $amb_atom = $2;
		    
		    unless ($amb_atom eq $first_name)
		    {
			$amb_atom =~ /(\w+)(\d)$/;
			my $amb_num = $2;
			my $amb_atm = $1;
			$first_name =~ /(\w+)(\d)$/;
			my $fir_num = $2;
			my $fir_atm = $1;

			if ($amb_num and $fir_num)
			{
			    if (($amb_atm eq $fir_atm) and ($amb_num != $fir_num))
			    {
				$first_name = $amb_atm."*";
			    }
			}
		    }
		}
		else
		{
		    $line =~ /^atom:\s*\d+\s*(\d+)\w-(\w+)\s*.*/;
		    my $amb_atom = $2;
		    
		    unless ($amb_atom eq $second_name)
		    {
			$amb_atom =~ /(\w+)(\d)$/;
			my $amb_num = $2;
			my $amb_atm = $1;
			$second_name =~ /(\w+)(\d)$/;
			my $sec_num = $2;
			my $sec_atm = $1;

			if ($amb_num and $sec_num)
			{
			    if (($amb_atm eq $sec_atm) and ($amb_num != $sec_num))
			    {
				$second_name = $amb_atm."*";
			    }
			}
		    }
		}
		$other_first_line = !$other_first_line;
		
	    }
	    
	    #store atom and number in number_to_atom array
	    $line =~ /^atom:\s*\d+\s+(\d+)([a-zA-Z])-.+/;
	    $atom[$1] = $2;
	}
	elsif ($line =~ /^::sumavg:\s*(.+)$/i)
	{
	    $distance = $1;
	    
	    $data[$section][$first]{$first_name}[$second]{$second_name} = $distance;
	    #print "parse_ql: data[$section][$first]{$first_name}[$second]{$second_name} = $distance\n";
	    

	    #print "data[$section][$ql_data{'1-residue number'}]{$ql_data{'1-atom name'}}[$ql_data{'2-residue number'}]{$ql_data{'2-atom name'}} = $ql_data{'distance'}\n";
	    #print "data[$section][$ql_data{'2-residue number'}]{$ql_data{'2-atom name'}}[$ql_data{'1-residue number'}]{$ql_data{'1-atom name'}} = $ql_data{'distance'}\n";

	    
	    #$data_count++;
	    $first_line_seen = 0;
	    $other_first_line = 1;
	    $pair_seen = 0;

	    
	}	
    }
    return (\@data, \@section_heading, \@atom);

}