python - 计算原子坐标之间的距离

标签 python perl biopython bioperl

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

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 进行此计算?

最佳答案

不要按空格分割

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

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

PDB.Residue 对象允许通过名称对 Atom 进行键控访问,并且 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!"

关于python - 计算原子坐标之间的距离,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/13645439/

相关文章:

arrays - 创建对哈希值的数组引用

python - 如何使用python编程将一组DNA序列转换成蛋白质序列?

python - gridfs "list"方法返回具有非空集合的空列表

python - 如何不使用 gmail 从 python 发送电子邮件?

Python ctypes 和没有足够的参数(缺少 4 个字节)

python - 如何解决 Pymongo 的这个错误?游标 ID 找不到 pymongo

Perl:将 STDERR 重定向到文件而不创建空文件?

linux - 如何允许服务器上的任何用户编辑特定文件?

python - 如何为多个文件运行Python脚本?

python - 从 SeqIO.index 生成的字典中删除一个项目