Protein Target Preparation (before Blastermaster)

Protein preparation for docking (prepare_protein.py)

prepare_protein.py is a Python 3 script that automates the receptor and ligand preparation procedure that is otherwise done by hand in Chimera + vim before running Blastermaster. Starting from the Maestro-minimized files (rec_minimized_final.pdb and xtal_minimized_final.pdb), it produces the cleaned, charged, correctly-named PDB files and the working/ directory that Blastermaster expects.

It folds in the logic of the two legacy helper scripts (replace_his_with_hie_hid_hip.py and 0000_remove_hydrogens_from_pdb.py), so those no longer need to be run separately.

Requirements

• Python 3 (no third-party packages; standard library only)
• Two input PDB files that have already been through the manual Chimera steps

 (open structure, add hydrogens / protonate with reduce, check termini and the
 protonation states of charged residues, then save):

• • the receptor, saved as rec_minimized_final.pdb
  • the crystallographic ligand, saved as xtal_minimized_final.pdb

The receptor file must still contain its hydrogens (including the polar HD1/HE2 on histidines) — these are required to assign the HIS protonation states. The script removes hydrogens itself at the correct stages.

Quick start

On Gimel, run in the directory that holds the two input files:

python3 /mnt/nfs/exa/work/ak87/UCSF/SCRIPTS/DOCKING/prepare_protein.py

Or specify paths explicitly:

python3 /mnt/nfs/exa/work/ak87/UCSF/SCRIPTS/DOCKING/prepare_protein.py \
    --rec  rec_minimized_final.pdb \
    --xtal xtal_minimized_final.pdb \
    --output-dir .

What it does

Ligand pipeline

1. Strips all header, REMARK, CONECT and END lines.
2. Deletes all hydrogens (sphere matching does not use them).
3. Writes xtal-lig.pdb.

Receptor pipeline

1. Checks + fixes residue numbering. If the order of residues in the file does not match their residue numbers, the residue lines are reordered into ascending numeric order (the residue number is treated as ground truth). See Residue numbering.
2. Removes insertion-code letters. A residue such as NMA 378A gets a free number between its neighbours when one exists (378, 378A, 380 → 379), so no other residue is renumbered. Following residues are shifted only when no number is free, and only until the next numbering gap.
3. Strips headers. Removes HEADER / TITLE / REMARK / ANISOU / CONECT / TER / END records, keeping only atoms.
4. "del HC". Deletes hydrogens bonded to carbon (the equivalent of Chimera's del HC); polar hydrogens on N/O/S are kept. See del HC.
5. Caps. Converts ACE and NMA cap records from HETATM to ATOM, and renames the NMA cap CA atom to CM.
6. Modified residues. A non-canonical residue written as HETATM but peptide-bonded into a chain (MSE, SEP, …) is converted to ATOM and treated as protein. See Non-canonical residues.
7. Hetero groups. Waters and hetero groups are dropped unless named with --keep-het; kept ions and cofactors are placed after the protein and written as ATOM records. See Chains and hetero groups.
8. Backbone amide H. Renames the backbone amide hydrogen (H1) of the residue following an ACE cap to H.
9. Disulfides. Detects disulfide-bonded cysteines and renames them CYS → CYX. See Disulfide detection.
10. Writes rec_noHC.pdb.
11. HIS protonation. Assigns HID / HIE / HIP from the reduce hydrogens (HD1 → delta, HE2 → epsilon, both → HIP). Writes rec_noHC.crg.pdb (this is rec.crg.pdb) and copies it into working/.
12. Removes all remaining hydrogens, writing rec_noH.pdb (this is rec.pdb).

Output layout

Inside --output-dir (default = current directory):

xtal-lig.pdb            # ligand, heavy atoms only
rec_noHC.pdb            # receptor, cleaned, carbon-H removed, caps/CYX/numbering fixed
rec_noHC.crg.pdb        # + HIS protonation states assigned  (== rec.crg.pdb)
rec_noH.pdb             # + all hydrogens removed             (== rec.pdb)
rec.pdb                 # final no-hydrogen receptor
working/
    rec.crg.pdb         # the ONLY file placed in working/ (charged receptor for Blastermaster)

None of the receptor files contains a HETATM record: everything kept is written as ATOM. xtal-lig.pdb is the exception — the crystallographic ligand keeps the record type it had in the input, since it is used for sphere matching rather than as receptor.

The two files you normally hand to Blastermaster are working/rec.crg.pdb (charged) and rec.pdb (no hydrogens), together with xtal-lig.pdb.

Command-line options

Option	Default	Description
--rec FILE	rec_minimized_final.pdb	Chimera-minimized receptor PDB.
--xtal FILE	xtal_minimized_final.pdb	Chimera-minimized ligand PDB.
--output-dir DIR	. (current dir)	Directory to write prepared files into.
--cyx "LIST"	(none)	CYS residues to force to CYX, e.g. "32 56 297 338" or chain-qualified "A32 B56". Combined with automatic detection unless --no-auto-disulfide is given.
--no-auto-disulfide	off	Disable automatic SG–SG disulfide detection; use only the residues supplied via --cyx.
--ss-cutoff Å	2.5	SG–SG distance (Å) below which a disulfide is called.
--keep-het "LIST"	(none)	HETATM residue names to keep, e.g. "CA ZN NAG", matched exactly as spelled in the residue-name column. Kept groups are written as ATOM records (Blastermaster deletes HETATM lines). Waters and unlisted hetero groups are dropped. ACE/NMA caps and modified residues bonded into a chain are always kept and converted to ATOM without being listed here.
--keep-carbon-h	off	Keep hydrogens bonded to carbon (skip "del HC"). Polar hydrogens are always kept regardless.
--relabel-chains	off	Relabel chain IDs to sequential letters A, B, C… in order of appearance (fixes gaps such as A, C → A, B).
--unify-chain-id ID	(none)	Give every chain the same ID in the final rec.pdb only; rec.crg.pdb keeps proper chains for Blastermaster.
--no-reorder	off	Do not reorder residues by number. Insertion-code letters are still stripped.

Feature details

Residue numbering

Sometimes the order in which residues appear in a file does not match their residue numbers (for example a residue numbered 14 physically appearing after residue 378). By default the script detects this and reorders the residue lines into ascending numeric order within each chain, treating the residue number as ground truth.

• Numbering gaps (missing residues) are preserved — they are legitimate.
• Insertion codes are then removed, and any residue whose number is not greater than the one before it is renumbered. The script tries to avoid renumbering other residues by picking a free number between the neighbours, in this order:
  • the residue's own base number, if free: 100, 102A, 103 → 100, 102, 103
  • the number right after the preceding residue: 378, 378A, 380 → 378, 379, 380
  • at the start of a chain, the number just below the following residue: 5A, 5, 6 → 4, 5, 6
• Only if no number is free are following residues pushed up, and only until the next numbering gap absorbs the shift: 378, 378A, 379, 380, 390 → 378, 379, 380, 381, 390. Chains are handled independently.
• Use --no-reorder to leave the physical order untouched (letters are still stripped).

The run log reports exactly what moved and what was renumbered, e.g.:

[check] residue numbering does NOT match sequence order in rec_minimized_final.pdb:
  * chain A: LEU 14 (line 8986) is out of order -- its number is not greater than the preceding residue NMA 378A (line 8980)
  -> residues will be reordered by number (see below)
  reordered to numeric position: A/LEU 14
  renumbered (insertion codes removed / free numbers reused; residues shifted only when no number was free):
    A/NMA 378A -> 379

Disulfide detection

Disulfides are found geometrically: the script measures the distance between every pair of cysteine SG atoms (across all chains) and renames both residues CYS → CYX when the distance is below the cutoff (default 2.5 Å). It does not rely on CONECT records, so it works even when connectivity is missing or incomplete.

For every cysteine it also prints the nearest-neighbour SG distance, so borderline cases are easy to inspect:

  disulfides (auto, SG-SG < 2.5 A):
    CYX A/32 -- CYX A/56  (2.05 A)
  all CYS SG nearest-neighbour distances:
    CYS A/32    nearest A/56    2.05 A -> CYX
    CYS A/297   nearest A/338   3.44 A

To override the geometry (e.g. to force a bridge, or to name residues yourself):

# add explicit residues to the auto-detected set
python3 prepare_protein.py --cyx "297 338"

# use ONLY the listed residues, chain-qualified, no geometry
python3 prepare_protein.py --cyx "A32 A56" --no-auto-disulfide

# loosen/tighten the distance threshold
python3 prepare_protein.py --ss-cutoff 2.3

HIS protonation

Histidine protonation is read directly from the reduce-added hydrogens present in the input:

Hydrogen present	Assigned residue name
HD1 only (on ND1)	HID
HE2 only (on NE2)	HIE
both HD1 and HE2	HIP

Because these markers are on nitrogen, they survive the "del HC" step, so the assignment is unaffected by hydrogen removal.

Removing carbon-bound hydrogens ("del HC")

For each hydrogen the script finds its bonded heavy atom (the nearest non-hydrogen atom in the same residue) and deletes the hydrogen when that atom is a carbon. Polar hydrogens on N/O/S — including the backbone amide H, the HIS HD1/HE2, and cap N–H atoms — are kept. This is geometry-based rather than name-based, so it correctly handles ambiguous names (such as HD1 on ND1 versus a carbon) and the renamed NMA CM. Pass --keep-carbon-h to skip this step.

Chains and hetero groups

• Renumbering, disulfide detection, caps, the H1→H rename, and HIS assignment all work per chain, so multi-chain receptors are supported.
• --relabel-chains renames chains to sequential letters in order of appearance (e.g. A, C → A, B).
• --unify-chain-id A gives every chain the same ID in the final rec.pdb only, while rec.crg.pdb keeps the proper chain names — this is the trick for getting Blastermaster to renumber a 2-chain receptor.
• Waters and other hetero groups are dropped by default. Keep specific ions or cofactors with --keep-het "CA ZN"; kept groups are placed after the protein and are not reordered or renumbered. They are carried through every stage, so they appear in rec_noHC.pdb, rec.crg.pdb and the final rec.pdb (minus their hydrogens).
• No receptor file contains HETATM records. Kept groups are rewritten as ATOM, because Blastermaster deletes HETATM lines — a group left as HETATM would silently vanish from the docking receptor. Only the residue name still identifies them (e.g. ATOM ... HEM A 401). Be aware that a calcium ion is then an ATOM record whose atom name is CA, like a protein alpha carbon; anything downstream that keys on the atom name alone, rather than the residue name or element column, can misread it.

Non-canonical residues

Non-canonical residues come in two kinds, and the script separates them automatically — there is no option to set and no list of residue names to maintain.

Modified residues that are part of the sequence (selenomethionine MSE, phosphoserine SEP, phosphothreonine TPO, …) are usually written as HETATM even though they are covalently built into the backbone. The script finds them geometrically: a HETATM residue is part of a chain when it has backbone N/CA/C atoms and its backbone N (or C) is within 2.0 Å of the C (or N) of a residue already in the chain. The search repeats until nothing new is found, so runs of consecutive modified residues are picked up. Such residues are converted to ATOM and then treated exactly like any other residue: kept in their place in the sequence, included in the numbering check, reordered and renumbered, and carried through to rec.crg.pdb and rec.pdb. The run log names them:

  modified residues bonded into a chain (HETATM -> ATOM, treated as protein): A/169 MSE

Free cofactors and ligands (HEM, metal ions, waters, and also a free amino acid sitting in a binding site) have no peptide bond, so the script does not treat them as part of the sequence: they are dropped by default, or kept with --keep-het and placed after the protein. Kept groups are written as ATOM records so that Blastermaster does not delete them. The log shows both decisions, so nothing disappears silently:

  modified residues bonded into a chain (HETATM -> ATOM, treated as protein): A/169 MSE
  kept HETATM, written as ATOM: CA, HEM
  dropped HETATM: GLY, HOH

In short: a modified residue needs no flag (it is part of the protein and is always kept), while a free cofactor or ion needs --keep-het (whether it belongs in the docking receptor is your decision).

Note that a cofactor such as HEM still needs charges and radii in the Blastermaster parameter tables; keeping it in the PDB does not by itself make it parameterised.

Examples

# Default run in the current directory
python3 prepare_protein.py

# Keep a catalytic calcium and a zinc
python3 prepare_protein.py --keep-het "CA ZN"

# Keep a haem cofactor (modified residues such as MSE need no flag)
python3 prepare_protein.py --keep-het "HEM"

# Two-chain receptor: relabel A,C -> A,B and unify IDs in rec.pdb
python3 prepare_protein.py --relabel-chains --unify-chain-id A

# Force a specific disulfide set, no geometric detection
python3 prepare_protein.py --cyx "A32 A56 A297 A338" --no-auto-disulfide

# Write outputs to a separate directory
python3 prepare_protein.py --output-dir prep_4x93

Verifying the output

The run log is the first sanity check. Confirm that:

• the numbering check passes (or reports the fix it applied);
• the detected disulfides match what you expect;
• the HIS states match your inspection in Chimera;
• the final lines read NMA CA->CM applied: True and NMA still CM: True;
• any modified residue you expect was reported as bonded into the chain, and any ion or cofactor you asked for appears under kept HETATM, written as ATOM rather than under dropped HETATM.

Then open rec.pdb and working/rec.crg.pdb in a viewer and check the termini, caps, disulfides and any kept ions before running Blastermaster.

See also

• Blastermaster
• DOCK 3.8
• Preparing a receptor for DOCK