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 is renumbered to the next integer (379), and every following residue in that chain is shifted up by one so nothing collides.
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. Backbone amide H. Renames the backbone amide hydrogen (H1) of the residue following an ACE cap to H.
7. Disulfides. Detects disulfide-bonded cysteines and renames them CYS → CYX. See Disulfide detection.
8. Writes rec_noHC.pdb.
9. 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/.
10. 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)

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". Waters and unlisted hetero groups are dropped. ACE/NMA caps are always kept and converted to ATOM.
--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: 378A becomes 379, and every following residue is shifted up by one to avoid a collision (a cascade). 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 / following residues shifted):
    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 HETATM

• 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 HETATM groups are placed after the protein and are not reordered or renumbered.

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"

# 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.

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