diff --git a/scripts/dist_mutation_to_na2.pl b/scripts/dist_mutation_to_na2.pl new file mode 100644 index 0000000..49bbd81 --- /dev/null +++ b/scripts/dist_mutation_to_na2.pl @@ -0,0 +1,257 @@ +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 +___________________________________________________________________________________\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 = ; +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"; +} +#----------------------------------------------------------------------------------------