在蛋白质结构文件中重新编号残基(pdb)
你好!
我现在正在做一个网站,目的是把所有关于乳头瘤病毒的信息集中到一个地方。
在这个过程中,我们正在整理所有在公共服务器上已知的文件(比如genbank)。
我遇到的一个问题是,很多(大约50%)已经解决的结构没有按照蛋白质的顺序编号。
比如说,一个亚域的晶体结构是由氨基酸310到450构成的,但晶体学家却把它标记成了残基1到140。
我想知道有没有人知道如何重新编号整个pdb文件。我找到了一些方法可以重新编号序列(通过seqres识别),但这并不能更新螺旋和片层的信息。
如果你有任何建议,我会非常感激……
谢谢!
3 个回答
7
我是pdb-tools的维护者,这个工具可能对你有帮助。
我最近对我应用中的residue-renumber
脚本进行了修改,让它变得更加灵活。现在它可以对hetatms
和特定链
进行重新编号,用户可以选择让残基编号连续,或者只是给所有残基加上一个自己指定的偏移量。
如果这个工具对你有帮助,请告诉我。
1
我也经常遇到这个问题。在放弃了我之前用的一个老旧的perl脚本后,我开始尝试用python来解决这个问题。这个方法假设你已经安装了Biopython、ProDy(http://www.csb.pitt.edu/ProDy/#prody)和EMBOSS(http://emboss.sourceforge.net/)。
我使用了这里的一个乳头瘤病毒的PDB条目。
from Bio import AlignIO,SeqIO,ExPASy,SwissProt
from Bio.Seq import Seq
from Bio.SeqRecord import SeqRecord
from Bio.Alphabet import IUPAC
from Bio.Emboss.Applications import NeedleCommandline
from prody.proteins.pdbfile import parsePDB, writePDB
import os
oneletter = {
'ASP':'D','GLU':'E','ASN':'N','GLN':'Q',
'ARG':'R','LYS':'K','PRO':'P','GLY':'G',
'CYS':'C','THR':'T','SER':'S','MET':'M',
'TRP':'W','PHE':'F','TYR':'Y','HIS':'H',
'ALA':'A','VAL':'V','LEU':'L','ILE':'I',
}
# Retrieve pdb to extract sequence
# Can probably be done with Bio.PDB but being able to use the vmd-like selection algebra is nice
pdbname="2kpl"
selection="chain A"
structure=parsePDB(pdbname)
pdbseq_str=''.join([oneletter[i] for i in structure.select("protein and name CA and %s"%selection).getResnames()])
alnPDBseq=SeqRecord(Seq(pdbseq_str,IUPAC.protein),id=pdbname)
SeqIO.write(alnPDBseq,"%s.fasta"%pdbname,"fasta")
# Retrieve reference sequence
accession="Q96QZ7"
handle = ExPASy.get_sprot_raw(accession)
swissseq = SwissProt.read(handle)
refseq=SeqRecord(Seq(swissseq.sequence,IUPAC.protein),id=accession)
SeqIO.write(refseq, "%s.fasta"%accession,"fasta")
# Do global alignment with needle from EMBOSS, stores entire sequences which makes numbering easier
needle_cli = NeedleCommandline(asequence="%s.fasta"%pdbname,bsequence="%s.fasta"%accession,gapopen=10,gapextend=0.5,outfile="needle.out")
needle_cli()
aln = AlignIO.read("needle.out", "emboss")
os.remove("needle.out")
os.remove("%s.fasta"%pdbname)
os.remove("%s.fasta"%accession)
alnPDBseq = aln[0]
alnREFseq = aln[1]
# Initialize per-letter annotation for pdb sequence record
alnPDBseq.letter_annotations["resnum"]=[None]*len(alnPDBseq)
# Initialize annotation for reference sequence, assume first residue is #1
alnREFseq.letter_annotations["resnum"]=range(1,len(alnREFseq)+1)
# Set new residue numbers in alnPDBseq based on alignment
reslist = [[i,alnREFseq.letter_annotations["resnum"][i]] for i in range(len(alnREFseq)) if alnPDBseq[i] != '-']
for [i,r] in reslist:
alnPDBseq.letter_annotations["resnum"][i]=r
# Set new residue numbers in the structure
newresnums=[i for i in alnPDBseq.letter_annotations["resnum"][:] if i != None]
resindices=structure.select("protein and name CA and %s"%selection).getResindices()
resmatrix = [[newresnums[i],resindices[i]] for i in range(len(newresnums)) ]
for [newresnum,resindex] in resmatrix:
structure.select("resindex %d"%resindex).setResnums(newresnum)
writePDB("%s.renumbered.pdb"%pdbname,structure)