如何從 PDB 文件中提取氨基酸序列
蛋白質(zhì)數(shù)據(jù)庫Protein Data Bank(PDB)是一個包含蛋白質(zhì)、核酸等生物大分子的結(jié)構(gòu)數(shù)據(jù)的數(shù)據(jù)庫，網(wǎng)址是http://www.rcsb.org。PDB可以經(jīng)由網(wǎng)絡(luò)免費訪問，是結(jié)構(gòu)生物學研究中的重要資源。為了確保PDB資料的完備與權(quán)威，各個主要的科學雜志、基金組織會要求科學家將自己的研究成果提交給PDB。PDB數(shù)據(jù)庫存儲結(jié)構(gòu)數(shù)據(jù)的文件是PDB文件，每一個蛋白質(zhì)或核酸都對應(yīng)著一個編號，即PDBID， 文件的擴展名為.pdb。PDB文件可以由各種3D結(jié)構(gòu)顯示軟件打開，比如pymol，Swiss-PDB viewer，VMD等。PDB文件里面的信息是有嚴格的格式的。各行數(shù)據(jù),如標識,原子名,原子序號,殘基名稱,殘基序號等,不僅要按照嚴格的順序書寫,而且各項所占的空符串長度,及其所處的各行的位置都是嚴格規(guī)定。更多信息請查看PDB文件中信息的格式。
寫在前面的話：
本人是一枚生物學的學生，由于對生物信息學特別感興趣，于是想自學生物信息學（新手莫怪）。了解到生物信息學要有編程基礎(chǔ)，尤其是要會一門編程語言，例如：R語言、Python、Perl等，還要熟悉Linux系統(tǒng)，作為生信小白，聽說Python挺簡單的，于是就自學了Python，花了兩天時間了解了Python的基礎(chǔ)語法后，今天想做個練習題試試手（實踐是檢驗真理的唯一標準）。
氨基酸殘基與核酸縮寫
眾所周知氨基酸標準縮寫是三個字母，而這次我們所用的PDB文件里面也是用三個字母表示的氨基酸，我們需要把它轉(zhuǎn)化為常見的fasta格式的氨基酸序列，所以我們首先應(yīng)該用字典將三個字母的氨基酸轉(zhuǎn)換為單字母代碼。
字典的鍵是氨基酸三字母代碼， 值是對應(yīng)的氨基酸的單字母代碼。使程序用該字典從 PDB 文件的 SEQRES 行讀取殘基名稱 (以三字母代碼的形式) (文件 1TDL 中) , 然后將它們轉(zhuǎn)換為一個字母代碼，連接這些字母，就得到了蛋白質(zhì)序列。最后，將該序列以 FASTA 的格式打印出來。
1TLD. pdb文件如圖所示：
要使用該程序，需要找到 PDB文檔 ，下載并保存 PDB文件。 此例使用了 PDB文件 1TLD. pdb，其中有牛自膜蛋白酶在1. 5A分辨率下的晶體結(jié)構(gòu)，該 PDB文件 的 SEQ阻S行列出了實驗用到的蛋白質(zhì)序列。每個 SEQRES行可分為 17 列。 第一列包含關(guān)鍵 字"SEQRES"; 第二列有序列的行號(從 1 開始) ;第三列是鏈 ID; 第四列是組成序列的殘基數(shù); 從第五列到最后一列是(三字母代碼形式的)殘基。列之間由一個空格進行分隔。
第一步：定義字典
aa_codes = {
'ALA':'A','CYS':'C','ASP':'D','GLU':'E',
'PHE':'F','GLY':'G','HIS':'H','LYS':'K',
'ILE':'I','LEU':'L','MET':'M','ASN':'N',
'PRO':'P','GLN':'Q','ARG':'R','SER':'S',
'THR':'T','VAL':'V','TYR':'Y','TRP':'W'}
第二步：讀取文件
首先讀取文件，以"SEQRES"定位到位置，并.split()去掉空格
seq = ''
for line in open("C:/Desktop/1tld.pdb"):
if line[0:6] =="SEQRES":
columns = line.split()
for resname in columns[4:]:
seq = seq + aa_codes[resname]
第三步：打印輸出*
i = 0
print(">1TLD")
while i < len(seq) :
print(seq[i:i+64])
i =i+64
最后附上完整代碼：
aa_codes = {
'ALA':'A','CYS':'C','ASP':'D','GLU':'E',
'PHE':'F','GLY':'G','HIS':'H','LYS':'K',
'ILE':'I','LEU':'L','MET':'M','ASN':'N',
'PRO':'P','GLN':'Q','ARG':'R','SER':'S',
'THR':'T','VAL':'V','TYR':'Y','TRP':'W'
}
seq = ''
for line in open("C:/Desktop/1tld.pdb"):
if line[0:6] =="SEQRES":
columns = line.split()
for resname in columns[4:]:
seq = seq + aa_codes[resname]
i = 0
print(">1TLD")
while i < len(seq) :
print(seq[i:i+64])
i =i+64
日常結(jié)尾：
雖然這是個小小的計算程序，但對于初學者的我來說每一次對原代碼的升級改造，哪怕是讀懂后的注釋都感覺是一次進步提升，總之代碼雖小，動手最重要！希望更多學習Python的愛好者不要像我一樣眼高手低，學習編程就是要，思考，敲碼，思考，敲碼，敲碼，再敲碼！
