Вычисление расстояния между координатами атомов

У меня есть текстовый файл, как показано ниже

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?


person user1866262    schedule 30.11.2012    source источник
comment
У меня такая же проблема, но я работаю с python, есть ли там подобное решение?   -  person Prelude    schedule 14.02.2014


Ответы (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/ вместо ' '.

person TLP    schedule 30.11.2012
comment
Большое спасибо за ваш ответ. Ваш код работает хорошо. Но я не получаю вывод построчно. Я хотел бы изменить внешний вид моего вывода, как ваш образец вывода. Как я могу изменить? Я ценю вашу помощь!! - person user1866262; 30.11.2012
comment
Ах да, забыл... Я запускал скрипт с флагом -l во время тестирования. perl -l script.pl input.txt. Я исправлю это. - person TLP; 30.11.2012
comment
Важное примечание: координаты атома PDB не обязательно разделены пробелами (см. спецификация). - person David Cain; 11.05.2013
comment
Вместо того, чтобы вычислять расстояние от одного ЦС до следующего и так далее... как мы можем изменить приведенный выше ответ, чтобы рассчитать расстояние между одним ЦС до всех других ЦС, а затем расстояние от следующего ЦС до остальных ЦС и так далее...? - person pradeep pant; 15.12.2016

НЕ разделять на пробелы

Другие ответы, приведенные здесь, делают ошибочное предположение - координаты будут ограничены пробелами. Согласно спецификации PDB ATOM, это не обязательно дело: значения записей 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 компании Biopython. Он будет анализировать структуру для вас в объекты 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!"
person David Cain    schedule 01.12.2012

Предполагая, что ваши данные находятся в «atoms.txt», это читает их построчно и разбивает записи на список:

import csv

with open("atoms.txt") as f:
  reader = csv.reader(f)
  for line, in reader:
      entries = line.split()
      print entries

Теперь для каждого списка извлеките нужные столбцы и рассчитайте расстояния и т. д. (Помните, что списки в python отсчитываются от нуля).

person ev-br    schedule 30.11.2012
comment
Вы тестировали этот код? Не обращайте внимания на ошибку с for line, in reader: - возвратом этого будет list, у которого не будет метода split... - person Jon Clements♦; 30.11.2012
comment
Именно то, что я сказал в первый раз... Вы получите ValueError для for line, in reader, и как только это будет исправлено, вы получите AttributeError: 'list' object has no attribute 'split' для entries = line.split() - person Jon Clements♦; 30.11.2012
comment
Нет, я этого не понимаю. Копирование и вставка первых пяти строк из текстового файла в вопросе и запуск кода из моего ответа распечатывает список по желанию. (здесь python 2.6, не знаю насчет python 3 --- отличается ли csv.reader?) - person ev-br; 30.11.2012
comment
Попробуйте поместить первые пять строк в atoms.txt и запустить код точно так, как вы опубликовали... - person Jon Clements♦; 30.11.2012
comment
Как я уже говорил, это именно то, что я делаю, и это печатает список токенов. - person ev-br; 30.11.2012
comment
for line, in reader: возвращает строку, а не список. Удаление запятой (for line in reader:) возвращает список с одним элементом, следовательно, распаковка. - person ev-br; 30.11.2012
comment
Хорошо, это работает, просто пропускает пустые строки... Хм, если вы собираетесь использовать csv.reader, то укажите delimiter=' ', и вам не нужно его разбивать, или просто не используйте ридер...? - person Jon Clements♦; 30.11.2012
comment
А, хорошо, значит, все-таки работает :-). Что касается пустых строк --- я не вижу упоминания о пустых строках в ОП. - person ev-br; 30.11.2012

В идеале вам следует использовать пакет MDAnalysis для питонического способа "нарезки" атомов и сегментов и расчета мер расстояния между ними. . Фактически, MDAnalysis поддерживает несколько форматов моделирования МД и химической структуры.

Более подробный пример см. также в следующей записи на Biostars.org.

person Paul Rigor    schedule 02.05.2013

Если вас интересует только одна пара, 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
person PhysicalChemist    schedule 03.09.2014

Простой код 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

Пожалуйста, дайте мне знать, если вы обнаружите какие-либо несоответствия или проблемы при использовании этого скрипта.

person Gaurav Singhal    schedule 14.12.2017