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

Question

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

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

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?

3

python perl biopython bioperl

Источник

user1866262 30 ноя '12 в 12:37

6 ответов

Решение

НЕ разбивать на пробелы

Другие ответы, приведенные здесь, делают ошибочное предположение - что координаты будут разделены пробелом. Согласно спецификации 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!"

6

Источник

user815632 01 дек '12 в 17:07

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

import csv

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

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

4

Источник

user143765 30 ноя '12 в 12:44

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

Для более подробного примера см. Также следующую запись на Biostars.org.

3

Источник

user565946 02 май '13 в 07:26

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

1

Источник

user3767980 03 сен '14 в 00:52

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

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

0

Источник

user3481574 14 дек '17 в 18:55

Другие вопросы по тегам python perl biopython bioperl

user725418 30 ноя '12 в 13:50 2012-11-30 13:50 · Accepted Answer · 2012-11-30 13:50

Если ваши данные разделены пробелами, простой 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/ вместо ' ',