MSBA Capstone

Introduction

Overview

In the biopharmaceutical industry, bringing a new drug to market represents a significant financial investment and commercial risk. The primary driver of rising R&D costs is the high attrition rate of drug candidates during clinical trials. While a compound may demonstrate efficacy in controlled laboratory environments, human genetic variation often introduces unforeseen complexities. Pathogenic mutations can alter a target protein’s structural integrity, rendering a drug ineffective or toxic across patient populations.

Currently, pharmaceutical companies invest substantial resources and time physically validating biological targets that may not work. This project addresses that bottleneck using predictive analytics. By calculating a clinical success probability score for human proteins before entering lab testing phase, biopharmaceutical stakeholders can avoid high-risk targets and strategically allocate R&D funds towards candidates with the highest statistical probability of structural success.

This project will integrate publicly available datasets using R programming language. The first dataset is AlphaMissense (via Bioconductor), developed by Google DeepMind, which contains predictions for over 71 million possible single amino acid substitutions across the human genome. Each variant is assigned a probability score between 0 (benign) and 1 (pathogenic) representing the likelihood that the mutation will disrupt protein function. The dataset is available through Bioconductor: https://bioconductor.org/packages/release/bioc/html/AlphaMissenseR.html

Another dataset that will be included is the Therapeutic Target Database (TTD), an open-access resource that tracks therapeutic targets, associated drugs, and clinical development outcomes. This data is available through the Therapeutic Target Database: https://ttd.idrblab.cn/full-data-download

UniProt data, central database for curated protein information that includes manually reviewed protein sequences and functional information, from https://www.uniprot.org/ will be used to connect and join the former datasets.

The technical workflow will be executed in R Studio, leveraging advanced data wrangling and modeling techniques. By connecting to the local AlphaMissense database, downloaded locally via the AlphaMissenseR package, the project will extract protein level features such as pathogenicity scores, and the proportion of highly volatile mutation regions. These risk indicators will then be merged with historical drug development outcomes from the Therapeutic Target Database using shared protein identifiers.

Following data integration and feature engineering, machine learning techniques such as random forest and neural network models will be applied to predict the likelihood of clinical drug target success or failure. The resulting predictive framework may provide valuable insights into how genetic variability influences drug development outcomes and serve as a decision tool for early-stage target prioritization.

This project will utilize skills acquired throughout the Business Analytics curriculum. Data integration, joins, and data manipulation techniques learned in BANA 7025 (Data Wrangling) will be used to prepare, clean, and merge the datasets. Predictive modeling methods, including random forests and classification techniques, will integrate concepts covered in BANA 7046 and BANA 7047 (Data Mining). Additionally, regression analysis and data visualization techniques acquired across multiple BANA courses will be applied to evaluate model performance and communicate results. The expected outcome of this project is a predictive model that quantifies the development risk associated with therapeutic protein targets based on genetic variation data and estimates their likelihood of clinical trial success. Such insights may help pharmaceutical organizations make more informed investment decisions, reduce costly late-stage failures, and improve the efficiency of the drug development pipeline.

Data Prep

Packages used:

library(tidyverse)
library(stringr)
library(AlphaMissenseR)
library(readxl)
library(knitr)
library(kableExtra)
library(randomForest)
library(GGally)
library(gbm)
library(pROC)
library(ROCR)
library(mgcv)
library(nnet)
library(neuralnet)
library(caret)
setwd("C:/Users/tanje/OneDrive/Desktop/uc grad/capstone 2026/")


Available datasets from Google DeepMind Alphamissense Bioconductor hosted on Zenodo.

am_available() |>
  kable() |>
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE) |>
  scroll_box(width = "100%")
record key size cached filename link
10813168 gene_hg38 253636 FALSE AlphaMissense_gene_hg38.tsv.gz https://zenodo.org/api/records/10813168/files/AlphaMissense_gene_hg38.tsv.gz/content
10813168 isoforms_hg38 1177361934 FALSE AlphaMissense_isoforms_hg38.tsv.gz https://zenodo.org/api/records/10813168/files/AlphaMissense_isoforms_hg38.tsv.gz/content
10813168 isoforms_aa_substitutions 2461351945 FALSE AlphaMissense_isoforms_aa_substitutions.tsv.gz https://zenodo.org/api/records/10813168/files/AlphaMissense_isoforms_aa_substitutions.tsv.gz/content
10813168 hg38 642961469 FALSE AlphaMissense_hg38.tsv.gz https://zenodo.org/api/records/10813168/files/AlphaMissense_hg38.tsv.gz/content
10813168 hg19 622293310 FALSE AlphaMissense_hg19.tsv.gz https://zenodo.org/api/records/10813168/files/AlphaMissense_hg19.tsv.gz/content
10813168 gene_hg19 243943 TRUE AlphaMissense_gene_hg19.tsv.gz https://zenodo.org/api/records/10813168/files/AlphaMissense_gene_hg19.tsv.gz/content
10813168 aa_substitutions 1207278510 TRUE AlphaMissense_aa_substitutions.tsv.gz https://zenodo.org/api/records/10813168/files/AlphaMissense_aa_substitutions.tsv.gz/content


Proceed using “aa_substitutions” file holding uniprot ID and pathogenicity. Note: Multi protein isoform not used.

raw_amdata <- am_data("aa_substitutions") #remote database connection with protein canonical versions 
kable(head(raw_amdata, 2))
uniprot_id protein_variant am_pathogenicity am_class
A0A024R1R8 M1A 0.4673 ambiguous
A0A024R1R8 M1C 0.3828 ambiguous


Data cleaning.

protein_risk_features <- raw_amdata |> 
  group_by(uniprot_id) |> 
  summarize(mean_pathogenicity = mean(am_pathogenicity, na.rm = TRUE), #average mutations
    pct_pathogenic_mutations = mean(as.integer(am_class == "pathogenic"), na.rm = TRUE),
    pathogenicity_sd = sd(am_pathogenicity, na.rm = TRUE), #risk variability across protein's structure.
    .groups = "drop") |>
  rename(accession_id = uniprot_id)
kable(head(protein_risk_features,2))
accession_id mean_pathogenicity pct_pathogenic_mutations pathogenicity_sd
A0A0G2JN01 0.3639831 0.2259307 0.2563501
A0A0G2JNJ9 0.5270716 0.4429825 0.2549995

Higher pathogenicity values show high likelihood of genetic variation.


Connect accession_id to uniprot_id with data from uniprot website. Cleaning uniprot data:

uniprot <- read_delim("ttd data/uniprotkb_organism_id_9606_AND_reviewed_2026_06_29.tsv", delim = "\t", show_col_types = FALSE)
kable(head(uniprot,2))
Entry Entry Name Protein names Gene Names Organism
A0A087X1C5 CP2D7_HUMAN Cytochrome P450 2D7 (EC 1.14.14.1) CYP2D7 Homo sapiens (Human)
A0A096LP01 SIM26_HUMAN Small integral membrane protein 26 SMIM26 LINC00493 Homo sapiens (Human) 
uniprot_cleaned <- uniprot |> 
  select(uniprot_id = `Entry Name`, accession_id = Entry)
kable(head(uniprot_cleaned,2))
uniprot_id accession_id
CP2D7_HUMAN A0A087X1C5
SIM26_HUMAN A0A096LP01


Join AlphaMissense table with accessionid / uniprotid. This will aid in connecting data with clinical trial data files.

protein_risk_features_local <- protein_risk_features |> 
  collect() 

pat_missense <- protein_risk_features_local |>
  left_join(uniprot_cleaned, by = "accession_id")

kable(head(pat_missense,2))
accession_id mean_pathogenicity pct_pathogenic_mutations pathogenicity_sd uniprot_id
A0A0C4DH43 0.4155238 0.2972136 0.2660014 HV70D_HUMAN
A0A0G2JLQ3 0.3591432 0.2096607 0.2457163 NA

Joining With Clinical Trial Data

The TTD_target_download file provides data on trial outcomes. It is joined with the Missense file for analysis.


Data preparation.

lines <- readLines("ttd data/P1-01-TTD_target_download.txt") # includes header
data_lines <- lines[32:length(lines)] 
data_lines <- data_lines[str_detect(data_lines, "\t")] # keep lines with tabs
head(data_lines)
## [1] "T47101\tTARGETID\tT47101"                                     
## [2] "T47101\tFORMERID\tTTDC00024"                                  
## [3] "T47101\tUNIPROID\tFGFR1_HUMAN"                                
## [4] "T47101\tTARGNAME\tFibroblast growth factor receptor 1 (FGFR1)"
## [5] "T47101\tGENENAME\tFGFR1"                                      
## [6] "T47101\tTARGTYPE\tSuccessful"
ttd_table <- tibble(line = data_lines) |>
  filter(str_detect(line, "TARGETID|FORMERID|UNIPROID|TARGNAME|TARGTYPE")) |>
  separate(line, into = c("TargetID", "Field", "Value"), sep = "\t") |> # split line into 3 col
  pivot_wider(names_from = Field, values_from = Value) # make fields col

kable(head(ttd_table,2))
TargetID TARGETID FORMERID UNIPROID TARGNAME TARGTYPE
T47101 T47101 TTDC00024 FGFR1_HUMAN Fibroblast growth factor receptor 1 (FGFR1) Successful
T59328 T59328 TTDS00355 EGFR_HUMAN Epidermal growth factor receptor (EGFR) Successful 
table(ttd_table$TARGTYPE)
## 
##      Clinical trial        Discontinued Literature-reported   Patented-recorded 
##                1327                  83                1845                 184 
##         Preclinical          Successful 
##                  57                 802


Recategorize target file into 4 categories: Preclinical, Terminated/Discontinued/Withdrawn, In Clinical Trial (Currently), Approved

ttd_cleaned <- ttd_table |>
  mutate(clinical_status = case_when(
    TARGTYPE == "Successful" ~ "Approved",
    TARGTYPE == "Clinical trial" ~ "In Clinical Trial (Currently)",
    TARGTYPE == "Discontinued" ~ "Terminated/Discontinued/Withdrawn",
    TARGTYPE %in% c("Preclinical", "Patented-recorded", "Literature-reported") ~ "Preclinical",
    TRUE ~ "Preclinical"),
    clinical_status = factor(clinical_status, levels = c("Preclinical", "Terminated/Discontinued/Withdrawn", "In Clinical Trial (Currently)", "Approved"), ordered = TRUE))
table(ttd_cleaned$clinical_status)
## 
##                       Preclinical Terminated/Discontinued/Withdrawn 
##                              2086                                83 
##     In Clinical Trial (Currently)                          Approved 
##                              1327                               802
kable(head(ttd_cleaned,2))
TargetID TARGETID FORMERID UNIPROID TARGNAME TARGTYPE clinical_status
T47101 T47101 TTDC00024 FGFR1_HUMAN Fibroblast growth factor receptor 1 (FGFR1) Successful Approved
T59328 T59328 TTDS00355 EGFR_HUMAN Epidermal growth factor receptor (EGFR) Successful Approved


Joining missense and trial data.

# unique IDs & remove NA 
ttd_unique <- ttd_cleaned |>
  filter(!is.na(UNIPROID)) |>
  distinct(UNIPROID, .keep_all = TRUE)

missense_unique <- pat_missense |>
  filter(!is.na(uniprot_id)) |>
  distinct(uniprot_id, .keep_all = TRUE)

final_ml_dataset <- missense_unique |>
  inner_join(ttd_unique, by = c("uniprot_id" = "UNIPROID"))

final_ml_dataset <- select(final_ml_dataset,-TARGETID) # remove duplicate col

final_ml_dataset |>
  head(5) |>
  kable() |>
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE) |>
  scroll_box(width = "100%", height = "300px")
accession_id mean_pathogenicity pct_pathogenic_mutations pathogenicity_sd uniprot_id TargetID FORMERID TARGNAME TARGTYPE clinical_status
O15164 0.6175176 0.5638596 0.3509092 TIF1A_HUMAN T56108 TTDI03589 Tripartite motif-containing 24 (TRIM24) Literature-reported Preclinical
O43526 0.6915094 0.6587086 0.3346199 KCNQ2_HUMAN T74483 TTDR00544 Voltage-gated potassium channel Kv7.2 (KCNQ2) Literature-reported Preclinical
O60706 0.6866594 0.6591349 0.3309775 ABCC9_HUMAN T02777 TTDS00332 ATP-binding cassette transporter C9 (ABCC9) Successful Approved
O75197 0.6085693 0.5586443 0.3408779 LRP5_HUMAN T86072 NA Low-density lipoprotein receptor-related protein 5 (LRP5) Clinical trial In Clinical Trial (Currently)
P07766 0.4918775 0.4080854 0.3152401 CD3E_HUMAN T87075 TTDS00473 T-cell surface glycoprotein CD3 epsilon (CD3E) Successful Approved

This table will be used for machine learning models and predictions.

Exploratory

Dataset Summary

final_ml_dataset  |> 
  collect() |>
  select(mean_pathogenicity,pct_pathogenic_mutations,pathogenicity_sd, clinical_status) |> 
  summary() |>
  kable()
mean_pathogenicity pct_pathogenic_mutations pathogenicity_sd clinical_status
Min. :0.2787 Min. :0.06599 Min. :0.1562 Preclinical :1348
1st Qu.:0.4663 1st Qu.:0.37840 1st Qu.:0.3126 Terminated/Discontinued/Withdrawn: 36
Median :0.5432 Median :0.47697 Median :0.3326 In Clinical Trial (Currently) : 892
Mean :0.5465 Mean :0.47771 Mean :0.3272 Approved : 530
3rd Qu.:0.6253 3rd Qu.:0.57905 3rd Qu.:0.3474 NA
Max. :0.9086 Max. :0.94645 Max. :0.3878 NA

Data Shape

dim(final_ml_dataset)
## [1] 2806   10

There are a total of 2806 observations across 10 variables. Mean pathogenicity shows symmetric distribution ranging from .28 to .91.


Distribution

 final_ml_dataset |> 
  ggplot(aes(x = mean_pathogenicity)) +
  geom_density(fill = "#E6E6FA", alpha = 0.7, color = NA) +
  theme_minimal() +
  labs(
    title = "Distribution of Mean Protein Pathogenicity",
    x = "Mean Pathogenicity",
    y = "Density"
  )

The plot shows a single peak distribution. Majority of the proteins have an average pathogenicity score concentrated between 0.40 and 0.65. Most proteins sit in the middle where specific mutations can break them, but many other mutations are tolerated. The left tail holds resilient mutation tolerant proteins, while right tail proteins are fragile. If TTD data shows that drugs targeting these specific right-tail proteins have a massive failure rate then the hypothesis stands.


Boxplot to compare pathogenicity distributions among trial status.

ggplot(final_ml_dataset, aes(x = clinical_status, y = mean_pathogenicity, fill = clinical_status)) +
  geom_boxplot(alpha = .8, size=.6) +
  scale_fill_brewer(palette = "Pastel1" , name = "Clinical Status") +
  theme_minimal() +
  labs(title = "Mean Pathogenicity by Clinical Status",
       x = "Clinical Status",
       y = "Mean Pathogenicity Score") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

Preclinical, in trial currently and approved show similar median for mean pathogenicity, but terminated group looks to be higher. This could indicate that proteins associated with clinical trial failures have higher mean pathogenicity scores than those that get approved.

Multicollinearity

final_ml_dataset |>
  select(mean_pathogenicity, pct_pathogenic_mutations, pathogenicity_sd) |>
  ggpairs(
    lower = list(continuous = wrap("points", color = "#A7C7E7", alpha = 0.5)),
    diag  = list(continuous = wrap("densityDiag", fill = "#A7D8F0", alpha = 0.4))
  )

Mean pathogenicity score and percentage of pathogenic mutations show multicollinearity. Only one variable will be used for a stable model.

Regression Model

Logistic Regression, taking approval / termination as binary outcomes.

historical_data <- final_ml_dataset |>
  filter(clinical_status %in% c("Approved", "Terminated/Discontinued/Withdrawn")) |>
  mutate(outcome = if_else(clinical_status == "Approved", 1, 0)) |> 
  mutate(outcome_factor = as.factor(outcome)) |> 
  arrange(accession_id)

dim(historical_data)
## [1] 566  12


Split training and testing data. 80% of the data is used to train the model, and the remaining 20% for testing the model.

set.seed(123)
index <- sample(nrow(historical_data), 0.8 * nrow(historical_data))
train <- historical_data[index, ]
test <- historical_data[-index, ]

logit_model <- glm(outcome ~ mean_pathogenicity + pathogenicity_sd, family = binomial, data = train)
summary(logit_model)
## 
## Call:
## glm(formula = outcome ~ mean_pathogenicity + pathogenicity_sd, 
##     family = binomial, data = train)
## 
## Coefficients:
##                    Estimate Std. Error z value Pr(>|z|)  
## (Intercept)           1.354      2.168   0.625   0.5323  
## mean_pathogenicity   -2.839      1.572  -1.806   0.0709 .
## pathogenicity_sd      8.603      6.160   1.397   0.1625  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 231.15  on 451  degrees of freedom
## Residual deviance: 226.78  on 449  degrees of freedom
## AIC: 232.78
## 
## Number of Fisher Scoring iterations: 5

The mean_pathogenicity variable is statistically significant at the 10% level with a p-value of 0.07. The negative estimate of 2.839 indicates that as a protein’s mutation risk (mean_pathogenicity) increases, its statistical probability of clinical trial success significantly decreases. 

pred_resp <- predict(logit_model, newdata = test, type = "response")
table(test$outcome, (pred_resp > 0.5)*1, dnn=c("Truth","Predicted"))
##      Predicted
## Truth   1
##     0   4
##     1 110

The model predicts 1/approved for all. There is no column 0 under predicted as the model did not classify any targets as 0/failure. This could be due to an imbalance. The set shows 110 successful targets and only 4 failed targets, meaning most of the historical test data consists of successful targets.

summary(pred_resp)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.8059  0.9261  0.9374  0.9322  0.9456  0.9597

The lowest probability the model assigned to a protein target is around 81%. Applying a standard cutoff of 0.5 means that observations get mapped to 1/Approved. The model is biasing predictions because the training dataset has many successful targets and contains few failures.

Addressing Imbalance

# 5 fold cv
set.seed(123)
folds <- createFolds(train$outcome_factor, k = 5, list = TRUE, returnTrain = FALSE)
# empty df to store out of fold predictions
oof_predictions <- data.frame()
for(i in 1:5) {
  val_indices <- folds[[i]]
  cv_train <- train[-val_indices, ]
  cv_val <- train[val_indices, ]
  fold_model <- glm(outcome ~ mean_pathogenicity + pathogenicity_sd, family = binomial, data = cv_train)
  fold_preds <- predict(fold_model, newdata = cv_val, type = "response")
  fold_results <- data.frame(Truth = cv_val$outcome, Predicted_Prob = fold_preds)
  oof_predictions <- rbind(oof_predictions, fold_results)
}

cv_roc <- roc(oof_predictions$Truth, oof_predictions$Predicted_Prob)
best_threshold <- coords(cv_roc, "best")$threshold
best_threshold
## [1] 0.9425365
# 80% training fit
final_logit_model <- glm(outcome ~ mean_pathogenicity + pathogenicity_sd, family = binomial, data = train)
plot(cv_roc, col = "#E6E6FA", main = "Cross-Validation ROC Curve", lwd = 2)
legend("bottomright", legend = paste("AUC =", round(auc(cv_roc), 3)), bty = "n", cex = 1.1)

# 20% test
test_resp <- predict(final_logit_model, newdata = test, type = "response")
table(test$outcome, (test_resp > best_threshold)*1, dnn=c("Truth","Predicted"))
##      Predicted
## Truth  0  1
##     0  3  1
##     1 71 39

Applying a strict threshold (0.9425) to the unseen test data successfully flags 3 out of 4 actual clinical failures. However, this cutoff misclassifies 71 out of 110 approvals as failures because their probabilities did not meet the cutoff requirement.

# test ROC
test_roc <- roc(test$outcome, test_resp)
plot(test_roc, col = "#E6E6FA", main = " Test ROC Curve", lwd = 2)
legend("bottomright", legend = paste("AUC =", round(auc(test_roc), 3)), bty = "n", cex = 1.1)

The test produces an AUC of 0.695, indicating a decent signal on the data split. However, the large gap between the 0.491 cross validation AUC and 0.695 test AUC shows instability in the data. The model is highly sensitive to how this data is partitioned. Performance metrics are volatile across data splits, likely because the data set only contains a handful of clinical failures.

GAM

A Generalized Additive Model (GAM) blends standard logistic regression with flexible smoothing splines to automatically detect and capture complex, non-linear relationships between the pathogenicity metrics and clinical outcomes.

protein_gam <- gam(outcome ~ s(mean_pathogenicity) + s(pathogenicity_sd), data = train, family = binomial())
summary(protein_gam)
## 
## Family: binomial 
## Link function: logit 
## 
## Formula:
## outcome ~ s(mean_pathogenicity) + s(pathogenicity_sd)
## 
## Parametric coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   2.6257     0.1906   13.77   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                       edf Ref.df Chi.sq p-value  
## s(mean_pathogenicity)   1      1  3.263  0.0709 .
## s(pathogenicity_sd)     1      1  1.950  0.1625  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.0104   Deviance explained = 1.89%
## UBRE = -0.48501  Scale est. = 1         n = 452

Effective degrees of freedom for all variables are 1. There are no non linear patterns or curves.

plot(protein_gam, shade=TRUE, seWithMean=TRUE, scale=0, pages = 1)

vis.gam(protein_gam, view=c("mean_pathogenicity","pathogenicity_sd"), theta= 140)

prob_test_gam <- predict(protein_gam, newdata = test, type = "response")
gam_roc <- roc(test$outcome, prob_test_gam)
auc(gam_roc)
## Area under the curve: 0.6955
gam_coords <- coords(gam_roc, "best", ret = c("threshold", "specificity", "sensitivity"))
table(test$outcome, (prob_test_gam > gam_coords$threshold)*1, dnn=c("Truth","Predicted"))
##      Predicted
## Truth  0  1
##     0  3  1
##     1 26 84

Generalized Additive Model (GAM) was used to find whether nonlinear relationships existed between pathogenicity features and clinical outcome. Both predictors had an effective degrees of freedom (EDF) of 1, indicating linear effects. Neither smooth term was statistically significant, and the model explained only 1.89% deviance. This suggests that GAM did not provide additional predictive value.

Random Forest

Random Forest is an ensemble machine learning method that builds a large collection of independent decision trees on random data subsets and aggregates their predictions to establish a robust classification boundary.

table(historical_data$outcome_factor)
## 
##   0   1 
##  36 530

The dataset contains 530 approvals and 36 terminations.

set.seed(123)
index <- sample(nrow(historical_data), 0.8 * nrow(historical_data))
train <- historical_data[index, ]
test <- historical_data[-index, ]
n_failures <- sum(train$outcome_factor == "0")

rf <- randomForest(outcome_factor ~ mean_pathogenicity + pathogenicity_sd, data = train, ntree = 500, sampsize = c(n_failures, n_failures), importance = TRUE)

print(rf)
## 
## Call:
##  randomForest(formula = outcome_factor ~ mean_pathogenicity +      pathogenicity_sd, data = train, ntree = 500, sampsize = c(n_failures,      n_failures), importance = TRUE) 
##                Type of random forest: classification
##                      Number of trees: 500
## No. of variables tried at each split: 1
## 
##         OOB estimate of  error rate: 37.39%
## Confusion matrix:
##     0   1 class.error
## 0   8  24   0.7500000
## 1 145 275   0.3452381

The model has an out of bag error of 37.39% with a 75% error rate on the terminated proteins.

rf_test_probs <- predict(rf, newdata = test, type = "prob")[, 2]
rf_test_roc <- roc(response = test$outcome, predictor = rf_test_probs, quiet = TRUE)
round(auc(rf_test_roc), 4)
## [1] 0.7136

The model achieved an AUC of 0.7136, indicating its moderate ability to discriminate between successful protein targets and clinical failures.

rf_test_preds <- ifelse(rf_test_probs > 0.5, 1, 0)
table(Truth = test$outcome, Predicted = rf_test_preds)
##      Predicted
## Truth  0  1
##     0  3  1
##     1 29 81

At a standard 0.5 probability threshold, the model flagged 3 out of 4 clinical failures in the test set, but misclassified 29 viable proteins as high-risk.

While the random forest model addressed imbalance by down sampling, its limited AUC of .7136 suggests that mean pathogenicity and variance alone are insufficient to define the complex threshold between clinical success and failure.

While the Random Forest is a sophisticated ensemble method, the Logistic Regression model is superior for this dataset because its linear approach provides a smoother boundary that better captures the relationship between protein pathogenicity and trial outcomes.

Boosting

Gradient Boosting is an ensemble machine learning method that improves prediction accuracy by combining many small decision trees. The outcome of this dataset is binary, whether a drug is approved or terminated, so a Bernoulli distribution is used. The model builds trees sequentially, with each new tree focusing on correcting the errors of the previous ones. As more trees are added, they work together to form a stronger predictive model that is generally more accurate than a single decision tree.

protein_boost <- gbm(outcome ~ mean_pathogenicity + pathogenicity_sd, data = train, distribution = "bernoulli", n.trees = 100, interaction.depth = 3, shrinkage = 0.05)       
summary(protein_boost)

##                                   var  rel.inf
## mean_pathogenicity mean_pathogenicity 52.76764
## pathogenicity_sd     pathogenicity_sd 47.23236

Chart shows how much each predictor contributes to the model’s ability to separate outcomes: approved vs termination.

par(mfrow = c(1, 2))
plot(protein_boost, i = "mean_pathogenicity")

plot(protein_boost, i = "pathogenicity_sd")

The jagged data is likely overfitted.


In this set the data imbalance is addressed.

# calculate 1s and 0s in the data.
majority_count <- sum(train$outcome == 1)
minority_count <- sum(train$outcome == 0)
weight_multiplier <- majority_count / minority_count

# assign 1 as weight for approval, and terminations for the heavier weight
train$case_weights <- ifelse(train$outcome == 0, weight_multiplier, 1)

set.seed(123)
protein_tuned_boost <- gbm(outcome ~ mean_pathogenicity + pathogenicity_sd, data = train, distribution = "bernoulli", n.trees = 500, interaction.depth = 1, shrinkage = 0.01, cv.folds = 5, weights = train$case_weights)

# find overfitting
par(mfrow = c(1, 1))
best_iter <- gbm.perf(protein_tuned_boost, method = "cv")

best_iter
## [1] 22

Plot shows how well model learns training data. Black line shows error drop as it goes through iterations. Green line shows performance on unseen data. Blue line shows that after around 22 trees, the model stops learning and memorizes noise. Stopping the model at 22 prevents overfitting from first run.

# partial dependence plot rerun
par(mfrow = c(1, 2))
plot(protein_tuned_boost, i = "mean_pathogenicity", n.trees = best_iter)

plot(protein_tuned_boost, i = "pathogenicity_sd", n.trees = best_iter)

Plot looks smoother on this run as model makes less complex / specific rules.

prob_test_tuned <- predict(protein_tuned_boost, newdata = test, n.trees = best_iter, type = "response")
tuned_roc <- roc(test$outcome, prob_test_tuned)
tuned_coords <- coords(tuned_roc, "best", ret = c("threshold", "specificity", "sensitivity"))
auc(tuned_roc)
## Area under the curve: 0.7295

Model achieved a strong AUC of .7

table(test$outcome, (prob_test_tuned > tuned_coords$threshold) * 1, dnn = c("Truth", "Predicted"))
##      Predicted
## Truth  0  1
##     0  3  1
##     1 40 70

The model identifies 3 out of 4 (75%) of the clinical terminations by acting conservatively, and results in 40 approved trials being flagged as false alarms. In a clinical research setting, the operational cost of manually reviewing these false alarms is small compared to the millions of dollars saved by preventing the wrong trials from advancing.

Neural Network

Neural Network is a computational model that passes input features through layers of interconnected nodes (neurons) to map complex, non-linear decision boundaries.

set.seed(123)
# random sample for 20% of rows in test set
test_index <- sample(nrow(historical_data), 0.20 * nrow(historical_data))
test <- historical_data[test_index, ]
remaining <- historical_data[-test_index, ]
# validation on 25% of the remaining 80% (20% of original data) 
val_index <- sample(nrow(remaining), 0.25 * nrow(remaining))
validation <- remaining[val_index, ]
train <- remaining[-val_index, ] # left over 60% for training
# training min max
maxs <- c(max(train$mean_pathogenicity), max(train$pathogenicity_sd))
mins <- c(min(train$mean_pathogenicity), min(train$pathogenicity_sd))
# training scale
train_scaled <- train
train_scaled$mean_pathogenicity <- (train$mean_pathogenicity - mins[1]) / (maxs[1] - mins[1])
train_scaled$pathogenicity_sd <- (train$pathogenicity_sd - mins[2]) / (maxs[2] - mins[2])
# validation scale
validation_scaled <- validation
validation_scaled$mean_pathogenicity <- (validation$mean_pathogenicity - mins[1]) / (maxs[1] - mins[1])
validation_scaled$pathogenicity_sd <- (validation$pathogenicity_sd - mins[2]) / (maxs[2] - mins[2])
# scale test set
test_scaled <- test
test_scaled$mean_pathogenicity <- (test$mean_pathogenicity - mins[1]) / (maxs[1] - mins[1])
test_scaled$pathogenicity_sd <- (test$pathogenicity_sd - mins[2]) / (maxs[2] - mins[2])

set.seed(123)
protein_nn <- neuralnet(outcome ~ mean_pathogenicity + pathogenicity_sd, data = train_scaled, hidden = c(3), algorithm = 'rprop+', linear.output = FALSE, likelihood = TRUE, stepmax = 1e6)

plot(protein_nn, rep = "best")

# validation
prob_nn <- compute(protein_nn, validation_scaled[, c("mean_pathogenicity", "pathogenicity_sd")])$net.result
val_roc <- roc(validation_scaled$outcome, prob_nn, quiet = TRUE)
pcut_nn <- coords(val_roc, "best", ret = "threshold")$threshold
# threshold
round(pcut_nn, 4)
## [1] 0.9543

Model showed a validation cutoff of .9543. Sample data is only classified as approval/1 when this strict cutoff is met.

# test data eval
prob_test <- compute(protein_nn, test_scaled[, c("mean_pathogenicity", "pathogenicity_sd")])$net.result
nn_roc <- roc(test_scaled$outcome, prob_test, quiet = TRUE)
# test predictions with validation threshold
pred_class <- ifelse(prob_test > pcut_nn, 1, 0)
table(Truth = test_scaled$outcome, Predicted = pred_class)
##      Predicted
## Truth  0  1
##     0  5  5
##     1 47 56
# oos test auc
auc(nn_roc)
## Area under the curve: 0.5291

The model achieved an out-of-sample test AUC of 0.53 and successfully identified 50% (5/10) of the clinical terminations. However, it flags 47 viable protein approvals as false alarms.

Practical Application

Based on the analysis, the Gradient Boosting Machine (GBM) was selected as the leading model for use. It demonstrated the highest overall discriminative power with an AUC of 0.7295, providing a strong baseline for evaluation.

active_trials <- final_ml_dataset |>
  filter(clinical_status == "In Clinical Trial (Currently)") |>
  arrange(accession_id)
# success probability prediction
active_trials$success_prob <- predict(protein_tuned_boost, newdata = active_trials, n.trees = best_iter, type = "response")
# success probability prediction
active_trials$predicted_status <- if_else(active_trials$success_prob > tuned_coords$threshold, "Expected Approval", "Termination/High Risk")

model <- active_trials |> 
  select(accession_id, mean_pathogenicity, pathogenicity_sd, success_prob, predicted_status) |> 
  arrange(success_prob)

kable(head(model))
accession_id mean_pathogenicity pathogenicity_sd success_prob predicted_status
P18146 0.6570661 0.2749205 0.4790503 Termination/High Risk
O14744 0.7018653 0.3306804 0.4825720 Termination/High Risk
O75581 0.6618153 0.3314799 0.4825720 Termination/High Risk
O75896 0.6728270 0.3271950 0.4825720 Termination/High Risk
P02538 0.6321627 0.3149018 0.4825720 Termination/High Risk
P05091 0.6627632 0.3216853 0.4825720 Termination/High Risk 
table(model$predicted_status)
## 
##     Expected Approval Termination/High Risk 
##                   496                   396

The target P18146 shows highest statistical risk in the active clinical trial dataset. This returned a predicted success probability of 47.91%. Based on the model, there were a total of 496 expected approvals and 396 expected termination and or high risk targets.

Strategic Value

Evaluating active pipelines prior to clinical trials provides an early stage screening ability for pharmaceutical companies and portfolio management. Flagging 396 targets as risky lets clinical development teams prioritize targets for stronger operational or biological testing. Identifying vulnerable targets helps mitigate inefficient expenses in late stage trials. This can allow biopharmaceutical organizations to optimize resource allocation toward pipelines that have a higher probability of success.

Reference