我也经常遇到这个问题.在放弃旧的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)