Training data was obtained from the 2019 DREAM Malaria Challenge ^{8, 9} . The training data consists of gene expression data of 5,540 genes of 30 isolates from the malaria parasite, Plasmodium falciparum. For each malaria parasite isolate, transcription data was collected at two time points [6 hours post invasion (hpi) and 24 hpi], with and without treatment of dihydroartemisinin (the metabolically active form of artemisinin), each with a biological replicate. This yields a total of at eight data points for each isolate. The initial form of the training dataset contains 272 rows and 5,546 columns, as shown in Table 1.

The transcription data was collected as described in Table 2. The transcription data set consists of 92 non-coding RNAs (denoted by gene IDs that begins with ’MAL’), while the rest are protein coding genes (denoted by gene IDs that start with ’PF3D7’). The feature to predict is DHA_ IC50.

We used Apache Spark ¹⁰ to pivot the dataset such that each isolate was its own row and each of the transcription values for each gene and attributes (i.e. timepoint, treatment, biological replicate) combination was its own column. This exercise transformed the training dataset from 272 rows and 5,546 columns to 30 rows and 44,343 columns, as shown in Table 3. We completed this pivot by slicing the data by each of the eight combinations of timepoint, treatment, and biological replicate, dynamically renaming the variables (genes) for each slice, and then joining all eight slices back together.

Example code shown below in the section labeled code 1. By using the massively parallel architecture of Spark, this transformation can be completed in a minimal amount of time on a relatively small cluster environment (e.g., <10 minutes using a 8-worker/36-core cluster with PySpark on Apache Spark 2.4.3).

1 ## Separate Dependent Variable 2 y = train . select(col( "Isolate" ), 3 col( "DHA_IC50" )) \ 4 . distinct() 5 6 ## Create Slice [Timepoint: 24HR, Treatment: DHA, BioRep: BRep1] 7 hr24_trDHA_br1 = train . drop( "Sample_Name" , "DHA_IC50" ) \ 8 . filter((col( "Timepoint" ) == "24HR" ) & 9 (col( "Treatment" ) == "DHA" ) & 10 (col( "BioRep" ) == "BRep1" )) 11 ## Rename Columns 12 column_list = hr24_trDHA_br1 . columns 13 prefix = "hr24_trDHA_br1_" 14 new_column_list = [prefix + s if s != "Isolate" else s for s in column_list] 15 16 column_mapping = [[o, n] for o, n in zip (column_list, new_column_list)] 17 18 hr24_trDHA_br1 = hr24_trDHA_br1 . select( list ( map ( lambda old, new: col(old) \ 19 . alias(new), * zip ( * column_mapping)))) 20 21 ## Join Slices Together 22 data = y . join(hr24_trDHA_br1, "Isolate" , how = 'left' ) \ 23 . join(hr24_trDHA_br2, "Isolate" , how = 'left' ) \ 24 . join(hr24_trUT_br1, "Isolate" , how = 'left' ) \ 25 . join(hr24_trUT_br2, "Isolate" , how = 'left' ) \ 26 . join(hr6_trDHA_br1, "Isolate" , how = 'left' ) \ 27 . join(hr6_trDHA_br2, "Isolate" , how = 'left' ) \ 28 . join(hr6_trUT_br1, "Isolate" , how = 'left' ) \ 29 . join(hr6_trUT_br2, "Isolate" , how = 'left' )

Lastly, the dataset is then vectorized using the Spark VectorAssembler, and converted into a Numpy ¹² -compatible array. Example code shown below in Code 1. Vectorization allows for highly scalable parallelization of the machine learning modeling in the next step.

1 ## Transform Data using VectorAssembler 2 assemblerInputs = numericalColumns 3 assembler = VectorAssembler(inputCols = assemblerInputs, outputCol = "features" ) \ 4 . setHandleInvalid( "keep" ) 5 stages += [assembler] 6 7 prepPipeline = Pipeline() . setStages(stages) 8 pipelineModel = prepPipeline . fit(data) 9 vectordata = pipelineModel . transform(data) \ 10 . select(col( "DHA_IC50" ), col( "features" )) \ 11 . withColumnRenamed( "DHA_IC50" , "label" ) 12 13 ## Convert to Numpy Array 14 import numpy as np 15 pddata = vectordata . toPandas() 16 seriesdata = pddata[ 'features' ] . apply( lambda x : np . array(x . toArray())) \ 17 . as_matrix() . reshape( -1 , 1 ) 18 X_train = np . apply_along_axis( lambda x : x[ 0 ], 1 , seriesdata) 19 y_train = pddata[ 'label' ] . values . reshape( -1 , 1 ) . ravel() 20 21 ## Example Output (X_train) After Vectorization 22 array([[ -0.62161893 , -0.60860881 , -1.11331369 , ... , -1.457377 , 23 -3.292903 , -1.869169 ], 24 [ -0.55719008 , -2.41660489 , -1.39244109 , ... , -1.770098 , 25 -3.698841 , -1.740082 ], 26 ... , 27 [ -0.17072536 , -2.32828532 , -1.08406554 , ... , -1.402658 , 28 -5.314896 , -1.328886 ], 29 [ -0.1923732 , -1.88763881 , -1.23867258 , ... , -1.971246 , 30 -3.567355 , -1.904116 ]])

We used the Microsoft Azure Machine Learning Service ¹³ as the tracking platform for retaining model performance metrics as the various models were generated. For this use case, 498 machine learning models were trained using various scaling techniques and algorithms. We then created two ensemble models of the individual models using Stack Ensemble and Voting ensemble methods. Scaling and normalization methods are shown in Table 14.

The Microsoft AutoML package ¹⁴ allows for the parallel creation and testing of various models, fitting based on a primary metric. For this use case, models were trained using Decision Tree, Elastic Net, Extreme Random Tree, Gradient Boosting, Lasso Lars, LightGBM, RandomForest, and Stochastic Gradient Decent algorithms along with various scaling methods from Maximum Absolute Scaler, Min/Max Scaler, Principal Component Analysis, Robust Scaler, Sparse Normalizer, Standard Scale Wrapper, Truncated Singular Value Decomposition Wrapper (as defined in Table 14). All of the machine learning algorithms are from the scikit-learn package ¹⁵ except for LightGBM, which is from the LightGBM package ¹⁶ . The settings for the model sweep are defined in Table 4. The ‘Preprocess Data?’ parameter enables the scaling and imputation of the features in the data.

Once the 498 individual models were trained, two ensemble models (voting ensemble and stack ensemble) were then created and tested. The voting ensemble method makes a prediction based on the weighted average of the previous models’ predicted regression outputs whereas the stacking ensemble method combines the previous models and trains a meta-model using the elastic net algorithm based on the output from the previous models. The model selection method used was the Caruana ensemble selection algorithm ¹⁷ .

The voting ensemble model (using soft voting) was selected as the best model, having the lowest normalized Root Mean Squared Error (RMSE), as shown in Table 5. The top 10 models trained are reported in Table 6. Having a normalized RMSE of only 0.1228 and a Mean Absolute Percentage Error (MAPE) of 24.27%, this model is expected to accurately predict IC ₅₀ in malaria isolates. See Figure 1 for a visualization of the experiment runs and Figure 2 for the distribution of residuals on the best model.

An in vivo transcription data set from Mok et al., (2015) Science ¹⁸ was used to predict the parasite clearance rate of malaria parasite isolates based on in vitro transcriptional profiles (see Table 8).

The training data consists of 1,043 isolates with 4,952 genes from the malaria parasite Plasmodium falciparum. For each malaria parasite isolate, transcription data was collected for various PF3D7 genes. The form of the training dataset contains 1,043 rows and 4,957 columns, as shown in Table 7. The feature to predict is ClearanceRate.

The training data for this use case did not require the same pivoting transformations as in the last use case as each record describes a single isolate. Thus, only the vectorization of the data was necessary, which was performed using the Spark VectorAssembler and then converted into a Numpy-compatible array ¹² . Example code is shown below. Note that this vectorization only kept the numerical columns, which excludes the Country, Kmeans_Grp, and Asexual_stage__hpi_ attributes as they are either absent or contain non-matching factors (i.e. different set of countries) in the testing data.

1 ## Transform Data using VectorAssembler 2 assemblerInputs = numericalColumns 3 assembler = VectorAssembler(inputCols = assemblerInputs, outputCol = "features" ) \ 4 . setHandleInvalid( "keep" ) 5 stages += [assembler] 6 7 prepPipeline = Pipeline() . setStages(stages) 8 pipelineModel = prepPipeline . fit(data) 9 vectordata = pipelineModel . transform(data) 10 11 ## Convert to Numpy Array 12 import numpy as np 13 pddata = vectordata . toPandas() 14 seriesdata = pddata[ 'features' ] . apply( lambda x : np . array(x . toArray())) \ 15 . as_matrix() . reshape( -1 , 1 ) 16 X_train = np . apply_along_axis( lambda x : x[ 0 ], 1 , seriesdata) 17 y_train = pddata[ 'label' ] . values . reshape( -1 , 1 ) . ravel() 18 19 ## Example Output (X_train) After Vectorization 20 array([[ 0.2263112 , -0.39682897 , -1.80458125 , ... , nan, 21 -1.30952803 , -0.64170958 ] , 22 [ 0.55442743 , 0.54200115 , -1.56157279 , ... , 1.83083869 , 23 0.21021662 , -1.06553341 ], 24 ... , 25 [ 1.24446867 , -0.09076431 , -1.62156926 , ... , 3.18060844 , 26 -0.43056353 , nan], 27 [ 0.99909549 , -1.47208829 , -1.91898139 , ... , 2.59463935 , 28 -1.21233458 , nan]])

Once the 98 individual models were trained, two ensemble models (voting ensemble and stack ensemble) were then created and tested as before. Model search parameters are shown in Table 9.

The voting ensemble model (using soft voting) was selected as the best model, having the highest area under the receiver operating characteristic curve (AUC), as shown in Table 11. The top 10 of the 100 models trained are reported in Table 10. Having a weighted AUC of 0.87 and a weighted F1 score of 0.80, this model is expected to accurately predict isolate clearance rates. A confusion matrix of the predicted results versus actuals is shown in Table 12. See Figure 3 for a visualization of the experiment runs and see Figure 4 and Figure 5 for the ROC and Precision-Recall curves on the best model.

subsectionFeature importance Feature importances were calculated using mimic-based model explanation of the ensemble model ¹⁹ . The mimic explainer works by training global surrogate models to mimic blackbox models. The surrogate model is an interpretable model, trained to approximate the predictions of a black box model as accurately as possible. See Figure 6 and Table 13.

The challenge datasets are available from Synapse ( https://www.synapse.org/; Synapse ID: syn18089524). Access to the data requires registration and agreement to the conditions for use at: https://www.synapse.org/#!Synapse: syn18089524.

Challenge documentation, including the detailed description of the Challenge design, data description, and overall results can be found at: https://www.synapse.org/#!Synapse:syn16924919/wiki/583955.

Whole genome expression profiling of artemsinin-resistant Plasmodium falciparum field isolates, Accession number GSE59099: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE59099.

Zenodo: colbyford/malaria_DREAM2019: Ensemble Machine Learning Modeling for the Prediction of Artemisinin Resistance in Malaria - Initial Code Release for Research Publication (F1000). https://doi.org/10.5281/zenodo.3590459 ²¹ .

This project contains the following underlying data:

Data are available under the terms of the Creative Commons Zero "No rights reserved" data waiver (CC0 1.0 Public domain dedication).