如何用awk，python或biopython從PDB文件獲取子單元的數量？

我有一個目錄中的PDB（文本）文件。我想打印每個PDB文件中的子單元的數量。如何用awk，python或biopython從PDB文件獲取子單元的數量？

閱讀PDB文件與ATOM
的ATOM行的第五列包含A，B，C，D等
如果它只包含A亞基的數量開始的所有行1.如果它包含A和B，亞基的數目是2.如果它包含A，B和C，亞基的數目是3.

1kg2.pdb文件

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文件

ATOM 1384 O ARG A 260  80.505 -20.450 15.420 1.00 22.10   O 
ATOM 1385 CB ARG A 260  78.980 -18.077 15.207 1.00 36.88   C 
ATOM 1399 SD MET B 196  34.003 -52.544 16.664 1.00 57.16   S 
ATOM 1401 N ASP C 197  34.781 -50.611 12.007 1.00 44.30   N

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做到這一點？

來源

2012-09-15 user1672711

歡迎的StackOverflow！你有什麼嘗試？你有沒有看過BioPython中的'Bio.PDB.PDBParser'類？有文檔[這裏]（http://biopython.org/DIST/docs/tutorial/Tutorial.html#htoc138），甚至有一篇關於它的論文（http://www.ncbi.nlm.nih.gov/pubmed）/14630660） –

您可以使用array來記錄第五列的所有可見值。

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

編輯：可以使用ENDFILE使用GAWK 4.x版來生成所需的輸出：

BEGIN { 
    print "pdb_id\t\tsubunits" 
    print 
} 

/^ATOM/ { 
    seen[$5] = 1 
} 

ENDFILE { 
    print FILENAME, "\t", length(seen) 
    delete seen 
}

結果：

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

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

來源

2012-09-15 15:04:40

重定向不是必需的。 –

我已將其刪除。 – 2012-09-15 17:34:35

非常感謝您的回答。代碼只適用於單一。我嘗試了以下代碼爲多個files.gawk'/^ATOM/{看到[$ 5] = 1}結束{打印長度（看到）}'文件/ *。pdb.But它沒有工作。如何更改多個文件的代碼？我還需要打印帶有子單元的文件名稱。 – user1672711

一個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/。後續問題應包括證據表明，用戶實際已經試圖解決這個問題。）

來源

2012-09-15 15:42:18 Dave

如何用awk，python或biopython從PDB文件獲取子單元的數量？

回答

相關問題