Running ChemSTEP
last update: oct 8 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].
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.
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.
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.
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.
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
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.
1. Source ChemSTEP virtual environment on Wynton
source /wynton/group/bks/work/shared/kholland/chemstep_0.3.1.5/bin/activate
2. Copy seed-set SDI file into your docking directory
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 .
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.
cp /wynton/group/bks/work/shared/kholland/chemstep_13B/130K_seeds.wynton.sdi . cp /wynton/group/bks/work/shared/kholland/chemstep_13B/1.3M_seeds.wynton.sdi . cp /wynton/group/bks/work/shared/kholland/chemstep_13B/13M_seeds.wynton.sdi . cp /wynton/group/bks/work/shared/kholland/chemstep_13B/26M_seeds.wynton.sdi .
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.
export MOLECULES_DIR_TO_BIND=[outermost folder containing the molecules to dock, just your base docking directory]
export DOCKFILES=[path to your dockfiles]
export INPUT_FOLDER=[the folder containing your .sdi file(s)]
export OUTPUT_FOLDER=[where you want the output ]
/wynton/group/bks/work/bwhall61/needs_github/super_dock3r.sh
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 (do this in a screen!!!):
cp /wynton/group/bks/work/shared/kholland/chemstep_v0p0/get_scores.py .
python get_scores.py 0
cp /wynton/group/bks/work/shared/chemstep_13B_scripts_tutorial/get_scores.py .
python get_scores.py 0
This script expects a directory named "output*" within the CWD. If your output from docking follows different organizational conventions, vim into get_scores.py and change the path. Successful output of the script 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.
4. Convert scores and molecule IDS into NumPY arrays for ChemSTEP recognition. Requires ChemSTEP venv
cp /wynton/group/bks/work/shared/kholland/chemstep_v0p0/convert_scores_to_npy.py .
python convert_scores_to_npy.py 0
cp /wynton/group/bks/work/shared/chemstep_13B_scripts_tutorial/convert_scores_to_npy.py .
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.
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.
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 .
cp /wynton/group/bks/work/shared/chemstep_13B_scripts_tutorial/params.txt .
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.
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.
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).
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.**
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.
max_n_rounds: if prospectively running ChemSTEP, as outlined in this wiki, no need to worry about this parameter.
There should be no need to edit run_chemstep.py or the SGE wrapper script. At this time, we strongly suggest running ChemSTEP as a job array with 16 or 32 CPU slots requested. The number of cores must be specified when calling CSAlgo (in run_chemstep.py) and the SGE wrapper(launch_chemstep_as_job.sh).
7. Run ChemSTEP with the following command:
qsub launch_chemstep_as_job.sh
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.
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.
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.
Prioritized molecules should be built, docked, and their scores fed back into ChemSTEP. More detailed instructions below:
8. Build prioritized molecules (DOCK 3.8). /taken from docs.docking.org
source /wynton/group/bks/soft/DOCK-3.8.5/env.sh
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
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:
find /path/to/your/building_output -name "bundle.db2.tgz" -type f > round_{}.wynton.sdi
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).
9. Set up for iterative rounds of ChemSTEP
cp /wynton/group/bks/work/shared/kholland/chemstep_v0p0/run_chemstep_iterative.py .
cp /wynton/group/bks/work/shared/kholland/chemstep_v0p0/launch_chemstep_iterative.sh .
cp /wynton/group/bks/work/shared/chemstep_13B_scripts_tutorial/run_chemstep_iterative.py .
cp /wynton/group/bks/work/shared/chemstep_13B_scripts_tutorial/launch_chemstep_iterative.sh .
10. Run ChemSTEP
qsub launch_chemstep_as_job.sh [round number]
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.
Helper scripts: a running list of scripts we have been using for analysis and visualization.
to get SMILES for beacons:
python /wynton/group/bks/work/bwhall61/needs_github/get_beacon_smi.py --beacon_df_path /path/to/your/chemstep/output/complete_info/beacons.df --n_workers 6 --outfile beacon_smiles.smi
- anecdotally true.