Расчет расстояния между атомными координатами
У меня есть текстовый файл, как показано ниже
ATOM 920 CA GLN A 203 39.292 -13.354 17.416 1.00 55.76 C
ATOM 929 CA HIS A 204 38.546 -15.963 14.792 1.00 29.53 C
ATOM 939 CA ASN A 205 39.443 -17.018 11.206 1.00 54.49 C
ATOM 947 CA GLU A 206 41.454 -13.901 10.155 1.00 26.32 C
ATOM 956 CA VAL A 207 43.664 -14.041 13.279 1.00 40.65 C
.
.
.
ATOM 963 CA GLU A 208 45.403 -17.443 13.188 1.00 40.25 C
Я хотел бы рассчитать расстояние между двумя альфа-атомами углерода, то есть рассчитать расстояние между первым и вторым атомом, а затем между вторым и третьим атомом и т. Д...... Расстояние между двумя атомами может быть выражено как:distance = sqrt((x1-x2)^2+(y1-y2)^2+(z1-z2)^2) .
Столбцы 7,8 и 9 представляют координаты x,y и z соответственно. Мне нужно напечатать расстояние и соответствующие пары остатков (столбец 4), как показано ниже.(Значения расстояния не являются реальными)
GLN-HIS 4.5
HIS-ASN 3.2
ASN-GLU 2.5
Как я могу сделать этот расчет с Perl или Python?
6 ответов
Если ваши данные разделены пробелами, простой split
могу сделать работу. Буферизуем строки, чтобы сравнить их друг с другом последовательно.
use strict;
use warnings;
my @line;
while (<>) {
push @line, $_; # add line to buffer
next if @line < 2; # skip unless buffer is full
print proc(@line), "\n"; # process and print
shift @line; # remove used line
}
sub proc {
my @a = split ' ', shift; # line 1
my @b = split ' ', shift; # line 2
my $x = ($a[6]-$b[6]); # calculate the diffs
my $y = ($a[7]-$b[7]);
my $z = ($a[8]-$b[8]);
my $dist = sprintf "%.1f", # format the number
sqrt($x**2+$y**2+$z**2); # do the calculation
return "$a[3]-$b[3]\t$dist"; # return the string for printing
}
Вывод (с образцами данных):
GLN-HIS 3.8
HIS-ASN 3.8
ASN-GLU 3.9
GLU-VAL 3.8
Если ваши данные разделены табуляцией, вы можете разделить на /\t/
вместо ' '
,
НЕ разбивать на пробелы
Другие ответы, приведенные здесь, делают ошибочное предположение - что координаты будут разделены пробелом. Согласно спецификации PDBATOM
Это не обязательно во всех случаях: значения записей PDB задаются индексами столбцов и могут переходить друг в друга. Например, ваш первый ATOM
запись гласит:
ATOM 920 CA GLN A 203 39.292 -13.354 17.416 1.00 55.76 C
Но это также совершенно справедливо:
ATOM 920 CA GLN A 203 39.292-13.3540 17.416 1.00 55.76 C
Лучший подход
Из-за указанных в столбцах индексов и ряда других проблем, которые могут возникнуть в файле PDB, вы не должны писать свой собственный анализатор. Формат PDB грязный, и есть много особых случаев и плохо отформатированных файлов для обработки. Вместо этого используйте парсер, который уже написан для вас.
Мне нравятся биопионыPDB.PDBParser
, Он разберет структуру для вас в объекты Python, в комплекте с удобными функциями. Если вы предпочитаете Perl, проверьте BioPerl.
PDB.Residue
объекты разрешают доступ к атомам по ключу по имени, и PDB.Atom
объекты перегружают -
оператор, чтобы вернуть расстояние между двумя атомами. Мы можем использовать это, чтобы написать чистый, краткий код:
Код
from Bio import PDB
parser = PDB.PDBParser()
# Parse the structure into a PDB.Structure object
pdb_code = "1exm"
pdb_path = "pdb1exm.ent"
struct = parser.get_structure(pdb_code, pdb_path)
# Grab the first two residues from the structure
residues = struct.get_residues()
res_one = residues.next()
res_two = residues.next()
try:
alpha_dist = res_one['CA'] - res_two['CA']
except KeyError:
print "Alpha carbon missing, computing distance impossible!"
Предполагая, что ваши данные находятся в "atom.txt", он читает их построчно и разбивает записи в список:
import csv
with open("atoms.txt") as f:
reader = csv.reader(f)
for line, in reader:
entries = line.split()
print entries
Теперь для каждого списка извлеките нужные вам столбцы, вычислите расстояния и т. Д. (Имейте в виду, что списки в python начинаются с нуля).
В идеале вы должны использовать пакет MDAnalysis для питонического способа "разрезания" атомов и сегментов и вычисления расстояний между ними. Фактически, MDAnalysis поддерживает несколько симуляций MD и форматы химической структуры.
Для более подробного примера см. Также следующую запись на Biostars.org.
Если вас интересует только одна пара, bash работает просто отлично. Это скрипт, который я использую, он настроен на повторный запуск в конце (отключите его, если хотите). Он подскажет вам, для какого атома. Файлы PDB могут иметь разные настроенные столбцы, поэтому для строки awk убедитесь, что столбцы совпадают. Выполните тестовый пример вручную перед использованием с новым файлом pdb. Это тривиально, но измените в скрипте мой файл pdb на ваш.
#!/usr/bin/env bash
echo " "
echo "Use one letter abbreviations. Case doesn't matter."
echo "Example: A 17 CA or n 162 cg"
echo " - - - - - - - - - - - - - - - - - -"
#------------First Selection------------
read -e -p "What first atom? " sel1
# echo $sel1
sel1caps=${sel1^^}
# echo "sel1caps="$sel1caps
arr1=($sel1caps)
# echo "arr1[0]= "${arr1[0]}
# echo "arr1[1]= "${arr1[1]}
# echo "arr1[2]= "${arr1[2]}
#To convert one to three letters
if [ ${arr1[0]} = A ] ; then
SF1=ALA
elif [ ${arr1[0]} = H ] ; then
SF1=HIS
elif [ ${arr1[0]} = R ] ; then
SF1=ARG
elif [ ${arr1[0]} = K ] ; then
SF1=LYS
elif [ ${arr1[0]} = I ] ; then
SF1=ILE
elif [ ${arr1[0]} = F ] ; then
SF1=PHE
elif [ ${arr1[0]} = L ] ; then
SF1=LEU
elif [ ${arr1[0]} = W ] ; then
SF1=TRP
elif [ ${arr1[0]} = M ] ; then
SF1=MET
elif [ ${arr1[0]} = P ] ; then
SF1=PRO
elif [ ${arr1[0]} = C ] ; then
SF1=CYS
elif [ ${arr1[0]} = N ] ; then
SF1=ASN
elif [ ${arr1[0]} = V ] ; then
SF1=VAL
elif [ ${arr1[0]} = G ] ; then
SF1=GLY
elif [ ${arr1[0]} = S ] ; then
SF1=SER
elif [ ${arr1[0]} = Q ] ; then
SF1=GLN
elif [ ${arr1[0]} = Y ] ; then
SF1=TYR
elif [ ${arr1[0]} = D ] ; then
SF1=ASP
elif [ ${arr1[0]} = E ] ; then
SF1=GLU
elif [ ${arr1[0]} = T ] ; then
SF1=THR
else
echo "use one letter codes"
echo "exiting"
exit
fi
# echo "SF1 ="$SF1
#If there is nothing printing for line 1, check the expression for your pdb file. The spaces may need adjustment at the end.
line1=$(grep -E "${arr1[2]} *?${SF1}(.*?) ${arr1[1]} " 1A23.pdb)
# echo $line1
ar_l1=($line1)
# echo "ar_l1[1]="${ar_l1[1]}
echo " - - - - - - - - - - - - - - - - - -"
#------------Second Selection------------
read -e -p "What second atom? " sel2
# echo $sel2
sel2caps=${sel2^^}
# echo "sel2caps="$sel2caps
arr2=($sel2caps)
# echo "arr2[0]= "${arr2[0]}
# echo "arr2[1]= "${arr2[1]}
# echo "arr2[2]= "${arr2[2]}
#To convert one to three letters
if [ ${arr2[0]} = A ] ; then
SF2=ALA
elif [ ${arr2[0]} = H ] ; then
SF2=HIS
elif [ ${arr2[0]} = R ] ; then
SF2=ARG
elif [ ${arr2[0]} = K ] ; then
SF2=LYS
elif [ ${arr2[0]} = I ] ; then
SF2=ILE
elif [ ${arr2[0]} = F ] ; then
SF2=PHE
elif [ ${arr2[0]} = L ] ; then
SF2=LEU
elif [ ${arr2[0]} = W ] ; then
SF2=TRP
elif [ ${arr2[0]} = M ] ; then
SF2=MET
elif [ ${arr2[0]} = P ] ; then
SF2=PRO
elif [ ${arr2[0]} = C ] ; then
SF2=CYS
elif [ ${arr2[0]} = N ] ; then
SF2=ASN
elif [ ${arr2[0]} = V ] ; then
SF2=VAL
elif [ ${arr2[0]} = G ] ; then
SF2=GLY
elif [ ${arr2[0]} = S ] ; then
SF2=SER
elif [ ${arr2[0]} = Q ] ; then
SF2=GLN
elif [ ${arr2[0]} = Y ] ; then
SF2=TYR
elif [ ${arr2[0]} = D ] ; then
SF2=ASP
elif [ ${arr2[0]} = E ] ; then
SF2=GLU
elif [ ${arr2[0]} = T ] ; then
SF2=THR
else
echo "use one letter codes"
echo "exiting"
exit
fi
# echo "SF2 ="$SF2
line2=$(grep -E "${arr2[2]} *?${SF2}(.*?) ${arr2[1]} " 1A23.pdb)
# echo $line2
ar_l2=($line2)
# echo "ar_l2[1]="${ar_l2[1]}
# echo "ar_l2[1]="${ar_l2[1]}
atom1=${ar_l1[1]}
atom2=${ar_l2[1]}
echo "==========================="
echo ${arr1[0]}${arr1[1]}${arr1[2]}" to "${arr2[0]}${arr2[1]}${arr2[2]}":"
# 6, 7, 8 are column numbers in the pdb file.
# If there are multiple molecules it should be 7, 8, 9.
awk '$2=='$atom1'{x1=$7;y1=$8;z1=$9} # get the ATOM 1
$2=='$atom2'{x2=$7;y2=$8;z2=$9} # get the ATOM 2
END{print sqrt((x1-x2)*(x1-x2) + (y1-y2)*(y1-y2) + (z1-z2)*(z1-z2))}' 1A23.pdb # calculate the distance.
echo "Angstroms"
echo "==========================="
echo " "
echo "-_-_-_-_Running Script Again_-_-_-_-"
./distance_soln.sh
Простой код Python может сделать эту работу. Я предположил, что все ваше содержимое находится в файле "input.txt".
def process(line1, line2):
content1 = line1.split()
content2 = line2.split()
x1, y1, z1 = float(content1[6]), float(content1[7]), float(content1[8])
x2, y2, z2 = float(content2[6]), float(content2[7]), float(content2[8])
distance = math.sqrt(math.pow(x1-x2, 2) + math.pow(y1-y2, 2) + math.pow(z1-z2, 2))
return content1[3]+"-"+content2[3]+" "+ "{0:.2f}".format(distance)
with open("input.txt") as f:
line1 = f.readline()
for line in f:
line2 = line
print(process(line1, line2))
line1 = line2
Пожалуйста, дайте мне знать, если вы обнаружите какие-либо расхождения или проблемы в использовании этого скрипта.