I have made some progress, but I cant figure it out why script doesnt work. I have wrote two loops and pushed values of coordinates into seperate arrays and ttried to print distances but it doesnt work. Can You please take a look cause I am trying to change it numerous times. This is the script :
#!/usr/local/bin/perl
use strict;
use warnings;
open(IN, $ARGV[0]) or die "$!";
my (@refer, @points);
my $part = 0;
my $dist;
while (my $line = <IN>) {
chomp($line);
if ($line =~ /^HETATM/) {
$part++;
next;
}
my @array = (substr($line, 30, 8),substr($line,38,8),substr($line,46,8));
### MD: use split to split by delimiter:
my @coords = split $line;
### MD: not sure what this is for. Do you need to distinguish between
### two types of atoms? I thought you wanted just
### the distances between all atoms?
### why not generate a simple nested list of all atom coords first?
# print "@array\n";
if ($part == 0) {
push @refer, [ @array ];
} elsif ($part ==1){
push @points, [ @array ];
}
}
### for calculation of a symmetric distance a nested for loop is more efficient
### than a simple foreach, because you can calculate the the lower triangle
### something like this should do the trick, not tested, and you need to define a sub dist
my $distMat = [];
for (my $i=0; $i < $M; $i++) {
$distMat->[$i] = [];
for (my $j=$i; $j < $N; $j++) {
$distMat->[$i][$j] = dist($coords[$i], $coords[$j]);
}
}
foreach my $ref(@refer) {
my ($x1, $y1, $z1) = @{$ref};
foreach my $atom(@points) {
my ($x, $y, $z) = @{$atom};
my $dist = sqrt( ($x-$x1)**2 + ($y-$y1)**2 + ($z-$z1)**2 );
print $dist;
}
}
• 0 views
•
link
updated
by
Michael Dondrup
•
written
by
srdjanmasirevic2 •
•
written
by
srdjanmasirevic2 •
This looks like homework, what have you tried?
It's some kind of personal homework, not for school but form me to make some dataset. So what I have done so far is wrote a script in perl that from the whole pdb dataset I kept only those pdb files that have ligand in its binding pocket (I included some cut off for number of heavy atoms to exclude ions and solvents, so ended up with pdb files that have ligands). I also extracted from each file protein and ligand coordinates into another txt file. So currently I am having a files that look like this: for example file. 2KRJ.txt:
The values for these coordinates are not the same of course. So I also found the formula fo calculating distances, but I am not sure how to calculate it: I want each atom (x,y,z coordinates) from protein (that starts with
^ATOM) to be calculated with each atom from ligand (that starts with^HETATM). If the number of each calculatation is more then 5 (I should assigned some cut off) then move the file or something like that.This is the formula: distance =
sqrt(($Ax - $Bx)**2 + ($Ay - $By)**2 + ($Az - $Bz)**2)I am not that experienced in programming so trying to figure out how to approach to this problem. Thanks
So basically I am trying to calculate
ATOM(1)with allHETATM(1,2,3,4)and the same for each atom of the protein. if I end up with result that is <= to five then I should keep this file. Overall from all files I am trying to keep only those in which the distance between ligand and some protein atoms is less or equal to 5 Armstrongs.Ok, that is totally feasible in Perl with two nested loops or maps, but why don't you use R on the file you made? read the file using
atoms <- read.delim(my.file, sep="\t", header=T, row.names=1) ; dist(atoms);that what be about 1/10 of the code required in Perl.