Differences

This shows you the differences between two versions of the page.

--- running_alphafold_at_scale [2026/03/10 11:19] – [Running AlphaFold on Perun] 172.20.54.200
+++ running_alphafold_at_scale [2026/03/11 11:51] (current) – 172.20.80.5
@@ Line 34: / Line 34: @@
-==== Running AlphaFold on Perun ====
+===== Running AlphaFold on Perun =====
 Currently we have AlphaFold3 installed on Perun.
-=== Get permission to use AlphaFold3 ===
+==== Get permission to use AlphaFold3 ====
 Before you can use it, you need to get a ~1GB file that contains the "model parameters". Each new user needs to get this file separately. You can get the file by submitting a form, essentially declaring that you are not going to share those parameters with the public. You must use your institutional e-mail address (i.e. your @dal.ca address) and also provide your gmail if you have it. Each model parameter contains a unique identifier linked to a user, so each model parameter file is unique. You can fill in and submit the form [[https://docs.google.com/forms/d/e/1FAIpQLSfWZAgo1aYk0O4MuAXZj8xRQ8DafeFJnldNOnh_13qAx2ceZw/viewform|here]].
@@ Line 44: / Line 44: @@
 Once you get approved, normally you would download the file and place it in your filesystem at the appropriate place so you can start using AlphaFold. For Perun however, you don't need to download that file, you only need to let Balagopal know that you have been approved and he will add you on a list, so that you can start using AlphaFold on Perun.
-=== Searching the protein sequence databases ===
+==== Searching the protein sequence databases ====
 As explained above, AlphaFold does two things. First it compares the query sequence to sequence databases, and then it actually predicts the 3D structure. The sequence database searching is the computational limiting step and takes a lot more time to run. The structure module is relatively quick.
+=== The sequence databases ===
 All the sequence databases ("Big Fantastic Database" (BFD)*, MGnify protein database (MGY), UniProt, SwissProt, PDB, Rfam, RNACentral, NT non-coding RNA) together take up a lot of space (currently ~432 GB on Perun). It is therefore recommended to have at least 1TB of free disk space. On Perun, the databases are stored at ''/db1/alphafold3/public_databases''.
@@ Line 52: / Line 54: @@
 *BFD was created by clustering 2.5 billion protein sequences from Uniprot/TrEMBL+Swissprot, Metaclust and Soil Reference Catalog Marine Eukaryotic Reference Catalog assembled by Plass.
-Since there is quite a bit of reading and writing going on in the sequence search phase, it is highly recommended to use a SSD drive. ``/db1/`` is indeed an SSD drive.
+Since there is quite a bit of reading and writing going on in the sequence search phase, it is highly recommended to use a SSD drive. ''/db1/'' is indeed an SSD drive.
+It is most efficient to run the sequence search phase (''run_alphafold.py --norun_inference'') using CPUs, and the structure prediction (''run_alphafold.py --norun_data_pipeline'') using a GPU.
+=== Preparing AlphaFold3 input query files ===
+AlphaFold3 is a little awkward with how it wants its input protein sequences formatted. Instead of a simple FASTA file, it requires a JSON file. Thankfully I found a relatively simple way of converting the FASTA into a JSON file using and editing slightly a python script I found online. It is called ''fasta2alphafoldjson.py'' and it is available on the [[https://github.com/RogerLab/gospel_of_andrew|Gospel of Andrew]] GitHub page. Run it simply as follows:
+<code>
+./fasta2alphafoldjson.py [input.fasta] > [output.json]
+</code>
+NOTE: Currently the script is hardcoded to make JSON entries of the 'alphafold3' dialect. This dialect is appropriate if you want to feed AlphaFold a single query protein. In my setup, this works because I will break up a multi query job into many single query jobs. If however, you want to run a multiquery job as is, you'll need the JSON entries to be of the 'alphafoldserver' dialect. So you'll have to edit the code slightly by replacing ''"dialect": "alphafold3"'' to ''"dialect": "alphafoldserver"''.
+If your ''[input.fasta]'' was a FASTA file with many protein sequences, your ''[output.json]'' will also contain all those protein sequences. Because we want to split our multi-query job into many individual-query jobs, we'll have to split up that JSON into many individual JSON files as well. We can do so as follows:
+<code>
+mkdir -p batch_split_jsons
+jq -c '.[]' [output.json] | awk '{print $0 > "batch_split_jsons/input_" NR ".json"}'
+</code>
+''jq'' splits up the JSON file, and then ''awk'' will give each new single JSON file a number. Each new JSON file now contains a single protein sequence and has been assigned a number.
+=== Submitting AlphaFold as an array job ===
+If you want to run AlphaFold on many proteins, it may be pragmatic to submit these sequence database searches as an **array job** on Perun. An array-job is a single job, that in turn schedules the submission of many other jobs. Here the idea is to submit a single general "AlphaFold sequence search job" which then schedules and submits all the individual search jobs for each query protein separately.
+<code>
+#!/bin/bash
+#$ -S /bin/bash
+#$ -cwd
+#$ -pe threaded 10
+#$ -q 768G-batch,2T-batch
+# $SGE_TASK_ID returns a number, which is its current task number
+# for example if you do
+# qsub -t 1-1000
+# then each job in the array will have a unique SGE_TASK_ID, varying from 1 to 1000
+# here, we are using that ID to get the corresponding file
+SINGLE_JSON_FILE="batch_split_jsons/input_${SGE_TASK_ID}.json"
+# setting up environment
+# first ensure that we are using the right python installation for run_alphafold.py
+# this is needed to make the 'source activate' below work
+export PATH=/scratch5/software/miniconda3/bin:$PATH
+# then append path to include the run_alphafold.py script
+# NOTE: run_alphafold.py was edited slightly to include the #!/usr/bin/env python shebang at the top
+export PATH=/scratch5/software/alphafold3:$PATH
+# activate environment
+source activate alphafold3
+# dir with all the sequence databases
+DB_DIR='/db1/alphafold3/public_databases'
+# where to find the model parameter file
+MODEL_DIR='/scratch5/db/alphafold3-model'
+# input protein or proteins
+#
+# if you want to submit a single protein, use a JSON of the 'alphafold3' dialect
+# if you want to submit multiple protein queries, you can use a JSON of the 'alphafoldserver' dialect
+# or you can specify an INPUT_DIR where each query is a single JSON file of the 'alphafold3' dialect
+# sequence search output directory
+SEQ_OUTPUT_DIR='batch_sequence_search_out'
+# threads to use
+THREADS=10
+mkdir -p $SEQ_OUTPUT_DIR
+# do sequence search with fast CPUs
+run_alphafold.py \
+    --norun_inference \
+    --json_path $SINGLE_JSON_FILE \
+    --output_dir $SEQ_OUTPUT_DIR \
+    --db_dir $DB_DIR \
+    --nhmmer_n_cpu $THREADS \
+    --jackhmmer_n_cpu $THREADS
+</code>
+What is this magic ''$SGE_TASK_ID'' variable??
+<code>
+SINGLE_JSON_FILE="batch_split_jsons/input_${SGE_TASK_ID}.json"
+</code>
+''$SGE_TASK_ID'' is a special variable that is useful when you submit an array job. Array jobs are submitted by passing the ''-t'' flag to ''qsub''. For example, ''qsub -t 1-1000'' tells qsub to submit an array job of a 1000 individual jobs. When let's say, the 52nd individual job is submitted by the array job, the ''$SGE_TASK_ID'' special variable gets assigned the number 52. So, in the example above, we are assigning the variable ''SINGLE_JSON_FILE'' to ''batch_split_jsons/input_52.json''. That file is then finally passed on to ''run_alphafold.py --norun_inference'' at the bottom of the script.
+So, to submit this part of AlphaFold as an array job, do:
+<code>
+qsub -t 1-1000 -tc 10 -p -128 run_alphafold_cpu.arrayjob.sh
+</code>
+  * Adjust the 1000 to the number of proteins that you wished to be analyzed.
+  * The ''-tc'' means it will submit a maximum of 10 individual jobs at any one time, to not overwhelm Perun
+  * The ''-p -128'' gives the individual jobs a lower priority. This ensures that other peoples jobs, which won't take as long, have a higher priority and thus don't have to wait for your whole array job to finish
+Once AlphaFold starts running, it should generate for each query protein a output directory with the name of the protein, under ''batch_sequence_search_out''. For example, ''batch_sequence_search_out/pyv62_000500''. In that directory, you will find an output JSON file, for example ''pyv62_000500_data.json''. This file looks something like:
+<code>
+{
+  "dialect": "alphafold3",
+  "version": 2,
+  "name": "PYV62_000500",
+  "sequences": [
+    {
+      "protein": {
+        "id": "A",
+        "sequence": "MAEVHRKLVSQVFPQTSIAMPLAGCGVQNVLRSQSAAVASVTKAAARTFSDMVTIDGNTAASHVAFHMSDLSCLFPITPATPMGENVDAWAAQGKKNIYGVVPGVYQLNSEAGASGSLHGALLAGALAASYTSSQGLLLMPPNMLKC>
+        "modifications": [],
+        "unpairedMsa": ">query\nMAEVHRKLVSQVFPQTSIAMPLAGCGVQNVLRSQSAAVASVTKAAARTFSDMVTIDGNTAASHVAFHMSDLSCLFPITPATPMGENVDAWAAQGKKNIYGVVPGVYQLNSEAGASGSLHGALLAGALAASYTSSQG>
+        "pairedMsa": ">query\nMAEVHRKLVSQVFPQTSIAMPLAGCGVQNVLRSQSAAVASVTKAAARTFSDMVTIDGNTAASHVAFHMSDLSCLFPITPATPMGENVDAWAAQGKKNIYGVVPGVYQLNSEAGASGSLHGALLAGALAASYTSSQGLL>
+        "templates": [
+          {
+            "mmcif": "data_6CIN\n#\n_entry.id 6CIN\n#\nloop_\n_chem_comp.formula\n_chem_comp.formula_weight\n_chem_comp.id\n_chem_comp.mon_nstd_flag\n_chem_comp.name\n_>
+            "queryIndices": [52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86,>
+            "templateIndices": [3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38,>
+          },
+          {
+            "mmcif": "data_1B0P\n#\n_entry.id 1B0P\n#\nloop_\n_chem_comp.formula\n_chem_comp.formula_weight\n_chem_comp.id\n_chem_comp.mon_nstd_flag\n_chem_comp.name\n_>
+            "queryIndices": [49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83,>
+            "templateIndices": [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 3>
+          },
+          {
+            "mmcif": "data_5C4I\n#\n_entry.id 5C4I\n#\nloop_\n_chem_comp.formula\n_chem_comp.formula_weight\n_chem_comp.id\n_chem_comp.mon_nstd_flag\n_chem_comp.name\n_>
+            "queryIndices": [51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85,>
+            "templateIndices": [3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39,>
+          },
+          {
+            "mmcif": "data_5C4I\n#\n_entry.id 5C4I\n#\nloop_\n_chem_comp.formula\n_chem_comp.formula_weight\n_chem_comp.id\n_chem_comp.mon_nstd_flag\n_chem_comp.name\n_>
+            "queryIndices": [1003, 1004, 1005, 1006, 1007, 1008, 1009, 1010, 1011, 1012, 1013, 1014, 1015, 1016, 1017, 1018, 1019, 1020, 1021, 1022, 1023, 1024, 1025, 1>
+            "templateIndices": [90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 100, 101, 102, 103, 104, 105, 106, 107, 108, 109, 110, 111, 112, 113, 114, 115, 116, 117, 118, 1>
+          }
+        ]
+</code>
+NOTE that this view of the file is truncated with >. The lines are much longer in the actual file. It would be impossible to show the entire file here in the wiki. These files can get quite big. The one I pasted above was 75 MB. I've seen up to 327 MB.
+As you can see it contains the MSA (although I don't yet understand what is meant with paired and unpaired MSA) and a bunch of Templates, which I'm not sure either on what those mean. In any case, it is this JSON file that will be the input for the Structure Module, i.e. the second phase of AlphaFold.
+=== AlphaFold uses jackhmmer under the hood ===
+''run_alphafold.py --norun_inference'' invokes **jackhmmer** of the HMMER suite of programs to search all the different sequence databases. **jackhmmer** is akin to PSI-BLAST, meaning it uses its search results to build an HMM profile, which is then used to search the databases again, which then updates the HMM profile, etc etc. However, if I understand correctly, jackhmmer is run with ''-N1'', meaning it actually does not build an HMM profile. It is therefore equivalent to **phmmer**, which simply compares a protein sequence to a protein sequence database and thus akin to BLASTP.
+Some technical details useful to know if you run into any errors:
+  * Per protein query, multiple independent jackhmmer processes are launched. One for each distinct protein sequence database
+  * The protein query sequence gets copied to the local disk of the node that is executing the jackhmmer process. These local disks are mounted on ''/tmp1/''. These disks are not massive space wise
+  * The multiple sequence alignment that AlphaFold needs to infer its co-evolving sites is also stored on these local disks (in Stockholm format). Sometimes these MSAs can grow extremely large and the local disk can run out of space.
+  * The sequence databases are located on the /db1 SSD. This means there is a lot of back-and-forth traffic over the network between the /db1 and the local disk
+  * jackhmmer does NOT load the entire sequence database into RAM. Instead, it is treated as a stream. The database is also not indexed, and it is also not indexed on the fly. Therefore, **I/O is often the bottleneck.** The CPU deals very quickly with fetched data, but then sits idle waiting for new data to arrive.
+  * I sometimes get out of memory problems, its still unclear to me why exactly
+=== Running the second phase of AlphaFold, actually predicting the 3D structure ===
+We ran the sequence search phase as an array job, and it was using CPUs, since the task of sequence searching has been optimized for CPUs. However, structure prediction is optimized for using GPUs. We have since recently acquired an NVIDIA GPU and it is now part of Perun as perun24.
+This part thankfully doesn't take as much time, but depending on the number of proteins, it can still take a couple of days. Currently it is pretty stupid, as in, its just doing the structure prediction for each query protein one by one in a for loop. There may be a way to parallalize the workload and thereby make it more time efficient, but I'm not sure yet on how to do that.
+''run_alphafold.py --norun_data_pipeline'' takes in the directory that contains the sequence search result JSON as input.
+<code>
+#!/bin/bash
+#$ -S /bin/bash
+#$ -cwd
+#$ -q GPU@perun24
+# setting up environment
+# first ensure that we are using the right python installation for run_alphafold.py
+# this is needed to make the 'source activate' below work
+export PATH=/scratch5/software/miniconda3/bin:$PATH
+# then append path to include the run_alphafold.py script
+# NOTE: run_alphafold.py was edited slightly to include the #!/usr/bin/env python shebang at the top
+export PATH=/scratch5/software/alphafold3:$PATH
+# activate environment
+source activate alphafold3
+# dir with all the sequence databases
+DB_DIR='/scratch5/software/alphafold3/public_databases'
+# where to find the model parameter file
+MODEL_DIR='/scratch5/db/alphafold3-model'
+SEQ_OUTPUT_DIR='batch_sequence_search_out/'
+# structure inference output directory
+STRUC_OUTPUT_DIR='structure_module_out'
+# do structure inference with GPU
+for QUERY in $SEQ_OUTPUT_DIR/*; do
+    run_alphafold.py \
+        --db_dir $DB_DIR \
+        --model_dir $MODEL_DIR \
+        --input_dir=$QUERY \
+        --output_dir=$STRUC_OUTPUT_DIR \
+        --norun_data_pipeline
+done
+</code>
+Simply submit the script with ''qsub run_alphafold_gpu.sh''. It will generate under ''$STRUC_OUTPUT_DIR'' an output directory for each query protein. That directory will contain many files and subdirectories, but the main file we're after is called ''pyv62_000500_model.cif''. This is your final predicted structure!
+=== Evaluating the final output ===
+== [ID]_model.cif ==
+CIF stands for Crystallographic Information File
+You can view the structures using softwares like PyMOL and ChimeraX
+== [ID]_summary_confidences.json ==
+Contains information regarding the expected overall accuracy of the predicted structure:
+The **ptm** or predicted Template Modeling score
+  * Between 0 and 1, with 1 being the perfect score.
+  * This is a measure of accuracy of the entire structure
+The **iptm** or interface pTM score.
+  * Also between 0 and 1. Null, if a monomer.
+  * This is a measure of confidence in all predicted interfaces between subunits in the multimer, or measure of accuracy of relative positions of subunits to one another
+**fraction disordered**
+  * Also between 0 and 1.
+  * What fraction of the structure is disordered?
+**has_clash**
+  * True or False
+  * True if >50% of atoms of a chain "clash"
+**ranking_score**
+  * Ranges from -100 to 1.5 ?
+  * Calculated as follows: 0.8 * ipTM + 2 * pTM + 0.5 * disorder - 100 * has_clash
+  * This calculation is then used to rank the multiple structure predictions
+There are more metrics to discuss, but I don't have the time right now to continue on them