LSHTM_analysis/scripts/dist_mutation_to_na2.pl

use strict;
use warnings;

sub trim;
sub distance;
sub res_cod1_to_res_cod3;
sub res_cod3_to_res_cod1;

# ____________________________________________________________________________________________________________________
# Input parameters
my $pdb = $ARGV[0];
my $mutation = $ARGV[1];
my $wt_chain = $ARGV[2];

if(scalar(@ARGV) != 3){
print "___________________________________________________________________________________
SINTAX:
	perl dist_mutation_to_na.pl <pdb> <mutation> <chain>
___________________________________________________________________________________\n";
	exit;
}

# ____________________________________________________________________________________________________________________

my $wild_res = substr($mutation, 0, 1);
my $wild_res_pos = substr($mutation, 1, length($mutation)-2);
my $mutated_res = substr($mutation, length($mutation)-1, 1);

open(PDB,"<$pdb") or die "$!Erro ao abrir: $pdb\n";
my @pdb = <PDB>;
close PDB;
	

my $k = 0;   
my @coord_x;            my @coord_y;        
my @coord_z;            my @res_num;
my @res_name;   	my @min_dist;        
my @chain;

# ==================================================================================================
foreach my $line (@pdb){
	if($line =~ /^ATOM|^HETATM/){
		if(	trim(substr($line,17,3)) eq "DA" or
			trim(substr($line,17,3)) eq "DG" or
			trim(substr($line,17,3)) eq "DC" or
			trim(substr($line,17,3)) eq "DT" or
			
			trim(substr($line,17,3)) eq "A" or
			trim(substr($line,17,3)) eq "G" or
			trim(substr($line,17,3)) eq "C" or
			trim(substr($line,17,3)) eq "U"
		){
			my $res_cod = res_cod3_to_res_cod1(trim(substr($line,17,3)));
			my $res_ind = trim(substr($line,22,4));
			my $x = trim(substr($line,30,8));
			my $y = trim(substr($line,38,8));
			my $z = trim(substr($line,46,8));
	
			$coord_x[$k] = $x;
			$coord_y[$k] = $y;
			$coord_z[$k] = $z;
	
			$res_num[$k] = $res_ind;
			$res_name[$k] = $res_cod;
		
			$chain[$k] = substr($line,21,1);
	
			$k++;
		}
	}
}

my $k2 = 0;
my @coord_x2;            my @coord_y2;
my @coord_z2;            my @res_num2;
my @res_name2;           my @min_dist2;
my @chain2;

foreach my $line (@pdb){
	if(trim(substr($line,0,6)) eq "ATOM"){
		my $res_cod = res_cod3_to_res_cod1(trim(substr($line,17,3)));
		my $res_ind = trim(substr($line,22,4));
		my $x = trim(substr($line,30,8));
		my $y = trim(substr($line,38,8));
		my $z = trim(substr($line,46,8));
		my $curr_chain = substr($line,21,1);

		if($wild_res_pos == $res_ind and $wt_chain eq $curr_chain){

			$coord_x2[$k2] = $x;
			$coord_y2[$k2] = $y;
			$coord_z2[$k2] = $z;

			$res_num2[$k2] = $res_ind;
			$res_name2[$k2] = $res_cod;

			$chain2[$k2] = substr($line,21,1);

			$k2++;
		}
	}
}

#print "$k2\t$k\n";
#print "Calculating distances\n";

# ==================================================================================================

my $min_dist = 999;
for(my $i=0; $i<$k2; $i++){
	for(my $j=0; $j<$k; $j++){

		my $dist = distance($coord_x2[$i],$coord_y2[$i],$coord_z2[$i],$coord_x[$j],$coord_y[$j],$coord_z[$j]);

		my $res_ind1 = $res_num2[$i];
		my $res_ind2 = $res_num[$j];

		if($min_dist > $dist){
			$min_dist = $dist;
		}
	}
}

printf '%.3f'."\n", $min_dist;
exit;
# ____________________________________________________________________________________________________________________
sub trim{
	my $string = shift;
	$string =~ s/^\s+//;
	$string =~ s/\s+$//;
	return $string;
}

sub distance{
	my($x1,$y1,$z1,$x2,$y2,$z2) = @_;
	my $distance;
	
	$distance = sqrt(($x1-$x2)**2 + ($y1-$y2)**2 + ($z1-$z2)**2);
	return $distance;
}

sub res_cod3_to_res_cod1{
	my $cod3 = shift;
		
	if($cod3 eq "ALA"){
		return "A";
	} elsif($cod3 eq "VAL"){
		return "V";
	} elsif($cod3 eq "LEU"){
		return "L";
	} elsif($cod3 eq "GLY"){
		return "G";
	} elsif($cod3 eq "SER"){
		return "S";
	} elsif($cod3 eq "TRP"){
		return "W";
	} elsif($cod3 eq "THR"){
		return "T";
	} elsif($cod3 eq "GLN"){
		return "Q";
	} elsif($cod3 eq "GLU"){
		return "E";
	} elsif($cod3 eq "CYS"){
		return "C";
	} elsif($cod3 eq "ARG"){
		return "R";
	} elsif($cod3 eq "PRO"){
		return "P";
	} elsif($cod3 eq "ASP"){
		return "D";
	} elsif($cod3 eq "PHE"){
		return "F";
	} elsif($cod3 eq "ILE"){
		return "I";
	} elsif($cod3 eq "HIS"){
		return "H";
	} elsif($cod3 eq "ASN"){
		return "N";
	} elsif($cod3 eq "MET"){
		return "M";
	} elsif($cod3 eq "TYR"){
		return "Y";
	} elsif($cod3 eq "LYS"){
		return "K";
	}
	return "ERRO";
}

#----------------------------------------------------------------------------------------
# Recebe codigo de residuo de um caractere e retorna o equivalente de tres
sub res_cod1_to_res_cod3($){
	
	my $cod1 = shift;
		
	if($cod1 eq "A"){
		return "ALA";
	}
	elsif($cod1 eq "V"){
		return "VAL";
	}
	elsif($cod1 eq "L"){
		return "LEU";
	}
	elsif($cod1 eq "G"){
		return "GLY";
	}
	elsif($cod1 eq "S"){
		return "SER";
	}
	elsif($cod1 eq "W"){
		return "TRP";
	}
	elsif($cod1 eq "T"){
		return "THR";
	}
	elsif($cod1 eq "Q"){
		return "GLN";
	}
	elsif($cod1 eq "E"){
		return "GLU";
	}
	elsif($cod1 eq "C"){
		return "CYS";
	}
	elsif($cod1 eq "R"){
		return "ARG";
	}
	elsif($cod1 eq "P"){
		return "PRO";
	}
	elsif($cod1 eq "D"){
		return "ASP";
	}
	elsif($cod1 eq "F"){
		return "PHE";
	}
	elsif($cod1 eq "I"){
		return "ILE";
	}
	elsif($cod1 eq "H"){
		return "HIS";
	}
	elsif($cod1 eq "N"){
		return "ASN";
	}
	elsif($cod1 eq "M"){
		return "MET";
	}
	elsif($cod1 eq "Y"){
		return "TYR";
	}
	elsif($cod1 eq "K"){
		return "LYS";
	}
	return "ERRO";
}
#----------------------------------------------------------------------------------------