001 从ChEMBL获取化合物数据

说明：涉及下载的代码在运行时，需要花费较多时间

Talktorial T001: This talktorial is part of the TeachOpenCADD pipeline described in the first TeachOpenCADD paper, comprising of talktorials T001-T010.

Aim of this talktorial

In this notebook, we will learn more about the ChEMBL database and how to extract data from ChEMBL, i.e. (compound, activity data) pairs for a target of interest. These data sets can be used for many cheminformatics tasks, such as similarity search, clustering or machine learning.

Our work here will include finding compounds that were tested against a certain target and filtering available bioactivity data.

Contents in Theory

• ChEMBL database
  • ChEMBL web services
  • ChEMBL web resource client
• Compound activity measures
  • IC50 measure
  • pIC50 value

Contents in Practical

Goal: Get a list of compounds with bioactivity data for a given target

• Connect to ChEMBL database
• Get target data (example: EGFR kinase)
  • Fetch and download target data
  • Select target ChEMBL ID
• Get bioactivity data
  • Fetch and download bioactivity data for targets
  • Preprocess and filter bioactivity data
• Get compound data
  • Fetch and download compound data
  • Preprocess and filter compound data
• Output bioactivity-compound data
  • Merge bioactivity and compound data, and add pIC50 values
  • Draw molecules with highest pIC50
  • Freeze bioactivity data to ChEMBL 27
  • Write output file

References

• ChEMBL bioactivity database: Gaulton et al., Nucleic Acids Res. (2017), 45(Database issue), D945–D954
• ChEMBL web services: Davies et al., Nucleic Acids Res. (2015), 43, 612-620
• ChEMBL web-interface
• GitHub ChEMBL web rescource client
• The EBI RDF platform: Jupp et al., Bioinformatics (2014), 30(9), 1338-9
• Info on half maximal inhibitory concentration: (p)IC50
• UniProt website

Theory

ChEMBL database

"ChEMBL is a manually curated database of bioactive molecules with drug-like properties. It brings together chemical, bioactivity and genomic data to aid the translation of genomic information into effective new drugs." (ChEMBL website)

• Open large-scale bioactivity database
• Current data content (as of 09.2020, ChEMBL 27):
  • >1.9 million distinct compounds
  • >16 million activity values
  • Assays are mapped to ~13,000 targets
• Data sources include scientific literature, PubChem bioassays, Drugs for Neglected Diseases Initiative (DNDi), BindingDB database, ...
• ChEMBL data can be accessed via a web-interface, the EBI-RDF platform and the ChEMBL web rescource client

ChEMBL web services

• RESTful web service
• ChEMBL web service version 2.x resource schema:
  

Figure 1: "ChEMBL web service schema diagram. The oval shapes represent ChEMBL web service resources and the line between two resources indicates that they share a common attribute. The arrow direction shows where the primary information about a resource type can be found. A dashed line indicates the relationship between two resources behaves differently. For example, the Image resource provides a graphical-based representation of a Molecule." Figure and description are taken from: Nucleic Acids Res. (2015), 43, 612-620.

ChEMBL web resource client

• Python client library for accessing ChEMBL data
• Handles interaction with the HTTPS protocol
• Lazy evaluation of results -> reduced number of network requests

Compound activity measures

IC50 measure

• Half maximal inhibitory concentration
• Indicates how much of a particular drug or other substance is needed to inhibit a given biological process by half

Figure 2: Visual demonstration of how to derive an IC50 value: (i) Arrange inhibition data on y-axis and log(concentration) on x-axis. (ii) Identify maximum and minimum inhibition. (iii) The IC50 is the concentration at which the curve passes through the 50% inhibition level. Figure Figure "Example IC50 curve demonstrating visually how IC50 is derived" by JesseAlanGordon is licensed under CC BY-SA 3.0. by JesseAlanGordon is licensed under CC BY-SA 3.0.

pIC50 value

• To facilitate the comparison of IC50 values, which have a large value range and are given in different units (M, nM, ...), often pIC50 values are used
• The pIC50 is the negative log of the IC50 value when converted to molar units: $pI{C}_{50}=-lo{g}_{10}(I{C}_{50})$, where $I{C}_{50}$ is specified in units of M
• Higher pIC50 values indicate exponentially greater potency of the drug
• Note that the conversion can be adapted to the respective IC50 unit, e.g. for nM: $pI{C}_{50}=-lo{g}_{10}(I{C}_{50}\ast 10^{-9})=9-lo{g}_{10}(I{C}_{50})$

Other activity measures:

Besides, IC50 and pIC50, other bioactivity measures are used, such as the equilibrium constant KI and the half-maximal effective concentration EC50.

Practical

In the following, we want to download all molecules that have been tested against our target of interest, the epidermal growth factor receptor (EGFR) kinase.

Connect to ChEMBL database

First, the ChEMBL web resource client as well as other Python libraries are imported.

Next, we create resource objects for API access.

Get target data (EGFR kinase)

• Get UniProt ID of the target of interest (EGFR kinase: P00533) from UniProt website
• Use UniProt ID to get target information

Select a different UniProt ID, if you are interested in another target.

Fetch target data from ChEMBL

Download target data from ChEMBL

The results of the query are stored in targets, a QuerySet, i.e. the results are not fetched from ChEMBL until we ask for it (here using pandas.DataFrame.from_records).

More information about the QuerySet datatype:

QuerySets are lazy – the act of creating a QuerySet does not involve any database activity. You can stack filters together all day long, and Django will actually not run the query until the QuerySet is evaluated. (querysets-are-lazy)

Select target (target ChEMBL ID)

After checking the entries, we select the first entry as our target of interest:

CHEMBL203: It is a single protein and represents the human Epidermal growth factor receptor (EGFR, also named erbB1)

Save selected ChEMBL ID.

Get bioactivity data

Now, we want to query bioactivity data for the target of interest.

Fetch bioactivity data for the target from ChEMBL

In this step, we fetch the bioactivity data and filter it to only consider

• human proteins,
• bioactivity type IC50,
• exact measurements (relation '='), and
• binding data (assay type 'B').

Each entry in our bioactivity set holds the following information:

Download bioactivity data from ChEMBL

Finally, we download the QuerySet in the form of a pandas DataFrame.

Note: This step should not take more than 2 minutes, if so try to rerun all cells starting from "Fetch bioactivity data for the target from ChEMBL" or read this message below:

# replace first line in cell below with this other line
bioactivities_df = pd.read_csv(DATA / "EGFR_bioactivities_CHEMBL27.csv.zip", index_col=0)

Note that the first two rows describe the same bioactivity entry; we will remove such artifacts later during the deduplication step. Note also that we have columns for standard_units/units and standard_values/values; in the following, we will use the standardized columns (standardization by ChEMBL), and thus, we drop the other two columns.

If we used the units and values columns, we would need to convert all values with many different units to nM:

Preprocess and filter bioactivity data

1. Convert standard_value's datatype from object to float
2. Delete entries with missing values
3. Keep only entries with standard_unit == nM
4. Delete duplicate molecules
5. Reset DataFrame index
6. Rename columns

1. Convert datatype of "standard_value" from "object" to "float"

The field standard_value holds standardized (here IC50) values. In order to make these values usable in calculations later on, convert values to floats.

2. Delete entries with missing values

Use the parameter inplace=True to drop values in the current DataFrame directly.

3. Keep only entries with "standard_unit == nM"

We only want to keep bioactivity entries in nM, thus we remove all entries with other units.

4. Delete duplicate molecules

Sometimes the same molecule (molecule_chembl_id) has been tested more than once, in this case, we only keep the first one.

Note other choices could be to keep the one with the best value or a mean value of all assay results for the respective compound.

5. Reset "DataFrame" index

Since we deleted some rows, but we want to iterate over the index later, we reset the index to be continuous.

6. Rename columns

We now have a set of 5575 molecule ids with respective IC50 values for our target kinase.

Get compound data

We have a DataFrame containing all molecules tested against EGFR (with the respective measured bioactivity).

Now, we want to get the molecular structures of the molecules that are linked to respective bioactivity ChEMBL IDs.

Fetch compound data from ChEMBL

Let's have a look at the compounds from ChEMBL which we have defined bioactivity data for: We fetch compound ChEMBL IDs and structures for the compounds linked to our filtered bioactivity data.

Download compound data from ChEMBL

Again, we want to export the QuerySet object into a pandas.DataFrame. Given the data volume, this can take some time. For that reason, we will first obtain the list of records through tqdm, so we get a nice progress bar and some ETAs. We can then pass the list of compounds to the DataFrame.

Preprocess and filter compound data

1. Remove entries with missing entries
2. Delete duplicate molecules (by molecule_chembl_id)
3. Get molecules with canonical SMILES

1. Remove entries with missing molecule structure entry

2. Delete duplicate molecules

3. Get molecules with canonical SMILES

So far, we have multiple different molecular structure representations. We only want to keep the canonical SMILES.

Sanity check: Remove all molecules without a canonical SMILES string.

Output (bioactivity-compound) data

Summary of compound and bioactivity data

Merge both datasets

Merge values of interest from bioactivities_df and compounds_df in an output_df based on the compounds' ChEMBL IDs (molecule_chembl_id), keeping the following columns:

• ChEMBL IDs: molecule_chembl_id
• SMILES: smiles
• units: units
• IC50: IC50

Add pIC50 values

As you can see the low IC50 values are difficult to read (values are distributed over multiple scales), which is why we convert the IC50 values to pIC50.

Draw compound data

Let's have a look at our collected data set.

First, we plot the pIC50 value distribution

In the next steps, we add a column for RDKit molecule objects to our DataFrame and look at the structures of the molecules with the highest pIC50 values.

Show the three most active molecules, i.e. molecules with the highest pIC50 values.

Freeze output data to ChEMBL 27

This is a technical step: Usually, we would continue to work with the dataset that we just created (latest dataset).

However, here on the TeachOpenCADD platform, we prefer to freeze the dataset to a certain ChEMBL releases (i.e. ChEMBL 27), so that this talktorial and other talktorials downstream in our CADD pipeline do not change in the future (helping us to maintain the talktorials).

Write output data to file

We want to use this bioactivity-compound dataset in the following talktorials, thus we save the data as csv file. Note that it is advisable to drop the molecule column (which only contains an image of the molecules) when saving the data.

Discussion

In this tutorial, we collected bioactivity data for our target of interest from the ChEMBL database. We filtered the data set in order to only contain molecules with measured IC50 bioactivity values.

Be aware that ChEMBL data originates from various sources. Compound data has been generated in different labs by different people all over the world. Therefore, we have to be cautious with the predictions we make using this data set. It is always important to consider the source of the data and consistency of data production assays when interpreting the results and determining how much confidence we have in our predictions.

In the next tutorials, we will filter our acquired data by Lipinski's rule of five and by unwanted substructures. Another important step would be to clean the molecular data. As this is not shown in any of our talktorials (yet), we would like to refer to the Standardiser library or MolVS as useful tools for this task.