【问题标题】：Calculating the distance between atomic coordinates计算原子坐标之间的距离
【发布时间】：2012-11-30 12:37:47
【问题描述】：

我有一个如下所示的文本文件

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

我想计算两个α碳原子之间的距离，即计算第一个和第二个原子之间的距离，然后计算第二个和第三个原子之间的距离等等.....两个原子之间的距离可以表示为：@ 987654322@

第 7,8 和 9 列分别代表 x,y 和 z 坐标。我需要打印距离和对应的残差对（第 4 列），如下所示。（距离的值不是真实的）

GLN-HIS   4.5
HIS-ASN   3.2
ASN-GLU   2.5

如何使用 perl 或 python 进行此计算？

【问题讨论】：

我有同样的问题，但是我正在使用python，那里也有类似的解决方案吗？

标签： python perl biopython bioperl

【解决方案1】：

不要在空格上分割

这里给出的其他答案做出了一个有缺陷的假设 - 坐标将以空格分隔。根据PDB specification of 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 格式很乱，有很多特殊情况和格式错误的文件需要处理。相反，请使用已经为您编写的解析器。

我喜欢Biopython 的PDB.PDBParser。它将为您将结构解析为 Python 对象，并具有方便的功能。如果您更喜欢 Perl，请查看 BioPerl。

PDB.Residue 对象允许按名称对 Atom 进行键控访问，PDB.Atom 对象重载 - 运算符以返回两个 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!"

【讨论】：

【解决方案2】：

如果您的数据由空格分隔，一个简单的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/ 而不是 ' ' 上拆分。

【讨论】：

非常感谢您的回答。您的代码运行良好。但我没有得到逐行的输出。我想像您的示例输出一样更改输出的外观。我该如何更改？感谢您的帮助！
哦，对了，忘记了。我在测试时使用-l 标志运行脚本。 perl -l script.pl input.txt。我会解决的。
重要提示：PDB 原子坐标不一定以空格分隔（请参阅specification）。
而不是计算从一个 CA 到下一个 CA 的距离等等.. 我们如何修改上面的答案来计算一个 CA 到所有其他 CA 之间的距离，然后计算从下一个 CA 到其余 CA 的距离CA 等等...？

【解决方案3】：

假设您的数据在“atoms.txt”中，这将逐行读取并将条目拆分为列表：

import csv

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

现在为每个列表提取您需要的列，并计算距离等（请记住，python 中的列表是从零开始的）。

【讨论】：

你测试过这段代码吗？不要介意for line, in reader: 的错误——返回的将是list，它没有split 方法...
和我第一次说的完全一样...您将收到 for line, in reader 的 ValueError，一旦修复 - 您将在 entries = line.split() 上收到 AttributeError: 'list' object has no attribute 'split'
不，我不明白。复制粘贴问题中文本文件的前五行，并根据需要运行我的答案中的代码会打印一个列表。（这里是python 2.6，不知道python 3 --- csv.reader 有什么不同吗？）
尝试将前五行放入atoms.txt 并完全按照您发布的方式运行代码...
就像我说的，这正是我所做的，它会打印令牌列表。

【解决方案4】：

理想情况下，您应该使用MDAnalysis package 来“切片”原子和段并计算它们之间的距离度量。事实上，MDAnalysis 支持多种 MD 模拟和化学结构格式。

有关更详细的示例，另请参阅following entry on Biostars.org。

【讨论】：

【解决方案5】：

如果您只对一对感兴趣，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

【讨论】：

【解决方案6】：

一个简单的 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

如果您在使用此脚本时发现任何差异或问题，请告诉我。

【讨论】：