如何使用awk、python或biopython从PDB文件中获取子单元的数量？

Question

我有一个目录中的PDB（文本）文件。我想从每个PDB文件打印子单元的数量。在

1. 读取pdb文件中以ATOM开头的所有行
2. ATOM行的第五列包含A、B、C、D等
3. 如果它只包含A，那么子单元的数目是1。如果它含有A和{}，那么亚基的数目是2。如果它含有A、B、和{}，那么亚基的数目是3。在

pd1kb2文件

ATOM   1363  N   ASN A 258      82.149 -23.468   9.733  1.00 57.80           N  
ATOM   1364  CA  ASN A 258      82.494 -22.084   9.356  1.00 62.98           C  
ATOM   1395  C   MET B 196      34.816 -51.911  11.750  1.00 49.79           C  
ATOM   1396  O   MET B 196      35.611 -52.439  10.963  1.00 47.65           O  

1uz3.pdb文件

^{pr2}$

2b69.pdb文件

ATOM   1393  N   MET B 196      33.300 -54.017  12.033  1.00 46.46           N  
ATOM   1394  CA  MET B 196      33.782 -52.714  12.566  1.00 49.99           C  

期望输出

pdb_id   subunits

 1kg2      2
 1uz3      3
 2b69      1

如何使用awk、python或Biopython来实现这一点？在

Answer 1

您可以使用array来记录第五列的所有可见值。在

$ gawk '/^ATOM/ {seen[$5] = 1} END {print length(seen)}' 1kg2.pdb
2

编辑：使用gawk 4.x，您可以使用^{}生成所需的输出：

^{pr2}$

结果是：

$ gawk -f pdb.awk 1kg2.pdb 1uz3.pdb 2b69.pdb
pdb_id          subunits

1kg2.pdb         2
1uz3.pdb         3
2b69.pdb         1

Answer 2

dictionary是计算唯一出现次数的一种方法。下面为每个子单元分配一个无意义的值（0），因为您只关心唯一子单元（字典键）的数量。在

import os

for fn in os.listdir():
    if ".pdb" in fn:
        sub = {}
        with open(fn, 'r') as f:
            for line in f:
                c = line.split()
                if len(c) > 5 and c[0] == "ATOM":
                    sub[c[4]] = 0
        print(fn, len(sub.keys()))

（一个全新的用户应该得到一个答案以及一个指向http://whathaveyoutried.com/的指针。后续问题应包括用户实际尝试解决问题的证据。）