sanjayk741/gcoupler repository overview

⁠Interactome Prediction using Gcoupler

⁠Introduction

Gcoupler leverages an integrative approach combining de novo ligand design, statistical methods, and Graph Neural Networks for rationally predicting high-affinity ligands. Gcoupler offers an efficient and comparatively faster route to explore endogenous allosteric sites of GPCRs, including the GPCR-Gα interface.

Note: Gcoupler, in principle can work flawlessly on almost all protein types; however, we only validated its functionality on orthostatic and allosteric cavities of GPCRs.

⁠How to use Gcoupler?

This is the repo of the official Docker image for Gcoupler

Pull the Gcoupler image from Docker Hub by running the following command in your terminal:

$ docker pull sanjayk741/gcoupler

Verify the new image has been created using the docker images command.

$ docker images

To access the terminal of a Docker image, you can use the docker run command with the -it option.

$ docker run -it <image-name> /bin/bash

Replace <image-name> with the name or ID of the Docker image of Gcoupler.

Find the ID of the currently running container for input and output.

$ docker ps -a

To start the container again and access its terminal.

$ docker start <container-ID>
$ docker exec -it <container-ID> bash

⁠Running Gcoupler

You can either run Gcoupler by python command line mode directly

$ python
>>> import Synthesizer
>>>

⁠OR

By Jupyter Notebook in following 2 step

STEP-I

Run the following command to get the IP address of your running container.

hostname -I

STEP-II

Launch Jupyter Notebook by running the following command:

$ jupyter-notebook --ip=<ip-address> --allow-root

Replace <ip-address> with the IP address that you noted down in STEP-I

This should start the Jupyter Notebook server and display a URL in the terminal that you can use to access the Jupyter Notebook interface.

⁠Default paths inside the container

>>> import Synthesizer as sz
>>> sz.Set_paths(LigBuilder_path='/LG/LigBuilderV3.0/',libiomp5_file='/opt/conda/envs/Gcoupler/lib/',Output_dir='/root/CDir/')

The option Output_dir can be changed in case user decides to make other directories inside the container.

⁠Input/Output

⁠Input

Find the ID of the currently running container, just executed using the docker ps -a command.

$ docker ps -a

To write a file to the container, use the docker cp command to copy the file from the host to the container.

$ docker cp file container_id:/root/CDir/

This command will copy the file file (pdb/tsv) from the host's current directory to the Gcoupler container with ID container_id at the CDir/ directory inside the container.

⁠Output

Find the ID of the currently running container, just executed using the docker ps -a command.

$ docker ps -a

To write a file from the container, use the docker cp command to copy the file from the container to the host.

$ docker cp container_id:/root/CDir/file-name .

This command will copy the file file (csv/pdf)from the Gcoupler container with ID container_id under the CDir/ directory inside the container to the host's current directory.

⁠Gcoupler Pipeline

Gcoupler supports four distinct modules:

1. Synthesizer
2. Authenticator
3. Generator
4. BioRanker

⁠Synthesizer

To identify the putative cavities on the protein surface with the potential to be an active or an allosteric site and perform de novo drug synthesis

>>> import Synthesizer as sz

Set paths for the installed third-party software and default output folder to collect the intermediate result files and plots

>>> sz.Set_paths(LigBuilder_path='path to LigBuilderV3.0/',libiomp5_file='path to libiomp5.so file containing folder/',Output_dir='path to default folder/')

Note: to find the folder containing libiomp5.h file

$ find / -name libiomp5.so

To submit the query protein file of interest in PDB format

>>> sz.input_structure(path='pre-set default Output folder/',pdb='path to pdbfile.pdb')
Cavity will output 16 cavity file(s)

The output shows the total number of cavities predicted (in this case, 16), which can be visualized by its integer identifier.

Users can view the cavity scores for all the detected cavities (Optional)

>>> sz.cavity_info(path='pre-set default Output folder/')

Users can also visualize a single cavity of choice in interactive mode (Optional)

>>> cavity=4 #To select cavity number 4 to view
>>> sz.cavity_view(path='pre-set default Output folder/',CvID=cavity)

Users can either directly choose a cavity number for the ligand synthesis.

>>> cavity=4 #To select cavity number 4 as the target cavity
>>> sz.compound_synthesis(path='pre-set default Output folder/',CavID=cavity)

Or the user can opt for cavity detection by submitting residue of interest in a TSV (Tab-separated) file.

>>> sz.cavity_detect(path='pre-set default Output folder/',res_list='Residue_list.tsv')

Residue_list.tsv (The first column contains the Single-letter amino acid codes, and the second column their respective positions)

$ head -5 Residue_list.tsv
E       305
T       306
I       310
Y       316
V       466

Note: No header should be provided in the TSV file

⁠Optional

User can specify the number of compounds to synthesize (Default: 500)

>>> cavity=4 #To select cavity number 4 as the target cavity
>>> lcount=800 #To synthesize 800 compounds 
>>> sz.compound_synthesis(path='pre-set default Output folder/',CavID=cavity, ligand_count=lcount)

OR

>>> lcount=800 #To synthesize 800 compounds 
>>> sz.cavity_detect(path='pre-set default Output folder/',res_list='Residue_list.tsv', ligand_count=lcount)

⁠Output folder

The output folder will contain the following files at the end of the successful execution of the Synthesizer module

Note: Synthesizer module with a single cavity for generation of ~500 synthetic binders (HABs/LABs) takes approximately one day

⁠Authenticator

To segregate the synthetic compounds into binary classes based on their actual interaction at the molecular level

>>> import Authenticator as au

To calculate the interaction (binding energy) of individual synthetic compounds with the target cavity (in which they are synthesized)

>>> au.synthetic_interaction(path='pre-set default Output folder/')

Additional arguments:

1. method: Statistical test to use for binding energy cutoff estimation

1. p_val: Significance cutoff for the statistical test (Default: 0.05)
2. plot: Plot to visualize the distribution of HABs and LABs at each qualified cutoff

Example

>>> au.synthetic_interaction(path='pre-set default Output folder/',method='KS-test',p_val=0.05,plot='Density')

To classify synthetic compounds into binary classes of HAB & LAB based on binding energy cutoff (Default: -7)

>>> cutoff = -9  #user decided binding energy cutoff for synthetic compound binary classification 
>>> au.synthetic_classify(path='pre-set default Output folder/',cf=cutoff)

In case the user want to opt for decoys as a negative class against HABs

>>> au.synthetic_decoys(path='pre-set default Output folder/')

Additional arguments:

Decoy_data.csv

$ head -5 Decoy_data.csv
SMILES,Annotation
O(CC1CO1)c2cccc3ccccc23,HAB
Nc1nc(N)nc2ncc(CN(C)c3ccc(cc3)C(=O)NC(CCC(=O)O)C(=O)O)nc12,HAB
N(N=C(C)CCCCCC)c1ccc(cc1[N+](=O)[O-])[N+](=O)[O-],Decoy
C(CCCCCC)CCC[Se]C#N,Decoy

Note: User should use either synthetic_classify or synthetic_decoys function to proceed with not both

⁠Output folder

The output folder will contain the following files at the end of the successful execution of Authenticator module

Note: Authenticator module without decoy generation takes approximately few minutes to complete. And it takes less than 5 minute with around ~300 decoy compound generation

⁠Generator

To build Graph-Neural Network-based classification models for large-scale screening of the user query compounds

>>> import Generator as ge

To pre-process the binary data (HAB & LAB) and test against four base models

1. GCM: GraphConv Model
2. AFP: Attentive FP
3. GCN: Graph Convolution Network
4. GAT: Graph Attention Network

>>> ge.multi_model_test(path='pre-set default Output folder/')

⁠Optional

Users can provide pre-compiled binary data for multi-model testing (limited to this function only)

>>> data='/home/username/cmp.csv' #cmp.csv file containing "SMILES" & "Status" column with SMILES of compounds and 1/0 as class information respectively 
>>> ge.multi_model_test(path='pre-set default Output folder/',fi=data)

Save the base model scoring metrics as Pandas data frame

>>> matrices = ge.multi_model_test(path='pre-set default Output folder/')

To select the best-performing model for the hyperparameter tuning (HPT) with K-Fold cross-validation

>>> ge.MD_kfold(path='pre-set default Output folder/',mdl='GCN',k=5)

Additional arguments:

Users can either opt for Gcoupler predefined hyperparameter grid for the selected model of interest for HPT.

OR

Users can also specify the range of each hyperparameter for the selected model of interest for HPT.

Modifiable Hyperparameter list:

⁠Optional

Users can opt for hyperparameter tuning of any selected model without K-Fold cross-validation (Only for testing different hyperparameter grid range)

>>> ge.model_hpt(path='pre-set default Output folder/',mdl='GCN', params={dictionary}) #dictionary=['parameter name':[parameter grid]]

Note: The above function does not create a model for large-scale screening.

⁠Output folder

The output folder will contain the following files at the end of the successful execution of the Generator module

Note: Generator module with simple hyper parameter tuning and 3Fold cross validation takes around 1hour to complete

⁠Large-scale screening

To predict the binding probability for individual query compounds

>>> import Generator as ge

Prepare a list of canonical SMILES (OpenBabel generated) strings of the query compounds

>>> smiles =  ['ClCC=C', 'C=CCOC(=O)CC(C)C', ...]

Run predictions on the pre-trained Graph-Neural Network model (Model selected for K-Fold cross on the last run)

>>> ge.MD_pred(path='pre-set default Output folder/',smi_list=smiles)

Save the result as Pandas data frame

>>> result = ge.MD_pred(path='pre-set default Output folder/',smi_list=smiles)

Note: Prediction for query compounds with a pre-trained model for a given receptor cavity in Gcoupler takes less than a second for each SMILE

⁠BioRanker

To evaluate and compare chemical compounds against the synthetic high-affinity binders for further screening

>>> import BioRanker as br

For detailed insights into the biological properties that the query compound may mimic

>>> br.analyse(path = 'pre-set default Output folder/',Qdf = smiles_list)   #list of query SMILES

Users can also use the prediction output from the Generator module as input for segregation of molecules

>>> import Generator as ge
>>> import BioRanker as br
>>> smiles_list =  ['c1cc2c(cc1Cl)n1c(c(n2)NC2CCCCC2)nnc1C(F)(F)F', 'c1ccc2c(c1)c(c(c(=O)o2)N)Nc1ccccc1C(=O)O', ...]
>>> result = ge.MD_pred(path='pre-set default Output folder/',smi_list=smiles)
>>> br.analyse(path = 'pre-set default Output folder/',Qdf = result)

The module by default, considers compounds with a probability of 0.9 or above for the processing User can provide their own thresholds using

>>> result
                                                 smiles    Probability 
0          c1cc2c(cc1Cl)n1c(c(n2)NC2CCCCC2)nnc1C(F)(F)F       0.824532 
1              c1ccc2c(c1)c(c(c(=O)o2)N)Nc1ccccc1C(=O)O       0.464561
...                                                 ...            ...
322                  c1cc(c2c(c1)[nH]c(n2)c1ccc(cc1)N)N       0.913462

>>> br.analyse(path = 'pre-set default Output folder/', Qdf = result, threshold=0.8)  

⁠Threshold Optimization

To opt for an optimal probability cutoff (instead of default) for a more precise selection of query compounds

>>> import BioRanker as br

Users can choose the model of interest to optimize the cutoff against

>>> br.threshold_optimize(path='pre-set default Output folder/',mdl='GCN',method='Optimisation method')

Options:

Note: It is recommended that users optimize the cutoff for the GNN model used for the prediction.

⁠Biological properties under BioRanker

BioRanker by default, analyses query molecules for the main properties

Users can also provide their custom list of main and/or sub-properties for analysis by the module

>>> property_list =['Chem', 'Bnd', 'ThrpA']
>>> br.analyse(path = 'pre-set default Output folder/', Qdf = result, Property = property_list)

Users can choose to analyze all the sub-properties in a single run

>>> br.analyse(path = 'pre-set default Output folder/', Qdf = result, Property = 25)

⁠Optional

The module by default, provides the top 15% of the compounds based on the property combination provided, along with a scaled heatmap of the rank score for each compound with respective properties.

The user can change the threshold for top hits

>>> out_df = br.analyse(path = 'pre-set default Output folder/', Qdf = result, top_hit = 20)   #for top 20%
>>> out_df
       QID                                              smiles 
0     Hit1        c1cc2c(cc1Cl)n1c(c(n2)NC2CCCCC2)nnc1C(F)(F)F 
1     Hit2                       c(c(c(=O)o2)N)Nc1ccccc1C(=O)O
...    ...                                                 ...
86   Hit52                  c1cc(c2c(c1)[nH]c(n2)c1ccc(cc1)N)N

The user can export the rank score matrix

>>> out_df = br.analyse(path = 'pre-set default Output folder/', Qdf = result, out_matrix=True)
>>> out_df
       QID        Chem         Bnd       ThrpA                                          smiles 
0     Hit1    0.626817    0.896565   -0.564356    c1cc2c(cc1Cl)n1c(c(n2)NC2CCCCC2)nnc1C(F)(F)F 
1     Hit2    0.735635    0.027846   -0.346736                   c(c(c(=O)o2)N)Nc1ccccc1C(=O)O
...    ...         ...         ...         ...                                             ...
52   Hit52   -0.178472    0.434677    0.246789              c1cc(c2c(c1)[nH]c(n2)c1ccc(cc1)N)N