Deep Learning Tutorial

In this tutorial we will demonstrate how to create a basic deep learning model using train, test, and validation sets. For the model a subset of the bladder dataset is used where transcriptomics and drug structures are the inputs and drug response values are the output. Throughout the tutorial we will be using functions included in coderdata. For more infomration on functions used from the package look to the coderdata/API reference.

Import coderdata, pandas, and numpy.

import pandas as pd
import numpy as np
import coderdata as cd

Download and View Datasets

To view the availabe datasets within the coderdata API use the following functions.

cd.list_datasets()

beataml: Beat acute myeloid leukemia (BeatAML) focuses on acute myeloid leukemia tumor data. Data includes drug response, proteomics, and transcriptomics datasets.
bladder: Tumor Evolution and Drug Response in Patient-Derived Organoid Models of Bladder Cancer. Data includes transcriptomics, mutations, copy number, and drug response data.
ccle: Cancer Cell Line Encyclopedia (CCLE).
cptac: The Clinical Proteomic Tumor Analysis Consortium (CPTAC) project is a collaborative network funded by the National Cancer Institute (NCI) focused on improving our understanding of cancer biology through the integration of transcriptomic, proteomic, and genomic data.
ctrpv2: Cancer Therapeutics Response Portal version 2 (CTRPv2)
fimm: Institute for Molecular Medicine Finland (FIMM) dataset.
hcmi: Human Cancer Models Initiative (HCMI) encompasses numerous cancer types and includes cell line, organoid, and tumor data. The models and the data are from the Human Cancer Models Initiative (HCMI) www.cancer.gov/ccg/research/functional-genomics/hcmi; dbGaP accession number phs001486.
mpnst: Malignant Peripheral Nerve Sheath Tumor is a rare, aggressive sarcoma that affects peripheral nerves throughout the body.
nci60: National Cancer Institute 60.
pancreas: Organoid Profiling Identifies Common Responders to Chemotherapy in Pancreatic Cancer. Data includes transcriptomics, mutations, copy number, and drug response data.
prism: Profiling Relative Inhibition Simultaneously in Mixtures.
sarcoma: The landscape of drug sensitivity and resistance in sarcoma. Data includes transcriptomics, mutations, and drug response data.
colorectal: Living organoid biobank of colorectal cancer patients.
liver: Pharmaco-proteogenomic characterization of liver cancer organoids for precision oncology.
novartis: Patient-derived tumor xenografts for drug response prediction.
gcsi: The Genentech Cell Line Screening Initiative (gCSI)
gdscv1: Genomics of Drug Sensitivity in Cancer version 1 (GDSCv1)
gdscv2: Genomics of Drug Sensitivity in Cancer version 2 (GDSCv2)

To download and load the needed dataset(s) use the function below. For this example we will be using bladder.

cd.download("bladder")

Example Output:

Downloaded ‘https://api.figshare.com/v2/file/download/69320275’ to ‘c:\Path\bladder_samples.csv’

Downloaded ‘https://api.figshare.com/v2/file/download/70636806’ to ‘c:\Path\bladder_copy_number.csv.gz’

Downloaded ‘https://api.figshare.com/v2/file/download/70636809’ to ‘c:\Path\bladder_drug_descriptors.tsv.gz’

Downloaded ‘https://api.figshare.com/v2/file/download/70636815’ to ‘c:\Path\bladder_drugs.tsv.gz’

Downloaded ‘https://api.figshare.com/v2/file/download/70636818’ to ‘c:\Path\bladder_experiments.tsv.gz’

Downloaded ‘https://api.figshare.com/v2/file/download/70636821’ to ‘c:\Path\bladder_mutations.csv.gz’

Downloaded ‘https://api.figshare.com/v2/file/download/70636824’ to ‘c:\Path\bladder_transcriptomics.csv.gz’

Downloaded ‘https://api.figshare.com/v2/file/download/70638963’ to ‘c:\Path\genes.csv.gz’

bladder=cd.load("bladder")

Example Output:

Importing raw data …

Importing ‘transcriptomics’ from c:\Path\bladder_transcriptomics.csv.gz … DONE

Importing ‘mutations’ from c:\Path\bladder_mutations.csv.gz … DONE

Importing ‘copy_number’ from c:\Path\bladder_copy_number.csv.gz … DONE

Importing ‘samples’ from c:\Path\bladder_samples.csv … DONE

Importing ‘drugs’ from c:\Path\bladder_drugs.tsv.gz … DONE

Importing ‘drug_descriptors’ from c:\Path\bladder_drug_descriptors.tsv.gz … DONE

Importing ‘experiments’ from c:\Path\bladder_experiments.tsv.gz … DONE

Importing ‘genes’ from c:\Path\genes.csv.gz … DONE

Importing raw data … DONE

The function type() can be used to view what the class of any object is as seen below, bladder is a dataset.

type(bladder)

coderdata.dataset.dataset.Dataset

To view all the datatypes included with a dataset object use the function below. The info() command allows us to see which datatypes are included in the bladder dataset object. We will be using the Transcriptomics, Drugs, and Experiments Data in our example. Note that each datatype has different format options and arguments. To view all options for each type navigate to the usage.

pd.DataFrame({"datatypes":bladder.types()})

	datatypes
0	transcriptomics
1	mutations
2	copy_number
3	samples
4	drugs
5	experiments
6	genes

Examples are shown below of how to view a regular dataset and the specificed datatype.

bladder.drugs

	improve_drug_id	chem_name	pubchem_id	canSMILES	InChIKey	formula	weight
0	SMI_2538	pf-00477736 [who-dd]	135565545	CN1C=C(C=N1)C2=C3C=NNC(=O)C4=C3C(=CC(=C4)NC(=O...	NDEXUOWTGYUVGA-LJQANCHMSA-N	C22H25N7O2	419.50
1	SMI_2538	cyclohexaneacetamide, alpha-amino-n-(5,6-dihyd...	135565545	CN1C=C(C=N1)C2=C3C=NNC(=O)C4=C3C(=CC(=C4)NC(=O...	NDEXUOWTGYUVGA-LJQANCHMSA-N	C22H25N7O2	419.50
2	SMI_2538	nsc800848	135565545	CN1C=C(C=N1)C2=C3C=NNC(=O)C4=C3C(=CC(=C4)NC(=O...	NDEXUOWTGYUVGA-LJQANCHMSA-N	C22H25N7O2	419.50
3	SMI_2538	pf 00477736	135565545	CN1C=C(C=N1)C2=C3C=NNC(=O)C4=C3C(=CC(=C4)NC(=O...	NDEXUOWTGYUVGA-LJQANCHMSA-N	C22H25N7O2	419.50
4	SMI_2538	pf-477736,pf-0044736	135565545	CN1C=C(C=N1)C2=C3C=NNC(=O)C4=C3C(=CC(=C4)NC(=O...	NDEXUOWTGYUVGA-LJQANCHMSA-N	C22H25N7O2	419.50
...	...	...	...	...	...	...	...
5272	SMI_56178	cpdc	5702198	N.N.Cl[Pt]Cl	LXZZYRPGZAFOLE-UHFFFAOYSA-L	Cl2H6N2Pt	300.05
5273	SMI_56178	trans--diamminedichloroplatinum(ii)	5702198	N.N.Cl[Pt]Cl	LXZZYRPGZAFOLE-UHFFFAOYSA-L	Cl2H6N2Pt	300.05
5274	SMI_56178	sr-05000001586	5702198	N.N.Cl[Pt]Cl	LXZZYRPGZAFOLE-UHFFFAOYSA-L	Cl2H6N2Pt	300.05
5275	SMI_56178	cisplatin, united states pharmacopeia (usp) re...	5702198	N.N.Cl[Pt]Cl	LXZZYRPGZAFOLE-UHFFFAOYSA-L	Cl2H6N2Pt	300.05
5276	SMI_56178	sy073091	5702198	N.N.Cl[Pt]Cl	LXZZYRPGZAFOLE-UHFFFAOYSA-L	Cl2H6N2Pt	300.05

5277 rows × 7 columns

bladder.transcriptomics

	entrez_id	transcriptomics	improve_sample_id	source	study
0	7105	69.277129	3370	Gene Expression Omnibus	Lee et al. 2018 Bladder PDOs
1	64102	0.027787	3370	Gene Expression Omnibus	Lee et al. 2018 Bladder PDOs
2	8813	81.403736	3370	Gene Expression Omnibus	Lee et al. 2018 Bladder PDOs
3	57147	20.024448	3370	Gene Expression Omnibus	Lee et al. 2018 Bladder PDOs
4	55732	29.819016	3370	Gene Expression Omnibus	Lee et al. 2018 Bladder PDOs
...	...	...	...	...	...
815970	140032	0.024479	3319	Gene Expression Omnibus	Lee et al. 2018 Bladder PDOs
815971	9782	44.949383	3319	Gene Expression Omnibus	Lee et al. 2018 Bladder PDOs
815972	102724737	0.024479	3319	Gene Expression Omnibus	Lee et al. 2018 Bladder PDOs
815973	101669762	2.203401	3319	Gene Expression Omnibus	Lee et al. 2018 Bladder PDOs
815974	100527949	0.024479	3319	Gene Expression Omnibus	Lee et al. 2018 Bladder PDOs

815975 rows × 5 columns

bladder.experiments

	source	improve_sample_id	improve_drug_id	study	time	time_unit	dose_response_metric	dose_response_value
0	Synapse	3290	SMI_19333	Lee etal 2018 Bladder PDOs	6	days	fit_auc	0.6428
1	Synapse	3290	SMI_29323	Lee etal 2018 Bladder PDOs	6	days	fit_auc	0.8841
2	Synapse	3290	SMI_32114	Lee etal 2018 Bladder PDOs	6	days	fit_auc	0.7491
3	Synapse	3290	SMI_34836	Lee etal 2018 Bladder PDOs	6	days	fit_auc	0.7746
4	Synapse	3290	SMI_41206	Lee etal 2018 Bladder PDOs	6	days	fit_auc	0.8367
...	...	...	...	...	...	...	...	...
32995	Synapse	3423	SMI_8050	Lee etal 2018 Bladder PDOs	6	days	dss	0.1300
32996	Synapse	3423	SMI_8361	Lee etal 2018 Bladder PDOs	6	days	dss	0.0893
32997	Synapse	3423	SMI_8781	Lee etal 2018 Bladder PDOs	6	days	dss	0.1837
32998	Synapse	3423	SMI_8995	Lee etal 2018 Bladder PDOs	6	days	dss	0.2239
32999	Synapse	3423	SMI_9053	Lee etal 2018 Bladder PDOs	6	days	dss	0.0331

33000 rows × 8 columns

Coderdata has a function format() that will reformat any dataset.datatype with the specificed arguments. Each datatype has different arguments specific to its contents and to view each datatypes argument options visit coderdata/usage. Here we will use the bladder.experiment to format. We want to only include auc metric values in our dataframe for this tutorial, however, to view other metric options for a dataset use the pandas unique function.

pd.DataFrame({"dose_response_metric":bladder.experiments.dose_response_metric.unique()})

	dose_response_metric
0	fit_auc
1	fit_ic50
2	fit_ec50
3	fit_r2
4	fit_ec50se
5	fit_einf
6	fit_hs
7	aac
8	auc
9	dss

Set bladder.experiments to only include samples where the dose_response_metric is auc.

bladder.experiments = bladder.format(data_type='experiments', metrics='auc')
bladder.experiments

	source	improve_sample_id	improve_drug_id	study	time	time_unit	dose_response_metric	dose_response_value
26400	Synapse	3290	SMI_19333	Lee etal 2018 Bladder PDOs	6	days	auc	0.7171
26401	Synapse	3290	SMI_29323	Lee etal 2018 Bladder PDOs	6	days	auc	0.8942
26402	Synapse	3290	SMI_32114	Lee etal 2018 Bladder PDOs	6	days	auc	0.7529
26403	Synapse	3290	SMI_34836	Lee etal 2018 Bladder PDOs	6	days	auc	0.7877
26404	Synapse	3290	SMI_41206	Lee etal 2018 Bladder PDOs	6	days	auc	0.8449
...	...	...	...	...	...	...	...	...
29695	Synapse	3423	SMI_8050	Lee etal 2018 Bladder PDOs	6	days	auc	0.7448
29696	Synapse	3423	SMI_8361	Lee etal 2018 Bladder PDOs	6	days	auc	0.6284
29697	Synapse	3423	SMI_8781	Lee etal 2018 Bladder PDOs	6	days	auc	0.8386
29698	Synapse	3423	SMI_8995	Lee etal 2018 Bladder PDOs	6	days	auc	0.9949
29699	Synapse	3423	SMI_9053	Lee etal 2018 Bladder PDOs	6	days	auc	0.4678

3300 rows × 8 columns

Develop Deep Learning Model

We will be using the keras functional API, for information on the keras API visit keras.io. This allows us to use multiple inputs in our model. You may need to install a couple of these dependencies to proceed.

Dependencies that are not included in coderdata:

• tensorflow (Note: tensorflow is only supported up to the python version 3.12)

• matplotlib

• rdkit

import keras
#from numpy import loadtxt
from tensorflow.keras.models import Sequential
from tensorflow.keras.layers import Dense
from sklearn.preprocessing import MinMaxScaler
#from sklearn.model_selection import train_test_split
from rdkit import Chem
from rdkit.Chem import AllChem
from rdkit.Chem.rdFingerprintGenerator import GetMorganGenerator
from keras.models import Model
from keras.layers import Input, Dense, concatenate
from keras.optimizers import Adam
from keras.losses import MeanSquaredError
from keras.metrics import MeanAbsoluteError
from keras import layers
from sklearn.metrics import mean_absolute_error, mean_squared_error, r2_score
import matplotlib.pyplot as plt
from sklearn.linear_model import LinearRegression

Coderdata has 2 split functions a two-way and three-way split that can be used with any dataset. Below are example of how to use both functions but for this tutorial we will be using the three-way splits for our model.

Quick tip: To view all available arguments for any function use the help() function below.

help(bladder.split_train_other)

Help on method split_train_other in module coderdata.dataset.dataset:

split_train_other(split_type: "Literal['mixed-set', 'drug-blind', 'cancer-blind']" = 'mixed-set', ratio: 'tuple[int, int]' = (8, 2), stratify_by: 'Optional[str]' = None, balance: 'bool' = False, random_state: 'Optional[Union[int, RandomState]]' = None, **kwargs: 'dict') -> 'TwoWaySplit' method of coderdata.dataset.dataset.Dataset instance
    Split the dataset into training and another subset (e.g., testing or validation).

    Parameters
    ----------
    split_type : {'mixed-set', 'drug-blind', 'cancer-blind'}, optional
        The type of split to perform, by default 'mixed-set'.
        - `mixed-set`: A random split, disregarding drug or cancer associations.
        - `drug-blind`: Ensures disjoint splits by drug ID.
        - `cancer-blind`: Ensures disjoint splits by sample or cancer association.
    ratio : tuple[int, int], optional
        The ratio of train to other subset sizes, by default (8, 2).
        For instance, (8, 2) translates to an 80%-20% split.
    stratify_by : str, optional
        The column used for stratification, if stratification is needed, by default None.
    balance : bool, optional
        Whether to adjust to balanced splits (equal representation of classes), by default False.
    random_state : int | RandomState | None, optional
        A seed for reproducibility of the random split, by default None.
    **kwargs : dict
        Additional arguments for advanced customization of the split.

    Returns
    -------
    TwoWaySplit
        An object containing the train and other subsets as separate datasets.

    Notes
    -----
    This method is a wrapper around the `split_train_other` utility function and
    ensures that the split configuration is applied to the dataset (self).

    Examples
    --------
    >>> split = dataset.split_train_other(split_type='cancer-blind', ratio=(7,3))
    >>> print(split.train, split.other)

Here is an example of a two-way split that is trained on bladder dataset dataframe. To view the definitions of different split_type and the associated arguments view the API reference. Note that the other split functions will use similar arguments.

two_way_split= bladder.split_train_other(split_type = 'drug-blind',
                                            ratio = [8,2],
                                            random_state = 42,
                                            )

The split function dataframe can be used with any datatype however here we will use the experiments datatype. Now to view both the train and other splits use the call below, you will see that that the train set is ~80% of the original bladder.experiments and the other set is ~20% aligning with the values input for ratio.

two_way_split.train.experiments.head()

	source	improve_sample_id	improve_drug_id	study	time	time_unit	dose_response_metric	dose_response_value
0	Synapse	3290	SMI_19333	Lee etal 2018 Bladder PDOs	6	days	auc	0.7171
1	Synapse	3290	SMI_32114	Lee etal 2018 Bladder PDOs	6	days	auc	0.7529
2	Synapse	3290	SMI_41206	Lee etal 2018 Bladder PDOs	6	days	auc	0.8449
3	Synapse	3290	SMI_43509	Lee etal 2018 Bladder PDOs	6	days	auc	0.6962
4	Synapse	3290	SMI_53355	Lee etal 2018 Bladder PDOs	6	days	auc	0.7154

two_way_split.other.experiments.head()

	source	improve_sample_id	improve_drug_id	study	time	time_unit	dose_response_metric	dose_response_value
0	Synapse	3290	SMI_29323	Lee etal 2018 Bladder PDOs	6	days	auc	0.8942
1	Synapse	3290	SMI_34836	Lee etal 2018 Bladder PDOs	6	days	auc	0.7877
2	Synapse	3291	SMI_29323	Lee etal 2018 Bladder PDOs	6	days	auc	0.8942
3	Synapse	3291	SMI_34836	Lee etal 2018 Bladder PDOs	6	days	auc	0.7877
4	Synapse	3292	SMI_29323	Lee etal 2018 Bladder PDOs	6	days	auc	0.8942

The three-way split uses similar input arguments as the two-way split. However, the three-way split provides a train, test, and validation set.

split = bladder.split_train_test_validate(split_type='mixed-set',
                                             ratio = [8,1,1],
                                             random_state = 42
                                             )

The CoderData package includes its own split function that is intended for only coderdata specific datasets. Here is an overview of what each split is and how it will be used.

Train (80%) : This gathers majority of the data and it allows the model to learn from the patterns and relationships of the data. This step involves iteratively adjusting its parameters to bettter fit the model.

Validate (10%): Gathering generally the smallest portion of data is used to refine the model’s hyperparameters. It evaluates the training performance by comparing it to validation performance based on overfitting or underfitting. Here it will follow each training iteration (epoch) and evaluating the performance based on the validation set.

Test (10%) : Finally to evaluate the final model performance the test data will be used to provide an estimate to how the model will perform on unseen data

Note: The percentage of data pulled for each step may vary but should typically follows this size comparison (Train > Test >= Validate). For a more in depth explanation of how train, test, and validate work visit mlu-explain.github.io.

split.train.experiments

	source	improve_sample_id	improve_drug_id	study	time	time_unit	dose_response_metric	dose_response_value
0	Synapse	3370	SMI_56178	Lee etal 2018 Bladder PDOs	6	days	auc	0.8081
1	Synapse	3373	SMI_55959	Lee etal 2018 Bladder PDOs	6	days	auc	0.9275
2	Synapse	3296	SMI_8361	Lee etal 2018 Bladder PDOs	6	days	auc	0.8972
3	Synapse	3383	SMI_29447	Lee etal 2018 Bladder PDOs	6	days	auc	0.3859
4	Synapse	3381	SMI_20241	Lee etal 2018 Bladder PDOs	6	days	auc	0.5096
...	...	...	...	...	...	...	...	...
2635	Synapse	3345	SMI_31810	Lee etal 2018 Bladder PDOs	6	days	auc	0.6982
2636	Synapse	3346	SMI_10344	Lee etal 2018 Bladder PDOs	6	days	auc	0.8777
2637	Synapse	3349	SMI_29447	Lee etal 2018 Bladder PDOs	6	days	auc	0.4976
2638	Synapse	3340	SMI_46822	Lee etal 2018 Bladder PDOs	6	days	auc	0.5837
2639	Synapse	3421	SMI_41910	Lee etal 2018 Bladder PDOs	6	days	auc	0.8155

2640 rows × 8 columns

split.test.experiments

	source	improve_sample_id	improve_drug_id	study	time	time_unit	dose_response_metric	dose_response_value
0	Synapse	3339	SMI_54374	Lee etal 2018 Bladder PDOs	6	days	auc	0.5098
1	Synapse	3361	SMI_2538	Lee etal 2018 Bladder PDOs	6	days	auc	0.9241
2	Synapse	3380	SMI_56178	Lee etal 2018 Bladder PDOs	6	days	auc	0.8509
3	Synapse	3356	SMI_2538	Lee etal 2018 Bladder PDOs	6	days	auc	0.8234
4	Synapse	3307	SMI_20241	Lee etal 2018 Bladder PDOs	6	days	auc	0.7549
...	...	...	...	...	...	...	...	...
325	Synapse	3301	SMI_9053	Lee etal 2018 Bladder PDOs	6	days	auc	0.4285
326	Synapse	3308	SMI_21256	Lee etal 2018 Bladder PDOs	6	days	auc	0.7320
327	Synapse	3349	SMI_22444	Lee etal 2018 Bladder PDOs	6	days	auc	0.9989
328	Synapse	3365	SMI_8361	Lee etal 2018 Bladder PDOs	6	days	auc	0.8503
329	Synapse	3307	SMI_56178	Lee etal 2018 Bladder PDOs	6	days	auc	0.8176

330 rows × 8 columns

split.validate.experiments

	source	improve_sample_id	improve_drug_id	study	time	time_unit	dose_response_metric	dose_response_value
0	Synapse	3420	SMI_47848	Lee etal 2018 Bladder PDOs	6	days	auc	0.3556
1	Synapse	3372	SMI_55959	Lee etal 2018 Bladder PDOs	6	days	auc	0.9275
2	Synapse	3310	SMI_43509	Lee etal 2018 Bladder PDOs	6	days	auc	0.9430
3	Synapse	3370	SMI_52290	Lee etal 2018 Bladder PDOs	6	days	auc	0.8663
4	Synapse	3306	SMI_32114	Lee etal 2018 Bladder PDOs	6	days	auc	0.4926
...	...	...	...	...	...	...	...	...
325	Synapse	3358	SMI_53355	Lee etal 2018 Bladder PDOs	6	days	auc	0.5340
326	Synapse	3354	SMI_29323	Lee etal 2018 Bladder PDOs	6	days	auc	0.9783
327	Synapse	3347	SMI_38445	Lee etal 2018 Bladder PDOs	6	days	auc	0.9706
328	Synapse	3359	SMI_36872	Lee etal 2018 Bladder PDOs	6	days	auc	0.9726
329	Synapse	3420	SMI_54351	Lee etal 2018 Bladder PDOs	6	days	auc	0.7900

330 rows × 8 columns

The coderdata function save() can be used to keep the split for later use.

Note: The save function can be used with any coderdata Dataset object. However only individual splits can be saved not the full split.

split.train.save(path='bladder_train.pickle')

Preparing the Data

Impute Missing Data

Now we will prepare the transcriptomic and data structure columns to ensure it is formated correclty to be merged with our drug response values and then used within our model.

First impute missing values in the bladder transcriptomics data. Here we take the global mean, however, this may not be the best method for you or needed for your datasets.

#Impute and transcriptomics if there are any based on global mean.
columns_to_fill = ['transcriptomics']

for i in columns_to_fill:
    if i in bladder.transcriptomics.columns[bladder.transcriptomics.isnull().any(axis=0)]:
        bladder.transcriptomics[i] = bladder.transcriptomics[i].fillna(bladder.transcriptomics[i].mean())

Scale Input Data

(Optional) Here use a sklearn function to scale the transcriptomics data in a small range of non-negative numbers. This step can be helpful when there is a wide range of values in the data.

scaler=MinMaxScaler()
bladder.transcriptomics[["transcriptomics"]] = scaler.fit_transform(bladder.transcriptomics[["transcriptomics"]])

Drop any duplicated transcriptomics data.

bladder.transcriptomics[['transcriptomics']].drop_duplicates().head()

	transcriptomics
0	1.572852e-03
1	3.745602e-07
2	1.848217e-03
3	4.544482e-04
4	6.768580e-04

From the bladder.drugs dataframe we are pulling the canSMILES column and using rdkit to convert them into Morgan fingerprints. The resulting fingerprint column will then be used as another input along with transcriptomics for the model.

# Convert SMILES to Morgan fingerprints
def smiles_to_fingerprint(smiles):
    try:
        mol = Chem.MolFromSmiles(smiles)
        if mol is None: # Handle cases where MolFromSmiles fails
                return None
        morgan_gen= GetMorganGenerator(radius=2, fpSize=1024)
        fingerprint = morgan_gen.GetFingerprint(mol)  # Adjust radius and nBits as needed
        fingerprint_array = np.array(fingerprint)
        return fingerprint_array
    except Exception as e:
        print(f"Error processing SMILES '{smiles}': {e}")
        return None


# smiles_to_fingerprint(merged_df.canSMILES[0])
# # Apply the function to your SMILES column
bladder.drugs['fingerprint']= bladder.drugs['canSMILES'].apply(smiles_to_fingerprint)

Merge Data

Group transcriptomics data into an array based on the improve_sample_id. This is in prepration for the model because the input values must be in an array. Then drop any duplicated transriptomics values with the same array.

group_df = bladder.transcriptomics.groupby('improve_sample_id')['transcriptomics'].apply(list).reset_index()
group_df['transcriptomics'] = group_df['transcriptomics'].apply(lambda x: list(set(x)))

group_df.head()

	improve_sample_id	transcriptomics
0	3290	[0.005498, 0.000412, 0.008812, 0.025414, 0.000...
1	3294	[0.035846, 0.000681, 0.006478, 0.001273, 0.004...
2	3299	[0.001088, 0.004633, 0.017393, 0.010687, 0.001...
3	3302	[0.001477, 0.003855, 0.000815, 0.019995, 0.008...
4	3305	[0.000516, 0.002067, 0.000333, 0.003617, 0.001...

We will merge the split dataframes with the bladder.drugs dataframe based on the improve_drug_id first and then in the next step with will merge the dataframes with the transcriptomics.

split.train.experiments = pd.merge(split.train.experiments[['improve_sample_id','improve_drug_id', 'dose_response_value']],
                                   bladder.drugs[['improve_drug_id','fingerprint']],
                                   on='improve_drug_id',
                                   how= 'left')

split.test.experiments = pd.merge(split.test.experiments[['improve_sample_id','improve_drug_id', 'dose_response_value']],
                                   bladder.drugs[['improve_drug_id','fingerprint']],
                                   on='improve_drug_id',
                                   how= 'left')

split.validate.experiments = pd.merge(split.validate.experiments[['improve_sample_id','improve_drug_id', 'dose_response_value']],
                                   bladder.drugs[['improve_drug_id','fingerprint']],
                                   on='improve_drug_id',
                                   how= 'left')

To the split bladder.experiments dataframes add transcriptomics arrays based on improve_sample_id. Then drop all columns besides our input values, transcriptomics and fingerprint, with the corresponding outputs dose_response_value. This is repeated for the train, test, and validate dataframes.

train_df = pd.merge(split.train.experiments, group_df, on='improve_sample_id', how = 'left')
train_df = train_df.dropna(subset=['transcriptomics'])
train_df = train_df.drop_duplicates(subset=['dose_response_value'])
train_df = train_df[['dose_response_value','transcriptomics', 'fingerprint' ]]

test_df = pd.merge(split.test.experiments, group_df, on='improve_sample_id', how = 'left')
test_df = test_df.dropna(subset=['transcriptomics'])
test_df = test_df.drop_duplicates(subset=['dose_response_value'])
test_df = test_df[['dose_response_value','transcriptomics','fingerprint']]

validate_df = pd.merge(split.validate.experiments, group_df, on='improve_sample_id', how = 'left')
validate_df = validate_df.dropna(subset=['transcriptomics'])
validate_df = validate_df.drop_duplicates(subset=['dose_response_value'])
validate_df = validate_df[['dose_response_value','transcriptomics','fingerprint']]

train_df.head()

	dose_response_value	transcriptomics	fingerprint
0	0.8081	[0.003186, 0.000027, 0.002066, 0.001147, 0.006...	[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...
50	0.9275	[0.02699, 0.003434, 0.001322, 0.001534, 0.0003...	[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, ...
489	0.3130	[0.003186, 0.000027, 0.002066, 0.001147, 0.006...	[0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...
872	0.9181	[0.044913, 0.000573, 0.002809, 0.000002, 0.005...	[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...
1036	0.4085	[0.0, 0.006175, 0.001342, 0.003682, 0.018603, ...	[0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...

test_df.head()

	dose_response_value	transcriptomics	fingerprint
164	0.8509	[0.000769, 0.002661, 0.002534, 0.000642, 0.006...	[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...
433	0.9100	[0.00862, 0.003253, 0.010242, 0.006115, 0.0015...	[0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...
1038	0.4285	[0.001088, 0.004633, 0.017393, 0.010687, 0.001...	[0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...
1240	0.4340	[0.003186, 0.000027, 0.002066, 0.001147, 0.006...	[0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...
2066	0.6813	[0.02699, 0.003434, 0.001322, 0.001534, 0.0003...	[0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, ...

validate_df.head()

	dose_response_value	transcriptomics	fingerprint
204	0.9430	[0.000768, 0.001443, 0.00352, 0.001669, 0.0018...	[0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, ...
293	0.8663	[0.003186, 0.000027, 0.002066, 0.001147, 0.006...	[0, 1, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...
588	0.3577	[0.00862, 0.003253, 0.010242, 0.006115, 0.0015...	[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...
677	0.9399	[0.002719, 0.002842, 0.002965, 0.002883, 0.002...	[0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, ...
766	0.7320	[0.000516, 0.002067, 0.000333, 0.003617, 0.001...	[0, 1, 0, 0, 1, 0, 1, 1, 0, 0, 0, 1, 0, 0, 0, ...

Quick Tip: View the length of the transcriptomics array to then later decide the shape that is input into the model. Here we have multiple size arrays so we will need to pad or truncate to ensure a consist shape. An example of that will be shown later in the tutorial.

pd.DataFrame({"Array Length":validate_df['transcriptomics'].apply(len).unique()})

	Array Length
0	3996
1	4741
2	4040
3	4766
4	3964
5	195
6	4207
7	4067
8	3723
9	4184
10	3698
11	4083
12	4236
13	2985
14	4193
15	4845
16	4491

Set the Model Parameters

Here we set the model parameters. You may modify the inputs, layers, activation functions, and outputs of the model. This is a general setup when using a keras model, but this may need to be adjusted to better suit your data. Here we use activation = 'relu' for the output layer, note that this activation is only intended when dealing with non-negative values (how the transcriptomic inputs were scaled previously).

#Define Inputs
transcriptomics_input = keras.Input(shape=(3000,), name="transcriptomics")
fingerprint_input = keras.Input(shape=(1024,), name="fingerprint")

# You may adjust these layers and activation functions to get better (or worse) results.
# Merge all available features into a single large vector via concatenation
x = layers.concatenate([transcriptomics_input, fingerprint_input])
x = layers.Dense(16, activation='relu')(x)
x = layers.Dense(64, activation="relu")(x)
x = layers.Dropout(0.5)(x)
x = layers.Dense(32, activation='relu')(x)
x = layers.Dense(12, activation='relu')(x)


# Priority prediction is here. You may add layers (and an output) after this too.
# This just allows for an early prediction in case you have a large datset and want preliminary results.
priority_pred = layers.Dense(1, name="priority", activation ='relu')(x)  # Using the activation relu.

# Instantiate an end-to-end model predicting both priority and department
model = keras.Model(
    inputs= [transcriptomics_input, fingerprint_input],
    outputs={"priority": priority_pred},
)

Compile the Model

Here we used a basic keras compile() input with the Adam optimizer. For other functions options visit keras.io/compile-method or optimizer options keras.io/optimizers. Then view the model summary, this can be helpful when adjusting parameters and seeing the model layers. For more information on the model summary visit keras.io/summary-method.

model.compile(
    optimizer=keras.optimizers.Adam(learning_rate=0.0001),
    loss={
        "priority": keras.losses.MeanSquaredError(),
    },
    metrics=[MeanAbsoluteError()],
    loss_weights={"priority": 1},
)

model.summary()

Model: "functional"

┏━━━━━━━━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━━━┳━━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━━━┓
┃ Layer (type)        ┃ Output Shape      ┃    Param # ┃ Connected to      ┃
┡━━━━━━━━━━━━━━━━━━━━━╇━━━━━━━━━━━━━━━━━━━╇━━━━━━━━━━━━╇━━━━━━━━━━━━━━━━━━━┩
│ transcriptomics     │ (None, 3000)      │          0 │ -                 │
│ (InputLayer)        │                   │            │                   │
├─────────────────────┼───────────────────┼────────────┼───────────────────┤
│ fingerprint         │ (None, 1024)      │          0 │ -                 │
│ (InputLayer)        │                   │            │                   │
├─────────────────────┼───────────────────┼────────────┼───────────────────┤
│ concatenate         │ (None, 4024)      │          0 │ transcriptomics[… │
│ (Concatenate)       │                   │            │ fingerprint[0][0] │
├─────────────────────┼───────────────────┼────────────┼───────────────────┤
│ dense (Dense)       │ (None, 16)        │     64,400 │ concatenate[0][0] │
├─────────────────────┼───────────────────┼────────────┼───────────────────┤
│ dense_1 (Dense)     │ (None, 64)        │      1,088 │ dense[0][0]       │
├─────────────────────┼───────────────────┼────────────┼───────────────────┤
│ dropout (Dropout)   │ (None, 64)        │          0 │ dense_1[0][0]     │
├─────────────────────┼───────────────────┼────────────┼───────────────────┤
│ dense_2 (Dense)     │ (None, 32)        │      2,080 │ dropout[0][0]     │
├─────────────────────┼───────────────────┼────────────┼───────────────────┤
│ dense_3 (Dense)     │ (None, 12)        │        396 │ dense_2[0][0]     │
├─────────────────────┼───────────────────┼────────────┼───────────────────┤
│ priority (Dense)    │ (None, 1)         │         13 │ dense_3[0][0]     │
└─────────────────────┴───────────────────┴────────────┴───────────────────┘

 Total params: 67,977 (265.54 KB)

 Trainable params: 67,977 (265.54 KB)

 Non-trainable params: 0 (0.00 B)

Run your Model

Before you can run the mdoel you may need to pad/truncate the input arrays. As mentioned earlier this is because the arrays must be the same length in order for the model to run correctly. Here we adjust both the transcriptomics and fingerprint inputs for all split dataframes.

# Set target_length as the maximum length
def pad_sequence(seq, target_length, pad_value=0):
    """Pad or truncate a sequence to a target length."""
    if len(seq) > target_length:  # Truncate if longer than target_length
        return seq[:target_length]
    elif len(seq) < target_length:  # Pad with zeros if shorter than target_length
        return seq + [pad_value] * (target_length - len(seq))
    else:
        return seq

# Example: Create fixed-length arrays
train_df['transcriptomics'] = train_df['transcriptomics'].apply(lambda x: pad_sequence(x, target_length=3000))
test_df['transcriptomics'] = test_df['transcriptomics'].apply(lambda x: pad_sequence(x, target_length=3000))
validate_df['transcriptomics'] = validate_df['transcriptomics'].apply(lambda x: pad_sequence(x, target_length=3000))

train_df['fingerprint'] = train_df['fingerprint'].apply(lambda x: pad_sequence(x, target_length=1024))
test_df['fingerprint'] = test_df['fingerprint'].apply(lambda x: pad_sequence(x, target_length=1024))
validate_df['fingerprint'] = validate_df['fingerprint'].apply(lambda x: pad_sequence(x, target_length=1024))

Define the output data, dose_response_value, in prepartion of running the model.

y_train = np.array(train_df['dose_response_value'])
y_test = np.array(test_df['dose_response_value'])
y_validate = np.array(validate_df['dose_response_value'])

Now use the keras fit() function with the train and validation data. For more information on this function visit kera.io/fit-method. Note you may adjust the epoch and bath_size to better fit your model.

# Train the model
history = model.fit(
    [np.array(train_df['transcriptomics'].tolist()),
     np.array(train_df['fingerprint'].tolist())
     ],
    y_train,
    validation_data=([np.array(validate_df['transcriptomics'].tolist()),
                      np.array(validate_df['fingerprint'].tolist())
                      ],
                     y_validate),
    epochs=100,  # Number of epochs
    batch_size=32,  # Batch size
    verbose=0 # Reduces the output (set equal to 1 if you want to see the output steps)
)

Evaluate your Model

Now we will use the test data to evaluate our model. Information on the evaluate() function, kera.io/evaluate-method.

losses = model.evaluate(
    [np.array(test_df['transcriptomics'].tolist()),
     np.array(test_df['fingerprint'].tolist())],
    y_test,
    return_dict=True)
print(losses)

3/3 ━━━━━━━━━━━━━━━━━━━━ 0s 4ms/step - loss: 0.0160 - mean_absolute_error: 0.1053
{'loss': 0.016004862263798714, 'mean_absolute_error': 0.10526344180107117}

Here is a site giving a guide on how to evaluate your model, visit machinelearning.com. The site gives examples of similary graphs shown below demonstrating what overfitting, underfitting, and correct models may output.

Plot Loss Function and Error of the Model

# summarize history for accuracy
plt.plot(history.history['val_mean_absolute_error'])
plt.plot(history.history['mean_absolute_error'])
plt.title('Mean Absolute Error of Deep Learning Model')
plt.ylabel('Mean Absolute Error')
plt.xlabel('epoch')
plt.legend(['Validation Set','Training Set',], loc='upper right')
plt.show()
# summarize history for loss
plt.plot(history.history['val_loss'])
plt.plot(history.history['loss'])
plt.title('Loss Function of Deep Learning Model')
plt.ylabel('Loss')
plt.xlabel('epoch')
plt.legend(['Validation Set','Training Set'], loc='upper right')
plt.show()

Prediction vs True Values

We will now use the test data to see how our model did in making predictions based off of the train and validation data it was given. Using graph and summary statistics we can see how well the model was able to reflect the given data.

# make a prediction
predict = model.predict([np.array(test_df['transcriptomics'].tolist()),
     np.array(test_df['fingerprint'].tolist())])

3/3 ━━━━━━━━━━━━━━━━━━━━ 0s 15ms/step

new_df = pd.DataFrame({
    'Predicted Values': predict["priority"].tolist(),
    'True Values': y_test
})
new_df['Predicted Values'] = new_df['Predicted Values'].apply(lambda x: x[0] if isinstance(x, list) and len(x) == 1 else x)
sorted_df = new_df.sort_values(by='True Values', ascending=True)  # Change 'ascending' to False for descending order
sorted_df.head()

	Predicted Values	True Values
18	0.548078	0.2845
83	0.668909	0.3339
41	0.378916	0.3577
62	0.373495	0.3597
19	0.465165	0.3733

true_values = np.array(sorted_df["True Values"])
predicted_values = np.array(sorted_df["Predicted Values"])

# Plot the true values and predicted values
plt.plot(true_values)
plt.plot(predicted_values)

# Fit trendlines
true_values_trendline = LinearRegression().fit(np.arange(len(true_values)).reshape(-1, 1), true_values.reshape(-1, 1)).predict(np.arange(len(true_values)).reshape(-1, 1))
predicted_values_trendline = LinearRegression().fit(np.arange(len(predicted_values)).reshape(-1, 1), predicted_values.reshape(-1, 1)).predict(np.arange(len(predicted_values)).reshape(-1, 1))

# Plot trendlines
plt.plot(true_values_trendline, linestyle='--', color='blue', alpha=0.3)
plt.plot(predicted_values_trendline, linestyle='--', color='orange', alpha=0.3)

# Customize plot
plt.title('Predicted AUC vs True AUC')
plt.ylabel('Drug Response AUC')
plt.xlabel('Observation Values')
plt.legend(['True Values', 'Predicted Values', 'True Values Trendline', 'Predicted Values Trendline'], loc='lower right')
plt.show()

true_values = np.array(new_df["True Values"])
predicted_values = np.array(new_df["Predicted Values"])

# Plot the true values and predicted values
plt.plot(true_values)
plt.plot(predicted_values)

# Fit trendlines
true_values_trendline = LinearRegression().fit(np.arange(len(true_values)).reshape(-1, 1), true_values.reshape(-1, 1)).predict(np.arange(len(true_values)).reshape(-1, 1))
predicted_values_trendline = LinearRegression().fit(np.arange(len(predicted_values)).reshape(-1, 1), predicted_values.reshape(-1, 1)).predict(np.arange(len(predicted_values)).reshape(-1, 1))


# Customize plot
plt.title('Predicted AUC vs True AUC')
plt.ylabel('Drug Response AUC')
plt.xlabel('Observation Values')
plt.legend(['True Values', 'Predicted Values', 'True Values Trendline', 'Predicted Values Trendline'], loc='lower right')
plt.show()

import seaborn as sns

# Visualizing the regression results
plt.figure(figsize=(10,6))
sns.scatterplot(x=y_test, y=predicted_values, alpha=0.6)
sns.lineplot(x=[y_test.min(), y_test.max()], y=[y_test.min(), y_test.max()], color='orange')
plt.title('True vs Predicted AUC Values')
plt.xlabel('True AUC Values')
plt.ylabel('Predicted AUC Values')
plt.show()

# Calculate Summary Statistics


# Calculate Mean Absolute Error (MAE)
mae = mean_absolute_error(new_df['True Values'], new_df['Predicted Values'])

# Calculate Root Mean Squared Error (RMSE)
rmse = np.sqrt(mean_squared_error(new_df['True Values'], new_df['Predicted Values']))

# Calculate R-squared (R2) score
r2 = r2_score(new_df['True Values'], new_df['Predicted Values'])

summary_statistics = new_df.describe()

# Print the statistics
print("Mean Absolute Error (MAE):", mae)
print("Root Mean Squared Error (RMSE):", rmse)
print("R-squared (R2) Score:", r2)
print("\n")
# Print summary statistics
print(summary_statistics)

Mean Absolute Error (MAE): 0.10526344083256191
Root Mean Squared Error (RMSE): 0.1265103241081565
R-squared (R2) Score: 0.5736720035806125


       Predicted Values  True Values
count         90.000000    90.000000
mean           0.661669     0.717819
std            0.139365     0.194841
min            0.349520     0.284500
25%            0.576000     0.580650
50%            0.671961     0.751550
75%            0.764141     0.896650
max            0.968248     0.998800

#Side by side comparison for first 50 values.
new_df[0:50].head()

	Predicted Values	True Values
0	0.691072	0.8509
1	0.968248	0.9100
2	0.542540	0.4285
3	0.449587	0.4340
4	0.605799	0.6813