2013-07-20 68 views
1

我正在解析PDB文件,我有一个链式名称列表以及格式(链,[坐标])中的XYZ坐标。我有很多坐标,但只有3个不同的链条。我想将所有来自同一链的坐标压缩到一个列表中,以便获得链= [坐标],[坐标],[坐标]等等。我查看了biopython文档,但我很难理解如何获得我想要的坐标,所以我决定手动提取坐标。这是我到目前为止的代码:分隔(X,Y)的列表

pdb_file = open('1adq.pdb') 
import numpy as np 

chainids = [] 
chainpos= [] 

for line in pdb_file: 
    if line.startswith("ATOM"): 
     # get x, y, z coordinates for Cas 
     chainid =str((line[20:22].strip())) 
     atomid = str((line[16:20].strip())) 
     pdbresn= int(line[23:26].strip()) 
     x = float(line[30:38].strip()) 
     y = float(line[38:46].strip()) 
     z = float(line[46:54].strip()) 
     if line[12:16].strip() == "CA": 
      chainpos.append((chainid,[x, y, z])) 
      chainids.append(chainid) 

allchainids = np.unique(chainids) 
print(chainpos) 

和一些输出:

[('A', [1.719, -25.217, 8.694]), ('A', [2.934, -21.997, 7.084]), ('A', [5.35, -19.779,  8.986]) 

我的理想输出将是:

A = ([1.719, -25.217, 8.694]), ([2.934, -21.997, 7.084]),(5.35, -19.779,8.986])... 

谢谢!

Here is a section of PDB file: 
ATOM  1 N PRO A 238  1.285 -26.367 7.882 0.00 25.30   N 
ATOM  2 CA PRO A 238  1.719 -25.217 8.694 0.00 25.30   C 
ATOM  3 C PRO A 238  2.599 -24.279 7.885 0.00 25.30   C 
ATOM  4 O PRO A 238  3.573 -24.716 7.275 0.00 25.30   O 
ATOM  5 CB PRO A 238  2.469 -25.791 9.881 0.00 25.30   C 

A是链名称有第4栏,我不知道是什么的链名是先验的,但因为我的行解析线,我与格式坐标粘链的名字我之前提到。现在我想把所有的坐标都放在一个“A”之前,并把它们放在一个名为“A”的列表中。我不能在“A”中硬编码,因为它不总是“A”。我也有“L”和“H”,但我想我可以让他们一旦我得到了理解驼峰..

+0

你想创建变量'A'。 ,或打印出“A =([1.719,-25.217,8.694]),([2.934,-21.997,7.084]),(5.35,-19.779,8.986))...” – inspectorG4dget

+0

此外,请张贴PDB文件相关部分的代表性样本,以便我们知道您要解析的内容 – inspectorG4dget

+0

有一个名为pdb_tools的密集型库,用Python编写。 – chupvl

回答

1

只是让元组的列表

>>> chainpos.append((chainid,x, y, z)) 
>>> chainpos 
[('A', 1.719, -25.217, 8.694), ('A', 2.934, -21.997, 7.084)] 
>>> import itertools 
>>> for id, coor in itertools.groupby(chainpos,lambda x:x[0]): 
...  print(id, [c[1:] for c in coor]) 
+0

这正是我一直在寻找的!我喜欢它维持秩序,因为它可能在未来某个时候变得重要。 – pioneer903

+0

有没有办法使用这个itertool将它们保存到自己的列表中? – pioneer903

+0

只需更改此行并将其附加到新列表中即可。 print(id,[c [1:] for c in coor] [0])。这是你想要的 – raton

0

您可以使用列表理解:

>>> print chainpos 
[('A', [1.719, -25.217, 8.694]), ('A', [2.934, -21.997, 7.084]), ('A', [5.35, -19.779,  8.986])] 
>>> print "A =", [ t[1] for t in chainpos] 
+0

如果它不是“A”呢?目前还不清楚我们是否知道该字母是什么,先验 – inspectorG4dget

+0

A只是象征原子所属的肽链,在这种情况下,它是脯氨酸。所以通常'A'不能保证。 – seth

+0

恰好Seth。在PDB文件中,链总是连续列出的,也就是说,一旦有一个字母来表示链,那个链字母就会被保留下来,直到下一个链开始。在这种情况下,它会去是这样的: 一个 一个 一个 一个 大号 大号 ^ h ^ h 在我的代码存储所有在allchainids链的名字,因此具有“[ 'A'“L ''H']在这里。 – pioneer903

1

你想要的东西像:

import numpy as np 

chain_dict = {} 

for line in open('input'): 
    if line.startswith("ATOM"): 
     line = line.split() 
     # get x, y, z coordinates for Cas 
     chainid = line[4] 
     atomid = line[2] 
     pdbresn= line[5] 
     xyz = [line[6],line[7],line[8]] 
     if chainid not in chain_dict: 
      chain_dict[chainid]=[xyz] 
     else: 
      chain_dict[chainid].append(xyz) 

其中,为您的数据。例如,给出:

>>> chain_dict 
{'A': [['1.285', '-26.367', '7.882'], ['1.719', '-25.217', '8.694'], ['2.599', '-24.279', '7.885'], ['3.573', '-24.716', '7.275'], ['2.469', '-25.791', '9.881']] 

,并因为它是一本字典,很明显,你可以d ○:

>>> chain_dict['A'] 
[['1.285', '-26.367', '7.882'], ['1.719', '-25.217', '8.694'], ['2.599', '-24.279', '7.885'], ['3.573', '-24.716', '7.275'], ['2.469', '-25.791', '9.881']] 

得到的只是你感兴趣链的XYZ COORDS

+0

这看起来很有趣。我将不得不花一点时间理解Python字典..它们比列表更复杂。 – pioneer903