The residue numbering in PDB (in my understanding) coincides with UniProtKB align. However, this might not be consistent with that in PDB fasta file. To complicate it further, UniProtKB align can be a few mutations away from the PDB sequence. Disordered/unidentified residues in structures also create gaps in numbering.

Is there a way/suggestion to create a matched residue numbering between sequence and structure? For example

12345678 (residue numbering)
ATPWQMSQ (sequence)
---W-MSQ (PDB residues) "-" for missing residue in structure

Unfortunately, PDB numbering can be one of different schemes. The main three:

  • Uniprot numbering (canonical isoform)
  • Construct numbering with or without tag removal
  • Start from 1 for what has density

However, the are many many corner cases. In Proteopedia there is a lovely article about weird PDB numbering.

The best way to find the details about what part of the sequence is covered is to use the Swissmodel API. It's actually better than either the PDBe and RCSB PDB APIs. Swissmodel is a homology modelling service that maintains close threaded models (>30% ID), but also has PDB data. Here is a Python3 method I use, (provider='pdb') is what is required here).

def get_data(uniprot: str, provider='swissmodel') -> List[dict]:
        Gets data from Swissmodel. Convert each entry with convert_SM_entry
    assert provider in ('swissmodel', 'pdb'), 'Provider has to be swissmodel or pdb'
        url = f'https://swissmodel.expasy.org/repository/uniprot/{uniprot.strip()}.json?provider={provider.strip()}'
        data = requests.get(url).json()
        return data['result']['structures']  # [0] # from to coordinates description

The returned dictionary contains a start and end of the uniprot <--> PDB. It does not however, say where the missing density residues are. But in deposited structures the PDB numbering skips where these residues ought to be... So might not be an important problem.

Another database of interest is SIFTS.

To make a multiple sequence alignment of a non-deposited PDB, you need to extract the sequence and do the alignment here is a snippet using pymol module (from conda):

import pymol2

pdbfile = '👾👾👾' #your pdb file name..
selector = 'chain A' # or whatever chain you want
with pymol2.PyMOL() as pymol:
    fasta = pymol.cmd.get_fastastr(selector)

seq = ''.join(fasta.replace('?','').split('\n')[1:]) # remove header and non AA residues

Biopython has a dozen sequence alignment tools. But the original one still works fine

from Bio import pairwise2
from Bio.Seq import Seq

alignment = pairwise2.align.globalxx(''.join(seq).replace('?',''), "ACQQG")[0]

alignments is a tuple of the MSA sequences.

If you want to do the latter without coding. Extract the sequence from the PDB using PyMOL the normal application and load the pdb and then print cmd.get_fastastr('chain A') in the command line part of the GUI. Copy this. Go to Muscle or other onliner aligner and align!

  • $\begingroup$ Thanks @MatteoFerla. The swissmodel API is an interesting one. However, since the rcsb fasta does not have a UniProtKB id, perhaps I have to do an additional alignment to obtain the UniProtKB id? I am surprised rcsb has the sequence data (with the alignment, UnitProtKB id, missing residues, etc), but there is no API/downloadable available. $\endgroup$ – Simon Aug 23 '20 at 16:26
  • 1
    $\begingroup$ To go from PDB 4-letter code to Uniprot accession use the SIFTS API or mapping table from their FTP site. Or look manually at the PDB entry. Do note that a uniprot record has an accession number, say for HRas P01112, and a entry name RASH_HUMAN. The "id" word is an ambiguous and can refer to either, but use to mean the latter. To confuse further, the entry name can change if it is moved from Trembl to Swissprot —say the genome of that organism was sequenced and before the protein was known from an EST read or something. $\endgroup$ – Matteo Ferla Aug 24 '20 at 10:14
  • 1
    $\begingroup$ The RCSB PDB has changed recently, dropping a few annotation from its protein viewer, which I am guessing will be re-added —once that happens I assume a sensible API will appear. How it used to work was that HTML/JS chunks were assembled in the backend, so there was no API for uniprot data and you could not reverse-engineer an ajax call into an API. $\endgroup$ – Matteo Ferla Aug 24 '20 at 10:18
  • $\begingroup$ Thank you so much. It is indeed a very comprehensive explanation. $\endgroup$ – Simon Aug 24 '20 at 16:58

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.