Running ChemSTEP: Difference between revisions

From DISI
Jump to navigation Jump to search
No edit summary
No edit summary
 
(19 intermediate revisions by the same user not shown)
Line 1: Line 1:
last update: oct 8 2025 katie. current ver = 0.3.1.5 (automatic resubmission of failed SGE jobs).  
last update: dec 16 2025 katie. current ver = 0.3.1.5 (automatic resubmission of failed SGE jobs).  


ChemSTEP (Chemical Space Traversal and Exploration Procedure) is an open-source, transparent acceleration algorithm for molecular docking capable of dealing with virtual libraries of several trillion compounds. This wiki page is a guide for BKS lab members to run ChemSTEP on Wynton HPC, using the current version of InifiSee XReal library (1.1T) or 13.3B library from Enamine REAL. '''for detailed instructions on the 13B space, see 'Running ChemSTEP on the 13B space' wiki page'''. For more general use directions, please refer to [ChemSTEP Read-the-Docs].  
ChemSTEP (Chemical Space Traversal and Exploration Procedure) is an open-source, transparent acceleration algorithm for molecular docking capable of dealing with virtual libraries of several trillion compounds. This wiki page is a guide for BKS lab members to run ChemSTEP on Wynton HPC, using a drug-like subset of ZINC (22B) or 13.3B library from Enamine REAL. '''for detailed instructions on the 13B space, see 'Running ChemSTEP on the 13B space' wiki page'''. For more general use directions, please refer to [ChemSTEP Read-the-Docs].  


At a high-level, ChemSTEP is an iterative process that identifies molecules from the larger virtual library to prioritize for docking. First, we identify a random sample of the total library (termed "seed set", round zero) and dock those molecules to the target of interest. From this seed set, we can calculate total-library pProp values (-log rank percentages) and the number of "virtual hits" in the total library (high-scoring molecules). ChemSTEP will identify a set of maximally diverse molecules that score above the desired pProp threshold ("beacons") from the seed set. These beacons guide prioritization, where molecules chosen and output by ChemSTEP are close in chemical space to the beacons. Prioritized molecules are then built, docked, and returned to ChemSTEP for a second round of prioritization. This process is iterated until you reach desired virtual hit recovery, or you are no longer recovering virtual hits.  
At a high-level, ChemSTEP is an iterative process that identifies molecules from the larger virtual library to prioritize for docking. First, we identify a random sample of the total library (termed "seed set", round zero) and dock those molecules to the target of interest. From this seed set, we can calculate total-library pProp values (-log rank percentages) and the number of "virtual hits" in the total library (high-scoring molecules). ChemSTEP will identify a set of maximally diverse molecules that score above the desired pProp threshold ("beacons") from the seed set. These beacons guide prioritization, where molecules chosen and output by ChemSTEP are close in chemical space to the beacons. Prioritized molecules are then built, docked, and returned to ChemSTEP for a second round of prioritization. This process is iterated until you reach desired virtual hit recovery, or you are no longer recovering virtual hits.  




Running the ChemSTEP algorithm on Wynton takes place in loosely three steps: (1) submission of the main ChemSTEP job, (2) ChainingLog sub-job array, and (3) gathering of SMILES for prioritization.
= Running ChemSTEP (Auto DOCK and Build) =


'''1. Main ChemSTEP job:''' reads input files (DOCK scores and line-matched indices NumPy arrays) and assigns beacons. Beacons are a set of N maximally diverse molecules that score above the pProp threshold specified in the parameter file. For the first round of ChemSTEP, the best-scoring molecule from the seed set is assigned as the first beacon. Subsequent beacons would be molecules that score well (above pProp thresh), and are as diverse as possible from already assigned beacons. Greedy max-diversity selection.  
ChemSTEP is configured to run on Wynton with libraries of '''13B''' and '''22B'''. This page covers the full workflow for running ChemSTEP with automatic submission of docking and building jobs.


During assignment, ECFP4 Tanimoto distances between the assigned beacon and all potential beacons are calculated. The second beacon is the molecule with the maximum Td from beacon 1. Beacon 3 will be the molecule that has the highest Td from beacon 1 and beacon 2. This continues until the number of beacons specified in the parameter file is reached. In iterative rounds of ChemSTEP, the beacon selection is based on diversity from ALL PREVIOUSLY assigned beacons, including those in earlier rounds.
__TOC__


== 1. Source Environment ==


'''2. ChainingLog sub-jobs:''' once beacons are assigned, ChemSTEP will launch a series of sub-jobs that calculate the Tanimoto distances of every molecule remaining in the library to the assigned beacons. Calculated distances to the NEAREST beacon (minimum-minimum Td) are updated in the mintddistrib_*.npy files within the /output directory.
<pre>
source /wynton/group/bks/work/shared/kholland/chemstep_auto_v02/bin/activate
</pre>


== 2. Dock the Seed Set ==


'''3. Gathering of SMILES for prioritization:''' When all Td calculations are completed, the algo will select N number of molecules for prioritization, the number of which is specified by the user in the parameter file. Molecules that are ''closest in chemical space'' to beacons are prioritized. I.e. if round size = 1 million, the 1 million molecules with the smallest min-Td to ANY beacon (current or previous rounds) are prioritized. chemstep_algo.log will provide a "max-minTd" for each round, which is the Tanimoto distance of the most dissimilar molecule prioritized to any one beacon per round. The prioritized molecules are output in /output/complete_info as smi_round_{}.smi
Copy the <code>.sdi</code> file for the library you want to use:


What the user need: DOCKFILES, directories for (1) docking (2) building and (3) running ChemSTEP. '''All iterative rounds of ChemSTEP should be run in the SAME directory.'''
{| class="wikitable"
! Library !! Path
|-
| 13B || <code>/wynton/group/bks/work/shared/kholland/chemstep_auto_v02/scripts/libraries/13B/13M_seeds.sdi</code>
|-
| 22B || <code>/wynton/group/bks/work/shared/kholland/chemstep_auto_v02/scripts/libraries/22B/22M_seeds.sdi</code>
|}


'''1. Source ChemSTEP virtual environment on Wynton'''
Then dock the seed set. See the '''Large-Scale Docking (LSD)''' directions.
    source /wynton/group/bks/work/shared/kholland/chemstep_0.3.1.5/bin/activate


== 3. Gather Scores for the Seed Set ==


'''2. Copy seed-set SDI file into your docking directory'''
Once docking is complete, run the following from the directory '''one level above''' your docking output (<code>MOLECULES_DIR_TO_BIND</code>).
This directory should already contain your dockfiles, with INDOCK parameters set to your liking. In this step, we are copying in a split database index file (SDI) containing paths to bundles of db2 files. This seed set contains 100 million molecules sampled randomly from the total virtual library, currently 1.1 trillion molecules.
    cp /wynton/group/bks/work/shared/kholland/chemstep_v0p0/XR_seed_set_v0p0.wynton.sdi .


'''22B library:'''
<pre>
python /wynton/group/bks/work/shared/kholland/chemstep_auto_v02/scripts/get_scores.py 0
</pre>


if you are interested in running on a library of 13.2B compounds from Enamine REAL, there are several seed set sizes available: 130k, 1.3M, 13M, and 26M.
'''13B library:'''
    cp /wynton/group/bks/work/shared/kholland/chemstep_13B/130K_seeds.wynton.sdi .
<pre>
    cp /wynton/group/bks/work/shared/kholland/chemstep_13B/1.3M_seeds.wynton.sdi .
python /wynton/group/bks/work/shared/kholland/chemstep_auto_v02/scripts/get_scores_13B.py 0 MOL
    cp /wynton/group/bks/work/shared/kholland/chemstep_13B/13M_seeds.wynton.sdi .
</pre>
    cp /wynton/group/bks/work/shared/kholland/chemstep_13B/26M_seeds.wynton.sdi .


{{Note|You must specify the molecule ID prefix for the 13B library (<code>MOL</code>).}}


'''3. Dock seed set to your receptor of interest using DOCK 3.8''' /docking directions taken from docs.docking.org. This is meant to be done as you would do a normal LSD.
Verify that <code>scores_round_0.txt</code> was correctly written:
    export MOLECULES_DIR_TO_BIND=[outermost folder containing the molecules to dock, just your base docking directory]
<pre>
    export DOCKFILES=[path to your dockfiles]
wc -l scores_round_0.txt
    export INPUT_FOLDER=[the folder containing your .sdi file(s)]
</pre>
    export OUTPUT_FOLDER=[where you want the output ]


    /wynton/group/bks/work/bwhall61/needs_github/super_dock3r.sh
== 4. Convert Scores to .npy Files ==


Wait for docking to complete. Next, you must extract all molecule IDs and corresponding DOCK scores from above. To do so, run the following commands in the base docking directory containing your docking files and output folder, while logged into a dev node (I suggest in a screen):
Convert scores to ChemSTEP-readable <code>.npy</code> files:
    cp /wynton/group/bks/work/shared/kholland/chemstep_v0p0/get_scores.py .
    python get_scores.py 0


<pre>
python /wynton/group/bks/work/shared/kholland/chemstep_auto_v02/scripts/convert_scores_to_npy.py 0 <mol_id_prefix>
</pre>


    cp /wynton/group/bks/work/shared/chemstep_13B_scripts_tutorial/get_scores.py
The <code>mol_id_prefix</code> should match the library:
    python get_scores.py 0


{| class="wikitable"
! Library !! Prefix
|-
| 22B / 72B || <code>CSLB</code>
|-
| 13B || <code>MOL</code>
|}


This script expects a directory named "output*" within the CWD. If your output from docking follows different naming conventions, vim into get_scores.py and change the path. The output will be a file named "scores_round_0.txt". For iterative rounds, pass increasing numbers into the command line. i.e. When docking the first round of prioritized molecules, pass [1] for scores_round_1.txt.
== 5. Set Up the ChemSTEP Run Directory ==


'''4. Convert scores and molecule IDS into NumPY arrays for ChemSTEP recognition. Requires ChemSTEP venv'''
Create and enter a new run directory, then copy in the necessary files:
    cp /wynton/group/bks/work/shared/kholland/chemstep_v0p0/convert_scores_to_npy.py .
    python convert_scores_to_npy.py 0


<pre>
mkdir chemstep_run
cd chemstep_run
chemstep-run-new
</pre>


    cp /wynton/group/bks/work/shared/chemstep_13B_scripts_tutorial/convert_scores_to_npy.py
This will populate the directory with <code>params.txt</code>, <code>run_chemstep.py</code>, and <code>launch_chemstep_as_job.sh</code>.
    python convert_scores_to_npy.py 0


For python speed purposes, ChemSTEP requires that DOCK scores and IDs be give in the form of a NumPY array. In this step, we are taking our readable scores file and converting them for ChemSTEP recognition. This script expects a txt file with Molecule IDS and DOCK scores. Use the same round number you used in step 3. As run above, this will output two files named scores_round_0.npy and indices_round_0.npy that contain line-matched indices (determined from Mol ID) and their respective docking scores.
=== Optional: Integrated IFP ===


'''5. Enter into or make a directory to run ChemSTEP in. Copy in necessary files for initiating ChemSTEP: params.txt, run_chemstep.py and launch_chemstep_as_job.sh. Copy in your score and indices numpy files as well, which should be in your docking directory.'''
If running with integrated IFP for beacon selection, also run:
    cp /wynton/group/bks/work/shared/kholland/chemstep_v0p0/params.txt .
    cp /wynton/group/bks/work/shared/kholland/chemstep_v0p0/run_chemstep.py .
    cp /wynton/group/bks/work/shared/kholland/chemstep_v0p0/launch_chemstep_as_job.sh .


<pre>
chemstep-run-ifp
</pre>


    cp /wynton/group/bks/work/shared/chemstep_13B_scripts_tutorial/params.txt .
This copies in the additional files <code>ifp_acceptance_criteria.txt</code> and <code>interactions.txt</code>.
    cp /wynton/group/bks/work/shared/chemstep_13B_scripts_tutorial/run_chemstep.py .
    cp /wynton/group/bks/work/shared/chemstep_13B_scripts_tutorial/launch_chemstep_as_job.sh .


'''6. Edit params.txt file''' This ONLY needs to be edited for the initial round of ChemSTEP. The parameters outlined here will be carried through all rounds of ChemSTEP.
== 6. Edit params.txt ==
    seed_indices_file: /absoulte/path/to/your/indices_round_0.npy
    seed_scores_file: /absolute/path/to/your/scores_round_0.npy
    hit_pprop: 5
    n_docked_per_round: 10000000
    max_beacons: 150
    max_n_rounds: 250


Within params.txt, add the absolute paths to the ChemSTEP-readable score and indices NumPY arrays. The rest of the values within the params.txt file are left to the discretion of the user, with some considerations below.  
Add the absolute paths to the ChemSTEP-readable score and indices <code>.npy</code> arrays generated in Step 4.


'''pProp:''' this value will define what is considered a "virtual hit" in your campaign. pProp is defined as the -log(rank %) of a molecule within the total library score-distribution. For example, a pProp of 4 in the 1.1T XReal space is equivalent the top 0.01% of the library, corresponding to the top 110 million molecules. pProp 5 = 0.001% = 11 million virtual hits, etc. From the seed set, ChemSTEP will estimate a DOCK score value in kcal/mol that associated with the lower-limit of your desired pProp zone. Any molecules that score better than the threshold is considered a virtual hit in your campaign. Be mindful of seed set size when choosing desired pProp. Our suggestion is that the size of the seed set should be at least 10^(pProp +2).
<pre>
seed_indices_file: /path/to/your/indices_round_0.npy
seed_scores_file:  /path/to/your/scores_round_0.npy
hit_pprop:          5.5
n_docked_per_round: 2000000
bundle_size:        1000
max_beacons:        100
max_n_rounds:      250
</pre>


'''n_docked_per_round:''' this number is the desired molecules for prioritization per round. When choosing this value, note that this number of molecules must be built and docked in between each round of ChemSTEP. Prioritizing many molecules will slow building and docking speeds, and coupled with few beacons (below) may lead to decreased diversity. Prioritizing too few molecules may result in slower virtual hit recovery. Round size does not significantly impact the algorithm running time.** 
=== Parameter Reference ===


'''max_beacons:''' this is the number of diverse, well-scoring molecules ChemSTEP will use to guide prioritization per round. All molecules that score better than the assigned pProp threshold may be considered for beacons. By default, ChemSTEP chooses each set of beacons to be as maximally diverse as possible. Choosing too many beacons may result in decreased diversity between beacons, but too few beacons could hinder space exploration. If not enough molecules score above your pProp threshold, ChemSTEP may assign fewer beacons than the assigned value.  
{| class="wikitable"
! Parameter !! Description
|-
| <code>hit_pprop</code> || Defines a "virtual hit." pProp = −log(rank%) within the total library score distribution. E.g., pProp 4 in 13B space ≈ top 0.01% (~1.3M molecules); pProp 5 ≈ 0.001% (~132K). The seed set should contain at least 10<sup>(pProp+2)</sup> molecules.
|-
| <code>n_docked_per_round</code> || Number of molecules prioritized per round. All must be built and docked between rounds. Too many slows throughput and may reduce diversity; too few slows virtual hit recovery.
|-
| <code>max_beacons</code> || Diverse, well-scoring molecules used to guide prioritization. All molecules above the pProp threshold are candidates. Too many reduces inter-beacon diversity; too few hinders space exploration. Fewer beacons than specified may be assigned if insufficient molecules clear the threshold.
|-
| <code>bundle_size</code> || In auto docking mode, number of molecules submitted as a single build job.
|-
| <code>max_n_rounds</code> || No adjustment needed when running ChemSTEP prospectively as described here.
|}


'''max_n_rounds:''' if prospectively running ChemSTEP, as outlined in this wiki, no need to worry about this parameter.
== 7. Edit run_chemstep.py ==


'''Note:''' All paths must be '''absolute paths'''.


There should be no need to edit run_chemstep.py or the SGE wrapper script. ''If using another virtual library, be sure to update the path in run_chemstep.py to point to your FP library.'' At this time, we strongly suggest running ChemSTEP as a job array with 16 CPU slots requested. The number of cores should be specified when calling CSAlgo (in run_chemstep.py) and the SGE wrapper (launch_chemstep_as_job.sh).
=== Required Settings ===


Set <code>lib_path</code> to the library pickle for your library:


'''7. Run ChemSTEP''' with the following command:
{| class="wikitable"
    qsub launch_chemstep_as_job.sh
! Library !! Path
|-
| 13B || <code>/wynton/group/bks/work/shared/kholland/chemstep_auto_v02/scripts/libraries/13B/boltz_fplib.pickle</code>
|-
| 22B || <code>/wynton/group/bks/work/shared/kholland/chemstep_auto_v02/scripts/libraries/22B/22B_fplib.pickle</code>
|}


<pre>
lib_path = '/full/path/to/library.pickle'
</pre>


This will launch the main ChemSTEP job to the SGE scheduler. Check your job status with the 'qstat' command. The algorithm will read the score and indices files provided, calculate pProp (for round 1) and assign beacons. As this job runs, it will launch a subsequent job array. These sub-jobs are the Chaining step of ChemSTEP, where each parallel worker calculates the Tanimoto distances of 100 million molecules from the virtual library to the beacons. These distances are then written to the files within the /output directory. When these jobs complete, the main job will read through all Tanimoto distances and pull molecules for prioritization.
Set <code>dockfiles_path</code>:


Running ChemSTEP successfully will result in the output of a SMILES file within output/complete_info/. You will also get (1) chemstep_algo.log and (2) chemstep_submission.log in the working directory after job submission. chemstep_algo.log contains the running output of the algorithm in a readable format, with timestamps, including the DOCK score associated with your desired pProp. chemstep_submission.log contains the output of ChemSTEP, as well as any information output by the SGE submission system (python errors, cluster issues, tracebacks if failed). These files will be updated with any information with iterative rounds of ChemSTEP run in the same directory. '''Visually inspect these files after each round of ChemSTEP to ensure the algorithm has picked your desired number of beacons and things ran smoothly.'''
<pre>
dockfiles_path="/full/path/to/dockfiles"
</pre>


'''Troubleshooting:''' if no jobs are running and there is no SMI file, check the chemstep_submission.log first. Any traceback error or SGE error should give you some idea of why ChemSTEP failed. Some errors may be due to SGE or Wynton issues. If directed, look at a few .out files within /output/jobs/ . If a ChemSTEP run fails on your FIRST run, delete the output and log files, fix what needs to be fixed, and rerun. Errors in subsequent rounds during the chaining step can potentially corrupt chaining files within the /output directory that are needed for prioritization. At the very least, you will definitely have duplicates of beacons and information written to the log files. At this time, if your run fails during iterative rounds, it's best to start from the beginning.
=== Optional: minTD Exclusion Zone ===


Prioritized molecules should be built, docked, and their scores fed back into ChemSTEP. More detailed instructions below:  
Molecules will not be prioritized from within a specified Tanimoto distance of beacons. Comment in the relevant lines and update the value. Consider also setting <code>enforce_n_docked_per_round = True</code> when using this option:


'''8. Build prioritized molecules (DOCK 3.8).''' /taken from docs.docking.org
<pre>
      source /wynton/group/bks/soft/DOCK-3.8.5/env.sh
min_td_search=0.5,
enforce_n_docked_per_round=True,
</pre>


      python /wynton/group/bks/soft/DOCK-3.8.5/DOCK3.8/zinc22-3d/submit/submit_building_docker.py --output_folder building_output --bundle_size 1000 --minutes_per_mol 5 --skip_name_check --scheduler sge --container_software apptainer --container_path_or_name /wynton/group/bks/soft/DOCK-3.8.5/building_pipeline.sif smi_round_{}.smi
=== Optional: Integrated IFP ===


When building has completed, you must write an SDI file with the complete paths to each built bundle and dock. you can do this with:
Only selects beacons that satisfy user-defined interaction criteria. Comment in the following lines and update the paths to the necessary files (copied in Step 5 if you ran <code>chemstep-run-ifp</code>):
    find /path/to/your/building_output -name "bundle.db2.tgz" -type f > round_{}.wynton.sdi


<pre>
use_IFP=True,
ifp_pdb_path='/full/path/to/rec.crg.pdb',
interactions_file='/full/path/to/interactions.txt',
ifp_acceptance_criteria_file='/full/path/to/ifp_acceptance_criteria.txt',
</pre>


'''Be sure to change the INDOCK file to save only poses that meet your score pProp score threshold calculated by ChemSTEP!''' Retrieve docking scores as convert to NumPy arrays as outlined above, updating the round number when running get_scores.py and convert_scores_to_npy.py. Copy new score and indices files into the directory you ran ChemSTEP in. If you are following along as a tutorial, you should have scores_round_1.npy and indices_round_1.npy from the previous step (from FIRST round of ChemSTEP prioritization).  
'''<code>interactions.txt</code>''' — one interaction per line, comma-separated. Format: <code>interaction_type, residue_name_and_number</code>. Example:


<pre>
Hydrogen bond, GLY19
Ionic, ASP149
</pre>


'''9. Set up for iterative rounds of ChemSTEP'''
Supported interaction types include: Proximal, Hydrogen bond, Ionic, Cation-pi, Hydrophobic, Halogen bond, and others. See LUNA and IFP documentation for the full list.


      cp /wynton/group/bks/work/shared/kholland/chemstep_v0p0/run_chemstep_iterative.py .
'''<code>ifp_acceptance_criteria.txt</code>''' — defines the number of unsatisfied donors/acceptors/specific interactions required for a molecule to pass IFP and be considered for beacon selection. Example:
      cp /wynton/group/bks/work/shared/kholland/chemstep_v0p0/launch_chemstep_iterative.sh .


<pre>
#_donors
#_acceptors
#_unstatisfied_donors == 0
#_unstatisfied_acceptors <= 4
Ionic/ASP-149 > 0
</pre>


      cp /wynton/group/bks/work/shared/chemstep_13B_scripts_tutorial/run_chemstep_iterative.py .
=== Example: AmpC on 22B with minTD=0.50, No IFP ===
      cp /wynton/group/bks/work/shared/chemstep_13B_scripts_tutorial/launch_chemstep_iterative.sh .


'''10. Run ChemSTEP'''
<pre>
      qsub launch_chemstep_as_job.sh [round number]
lib_path = '/wynton/group/bks/work/shared/kholland/chemstep_auto_v02/scripts/libraries/22B/22B_fplib.pickle'
lib = load_library_from_pickle(lib_path)
algo = CSAlgo(lib, 'params.txt', 'output', 16, verbose=True,
    scheduler='sge', smi_id_prefix='CSLB',
    python_exec="/wynton/group/bks/work/shared/kholland/chemstep_auto_v02/bin/python",
    dockfiles_path="/wynton/group/bks/work/kholland/chemstep_ampc_22B/seed_docking/dockfiles",
    min_td_search=0.5,
    enforce_n_docked_per_round=True,
    #use_IFP=True,
    #ifp_pdb_path='/path/to/your/reference/rec.crg.pdb',
    #interactions_file='/path/to/your/interactions.txt',
    #ifp_acceptance_criteria_file='/path/to/your/ifp_acceptance_criteria.txt',
    docking_method="auto", track_beacon_orig=True)
</pre>


For the first iterative round, the round number is [2], and should increase by one for every subsequent round of ChemSTEP. The output will be smi_round_****.smi file. Repeat steps 8-10 for as many rounds as needed. The performance is reported in outputy/complete_info/run_summary.df, which contains the number of beacons selected, the number of molecules docked, the number of hits found, the distance threshold for the selected molecules to dock, and the last added beacon's distance to all previous beacons.
== 8. Launch the Job ==


Submit the main ChemSTEP job:


** anecdotally true.
<pre>
qsub launch_chemstep_as_job.sh
</pre>
 
== 9. Monitor Job Status ==
 
Check job status with <code>qstat</code>. The main job will run for up to '''2 weeks''' given no errors. ChemSTEP will launch search, building, and docking jobs in successive rounds.
 
'''Note:''' If any building or docking subjobs hang, the main job will not proceed until those are canceled or finished. Monitor job statuses regularly and occasionally verify that docking output files (<code>scores_round_*.txt</code>) are being populated.
 
== 10. View Beacon SMILES and IDs ==
 
From the ChemSTEP running directory, run the following in a '''screen session on a dev node''':
 
<pre>
python /wynton/group/bks/work/shared/kholland/chemstep_auto_v02/scripts/get_beacon_smiles.py /path/to/library/pickle chemstep_algo.log
</pre>
 
Use the library pickle path from [[#7. Edit run_chemstep.py|Step 7]].
 
== 11. Get Poses After Docking ==
 
Make a list of <code>test.mol2.gz.0</code> files from docking:
 
<pre>
find /round_*_docking/bundle_paths -maxdepth 2 -name "test.mol2.gz.0" > docked_poses.txt
</pre>
 
Then extract top poses:
 
<pre>
python /wynton/group/bks/work/bwhall61/for_beau/top_poses.py \
    -t <pProp_threshold> \
    -s <num_poses_per_file> \
    -dock_results_path docked_poses.txt
</pre>

Latest revision as of 16:53, 16 March 2026

last update: dec 16 2025 katie. current ver = 0.3.1.5 (automatic resubmission of failed SGE jobs).

ChemSTEP (Chemical Space Traversal and Exploration Procedure) is an open-source, transparent acceleration algorithm for molecular docking capable of dealing with virtual libraries of several trillion compounds. This wiki page is a guide for BKS lab members to run ChemSTEP on Wynton HPC, using a drug-like subset of ZINC (22B) or 13.3B library from Enamine REAL. for detailed instructions on the 13B space, see 'Running ChemSTEP on the 13B space' wiki page. For more general use directions, please refer to [ChemSTEP Read-the-Docs].

At a high-level, ChemSTEP is an iterative process that identifies molecules from the larger virtual library to prioritize for docking. First, we identify a random sample of the total library (termed "seed set", round zero) and dock those molecules to the target of interest. From this seed set, we can calculate total-library pProp values (-log rank percentages) and the number of "virtual hits" in the total library (high-scoring molecules). ChemSTEP will identify a set of maximally diverse molecules that score above the desired pProp threshold ("beacons") from the seed set. These beacons guide prioritization, where molecules chosen and output by ChemSTEP are close in chemical space to the beacons. Prioritized molecules are then built, docked, and returned to ChemSTEP for a second round of prioritization. This process is iterated until you reach desired virtual hit recovery, or you are no longer recovering virtual hits.


Running ChemSTEP (Auto DOCK and Build)

ChemSTEP is configured to run on Wynton with libraries of 13B and 22B. This page covers the full workflow for running ChemSTEP with automatic submission of docking and building jobs.

1. Source Environment

source /wynton/group/bks/work/shared/kholland/chemstep_auto_v02/bin/activate

2. Dock the Seed Set

Copy the .sdi file for the library you want to use:

Library Path
13B /wynton/group/bks/work/shared/kholland/chemstep_auto_v02/scripts/libraries/13B/13M_seeds.sdi
22B /wynton/group/bks/work/shared/kholland/chemstep_auto_v02/scripts/libraries/22B/22M_seeds.sdi

Then dock the seed set. See the Large-Scale Docking (LSD) directions.

3. Gather Scores for the Seed Set

Once docking is complete, run the following from the directory one level above your docking output (MOLECULES_DIR_TO_BIND).

22B library: 

python /wynton/group/bks/work/shared/kholland/chemstep_auto_v02/scripts/get_scores.py 0

13B library: 

python /wynton/group/bks/work/shared/kholland/chemstep_auto_v02/scripts/get_scores_13B.py 0 MOL

Template:Note

Verify that scores_round_0.txt was correctly written: 

wc -l scores_round_0.txt

4. Convert Scores to .npy Files

Convert scores to ChemSTEP-readable .npy files: 

python /wynton/group/bks/work/shared/kholland/chemstep_auto_v02/scripts/convert_scores_to_npy.py 0 <mol_id_prefix>

The mol_id_prefix should match the library:

Library Prefix
22B / 72B CSLB
13B MOL

5. Set Up the ChemSTEP Run Directory

Create and enter a new run directory, then copy in the necessary files: 

mkdir chemstep_run
cd chemstep_run
chemstep-run-new

This will populate the directory with params.txt, run_chemstep.py, and launch_chemstep_as_job.sh.

Optional: Integrated IFP

If running with integrated IFP for beacon selection, also run: 

chemstep-run-ifp

This copies in the additional files ifp_acceptance_criteria.txt and interactions.txt.

6. Edit params.txt

Add the absolute paths to the ChemSTEP-readable score and indices .npy arrays generated in Step 4. 

seed_indices_file:  /path/to/your/indices_round_0.npy
seed_scores_file:   /path/to/your/scores_round_0.npy
hit_pprop:          5.5
n_docked_per_round: 2000000
bundle_size:        1000
max_beacons:        100
max_n_rounds:       250

Parameter Reference

Parameter Description
hit_pprop Defines a "virtual hit." pProp = −log(rank%) within the total library score distribution. E.g., pProp 4 in 13B space ≈ top 0.01% (~1.3M molecules); pProp 5 ≈ 0.001% (~132K). The seed set should contain at least 10(pProp+2) molecules.
n_docked_per_round Number of molecules prioritized per round. All must be built and docked between rounds. Too many slows throughput and may reduce diversity; too few slows virtual hit recovery.
max_beacons Diverse, well-scoring molecules used to guide prioritization. All molecules above the pProp threshold are candidates. Too many reduces inter-beacon diversity; too few hinders space exploration. Fewer beacons than specified may be assigned if insufficient molecules clear the threshold.
bundle_size In auto docking mode, number of molecules submitted as a single build job.
max_n_rounds No adjustment needed when running ChemSTEP prospectively as described here.

7. Edit run_chemstep.py

Note: All paths must be absolute paths.

Required Settings

Set lib_path to the library pickle for your library:

Library Path
13B /wynton/group/bks/work/shared/kholland/chemstep_auto_v02/scripts/libraries/13B/boltz_fplib.pickle
22B /wynton/group/bks/work/shared/kholland/chemstep_auto_v02/scripts/libraries/22B/22B_fplib.pickle 
lib_path = '/full/path/to/library.pickle'

Set dockfiles_path: 

dockfiles_path="/full/path/to/dockfiles"

Optional: minTD Exclusion Zone

Molecules will not be prioritized from within a specified Tanimoto distance of beacons. Comment in the relevant lines and update the value. Consider also setting enforce_n_docked_per_round = True when using this option: 

min_td_search=0.5,
enforce_n_docked_per_round=True,

Optional: Integrated IFP

Only selects beacons that satisfy user-defined interaction criteria. Comment in the following lines and update the paths to the necessary files (copied in Step 5 if you ran chemstep-run-ifp): 

use_IFP=True,
ifp_pdb_path='/full/path/to/rec.crg.pdb',
interactions_file='/full/path/to/interactions.txt',
ifp_acceptance_criteria_file='/full/path/to/ifp_acceptance_criteria.txt',

interactions.txt — one interaction per line, comma-separated. Format: interaction_type, residue_name_and_number. Example: 

Hydrogen bond, GLY19
Ionic, ASP149

Supported interaction types include: Proximal, Hydrogen bond, Ionic, Cation-pi, Hydrophobic, Halogen bond, and others. See LUNA and IFP documentation for the full list.

ifp_acceptance_criteria.txt — defines the number of unsatisfied donors/acceptors/specific interactions required for a molecule to pass IFP and be considered for beacon selection. Example: 

#_donors
#_acceptors
#_unstatisfied_donors == 0
#_unstatisfied_acceptors <= 4
Ionic/ASP-149 > 0

Example: AmpC on 22B with minTD=0.50, No IFP

lib_path = '/wynton/group/bks/work/shared/kholland/chemstep_auto_v02/scripts/libraries/22B/22B_fplib.pickle'
lib = load_library_from_pickle(lib_path)
algo = CSAlgo(lib, 'params.txt', 'output', 16, verbose=True,
    scheduler='sge', smi_id_prefix='CSLB',
    python_exec="/wynton/group/bks/work/shared/kholland/chemstep_auto_v02/bin/python",
    dockfiles_path="/wynton/group/bks/work/kholland/chemstep_ampc_22B/seed_docking/dockfiles",
    min_td_search=0.5,
    enforce_n_docked_per_round=True,
    #use_IFP=True,
    #ifp_pdb_path='/path/to/your/reference/rec.crg.pdb',
    #interactions_file='/path/to/your/interactions.txt',
    #ifp_acceptance_criteria_file='/path/to/your/ifp_acceptance_criteria.txt',
    docking_method="auto", track_beacon_orig=True)

8. Launch the Job

Submit the main ChemSTEP job: 

qsub launch_chemstep_as_job.sh

9. Monitor Job Status

Check job status with qstat. The main job will run for up to 2 weeks given no errors. ChemSTEP will launch search, building, and docking jobs in successive rounds.

Note: If any building or docking subjobs hang, the main job will not proceed until those are canceled or finished. Monitor job statuses regularly and occasionally verify that docking output files (scores_round_*.txt) are being populated.

10. View Beacon SMILES and IDs

From the ChemSTEP running directory, run the following in a screen session on a dev node: 

python /wynton/group/bks/work/shared/kholland/chemstep_auto_v02/scripts/get_beacon_smiles.py /path/to/library/pickle chemstep_algo.log

Use the library pickle path from Step 7.

11. Get Poses After Docking

Make a list of test.mol2.gz.0 files from docking: 

find /round_*_docking/bundle_paths -maxdepth 2 -name "test.mol2.gz.0" > docked_poses.txt

Then extract top poses: 

python /wynton/group/bks/work/bwhall61/for_beau/top_poses.py \
    -t <pProp_threshold> \
    -s <num_poses_per_file> \
    -dock_results_path docked_poses.txt
Retrieved from "http://wiki.docking.org/index.php?title=Running_ChemSTEP&oldid=17120"