Filling missing loops from a homologous structure#

You will learn: how to fill unresolved (missing-residue) regions of a crystal structure by grafting the missing residues from a more complete donor structure of the same protein, using spliceMissingResidues(), and how to find and check those gaps with resolveFullSequences(), detectSequenceGaps(), and detectSplicedClashes().

Prerequisites:

Setup#

from moleculekit.molecule import Molecule
from moleculekit.rcsb import resolveFullSequences
from moleculekit.tools.modelling import (
    detectSequenceGaps,
    spliceMissingResidues,
    detectSplicedClashes,
)

Step 1 - Load a structure with unresolved regions#

We use 5VQ2, a GTP/Mg-bound KRAS homodimer (chains A and B) whose crystal has several unresolved regions per chain. We load it as deposited, keeping its cofactors and crystallographic waters:

mol = Molecule("5VQ2")

moleculekit.molecule - WARNING - Alternative atom locations detected. Only altloc A was kept. If you prefer to keep all use the keepaltloc="all" option when reading the file.
moleculekit.molecule - INFO - Removed 62 atoms. 2675 atoms remaining in the molecule.

Step 2 - Resolve the full sequence and find the gaps#

resolveFullSequences() returns the full deposited sequence for each protein chain (keyed by chain, each value a dict with sequence, source, and identity). detectSequenceGaps() takes a plain {chain: sequence} mapping, so we unwrap the sequence field first:

resolved = resolveFullSequences(mol, "5VQ2")
sequences = {ch: v["sequence"] for ch, v in resolved.items()}

gaps, skipped = detectSequenceGaps(mol, sequences)
for g in gaps:
    print(g)

{'chain': 'A', 'after_resid': 33, 'before_resid': 37, 'missing_seq': 'PTI', 'is_terminal': False}
{'chain': 'A', 'after_resid': 58, 'before_resid': 69, 'missing_seq': 'AGQEEYSAMR', 'is_terminal': False}
{'chain': 'A', 'after_resid': 168, 'before_resid': None, 'missing_seq': 'K', 'is_terminal': True}
{'chain': 'B', 'after_resid': 33, 'before_resid': 37, 'missing_seq': 'PTI', 'is_terminal': False}
{'chain': 'B', 'after_resid': 58, 'before_resid': 69, 'missing_seq': 'AGQEEYSAMR', 'is_terminal': False}
{'chain': 'B', 'after_resid': 168, 'before_resid': None, 'missing_seq': 'K', 'is_terminal': True}

Each gap dict gives the chain, the observed resid just before (after_resid) and after (before_resid) the gap, the missing subsequence (missing_seq), and whether it is a chain terminus (is_terminal). Chains A and B show the identical set of three gaps: a 3-residue internal loop between resid 33 and 37 (PTI), a 10-residue internal loop between resid 58 and 69 (AGQEEYSAMR, the switch-II region), and a single unresolved residue after resid 168 (K) with no observed residue following it, so before_resid is None and is_terminal is True.

Step 3 - Graft the missing residues from a donor crystal#

4DSO is the same protein in a more complete crystal. spliceMissingResidues() pairs each original chain to its best-matching donor chain by sequence identity and inserts only the residues present in the donor but missing from our structure. All original atoms keep their deposited coordinates except the graft_flanks residues on each side of a filled gap, which are taken from the donor to seat the insertion:

donor = Molecule("4DSO")
sp1, new1 = spliceMissingResidues(mol, donor, graft_flanks=1)

moleculekit.tools.modelling - WARNING - Skipping 12 modelled residue(s) in chain A after resid 168: the donor could not be attached with a physical backbone junction within graft_flanks=1. Re-run spliceMissingResidues on the result with a larger graft_flanks to graft deeper and retry the remaining gaps (a deeper graft may introduce sidechain clashes to relax downstream).
moleculekit.tools.modelling - WARNING - Skipping 12 modelled residue(s) in chain B after resid 168: the donor could not be attached with a physical backbone junction within graft_flanks=1. Re-run spliceMissingResidues on the result with a larger graft_flanks to graft deeper and retry the remaining gaps (a deeper graft may introduce sidechain clashes to relax downstream).

With graft_flanks=1 the loops whose flanking residues already agree between the two crystals splice cleanly, but the flexible switch-II / C-terminal region diverges too far for a single-residue graft and is skipped with a warning (look for a Skipping ... after resid 168 message above). spliceMissingResidues returns the new Molecule and a boolean mask marking the newly added atoms.

Step 4 - Re-run with a deeper graft to fill the rest#

Re-running on the result with a larger graft_flanks grafts deeper into each flank, which closes the region the minimal graft skipped:

sp2, new2 = spliceMissingResidues(sp1, donor, graft_flanks=2)

Both chains are now full length and continuous. The two GTP molecules and both Mg ions are carried through untouched:

for ch in ("A", "B"):
    n = len({int(r) for r in sp2.resid[sp2.atomselect("protein") & (sp2.chain == ch)]})
    print(f"chain {ch} residues:", n)
print("GTP atoms:", int((sp2.resname == "GTP").sum()), " MG atoms:", int((sp2.resname == "MG").sum()))

chain A residues: 181
chain B residues: 181
GTP atoms: 64  MG atoms: 2

Step 5 - Check the modelled residues for clashes#

A grafted loop can land in a cofactor, ion, or water site. detectSplicedClashes() checks the newly added atoms (the mask returned by the splice) against the non-protein components, which here include the crystallographic waters:

clashes = detectSplicedClashes(sp2, new2)
print(clashes)

[{'new_chain': 'A', 'new_resid': 169, 'target_resname': 'HOH', 'target_resid': 308, 'target_segid': '4', 'min_distance': 1.3294672966003418}, {'new_chain': 'A', 'new_resid': 172, 'target_resname': 'HOH', 'target_resid': 343, 'target_segid': '4', 'min_distance': 1.9198095798492432}, {'new_chain': 'B', 'new_resid': 169, 'target_resname': 'HOH', 'target_resid': 305, 'target_segid': '4', 'min_distance': 1.3065155744552612}, {'new_chain': 'B', 'new_resid': 172, 'target_resname': 'HOH', 'target_resid': 335, 'target_segid': '4', 'min_distance': 1.7084144353866577}]

Here the run reports four overlaps: the deeper graft placed the modelled C-terminal residues 169 and 172 (in both chains) within the 2 Angstrom cutoff of crystallographic waters. Each entry names the offending grafted residue and the atom it overlaps, with the minimum distance. Overlaps with waters are expected when loops are grafted into a structure that still carries its deposited solvent; they are resolved downstream by removing or re-placing the affected waters, or by a short minimization, and finding them is exactly why this check runs before you proceed. To restrict the check to cofactors and ions only, pass targets="resname GTP MG".

The filled structure is now ready to hand to systemPrepare() (see the Basic protonation tutorial); sidechain strain introduced by a deeper graft is left for downstream minimization.

Recap#

  • resolveFullSequences() gets the full deposited sequence per chain; unwrap its sequence field before passing a {chain: sequence} map to detectSequenceGaps().

  • detectSequenceGaps() returns (gaps, skipped_chains); each gap locates a missing-residue run between two observed resids.

  • spliceMissingResidues() grafts the residues missing from your structure out of a more complete donor of the same protein, returning (mol, new_mask); increase graft_flanks and re-run to close regions a minimal graft skips.

  • detectSplicedClashes() flags newly modelled atoms that overlap cofactors, ions, or ligands.

Next#