2  The Core Workflow of ML & Biological Pre-processing

In the transition from classical biostatistics to Machine Learning, the most critical shift is moving from explaining the past to predicting the future.

This means that, unlike traditional hypothesis testing where we focus on \(p\)-values and confidence intervals, in ML we prioritize out-of-sample performance, or well can we generalize to new, unseen data.

While a \(p\)-value tells us how unlikely a result is under the null hypothesis, an ML model’s utility is measured by its Generalization Error: how well it performs on a sample it has never seen before.

2.1 The ML Workflow Overview

A robust ML pipeline for biological data follows these non-negotiable steps:

1. Data Cleaning: Handling NAs and filtering low-variance features (e.g., genes with near-zero counts).
2. Data Spending: Splitting the data into training and testing sets to prevent data leakage.
3. Biological Pre-processing: Normalizing, filtering, and transforming biological data
4. Model Training: Fitting the algorithm to the training data.
5. Evaluation: Testing the “unseen” data to see if the model actually learned biology or just memorized noise.

2.2 {tidymodels}: A Unified Framework for ML in R

The {tidymodels} ecosystem provides a cohesive set of packages that streamline the entire ML workflow, from data pre-processing to model evaluation. It emphasizes a tidy, consistent syntax that integrates well with the broader {tidyverse}.

At its core, {tidymodels} is a unified framework for machine learning in R. Instead of learning different syntax for every single algorithm, you use one consistent language.

Think of it as a Lego set for data science: you snap different pieces together to build a complete pipeline.

To understand {tidymodels}, we just need to know the roles of its four main “workers”:

2.2.1 1. {rsample} (The Splitter)

Before you do anything, you need to set aside data to test your model later. {rsample} handles the “Training vs. Testing” split.

• The Goal: To make sure the model doesn’t just “memorize” the data it has already seen.

2.2.2 2. {recipes} (The Chef)

This is where you define your preprocessing steps. Just like a cooking recipe, it’s a list of instructions to get your raw data ready for the model.

• Common Steps: Filling in missing values (imputation), scaling numbers so they are on the same range (0 to 1), or converting text categories into numbers (dummy variables).

• Crucial Detail: The recipe defines the steps, but it doesn’t “cook” them until the model is actually run.

2.2.3 3. {parsnip} (The Interface)

In R, different models (like Random Forests vs. Linear Regression) often require completely different code styles. {parsnip} solves this by providing a unified interface.

• Example: You tell {parsnip} you want a “Random Forest,” and it handles the background translation to whatever specific R package (the “engine”) you want to use. You only have to learn one way to write it.

2.2.4 4. {workflows} (The Glue)

A Workflow bundles your {recipe} and your {parsnip} model together into a single object.

• It treats the entire process as one unit. When you want to predict new data, you just pass it to the workflow, and it automatically applies the same preprocessing steps before running the model.

2.2.5 Example

library(tidymodels)

# 1. Split the data
data_split <- initial_split(mtcars, prop = 0.8)
train_data <- training(data_split)

# 2. Define the Recipe (The Prep)
# "Predict mpg using all other variables, then scale them"
simple_recipe <- recipe(mpg ~ ., data = train_data) %>%
  step_normalize(all_predictors())

# 3. Define the Model (The Engine)
# "I want a linear regression model using the 'lm' engine"
lm_spec <- linear_reg() %>%
  set_engine("lm")

# 4. Create the Workflow (The Glue)
simple_workflow <- workflow() %>%
  add_recipe(simple_recipe) %>%
  add_model(lm_spec)

# 5. Fit the whole thing
final_model <- fit(simple_workflow, 
                   data = train_data)

2.3 TCGA Data

For our course, we will use the TCGA Cholangiocarcinoma (CHOL) cohort. We will use the {TCGAbiolinks} package to download and prepare the data.

###########
# Prepare our Data for Machine Learning with 
# TCGA Cholangiocarcinoma (CHOL) Cohort
library(TCGAbiolinks)
library(SummarizedExperiment)
library(tidyverse)

# 1. Query only the Cholangiocarcinoma cohort
query_chol <- GDCquery(
  project = "TCGA-CHOL",
  data.category = "Transcriptome Profiling",
  data.type = "Gene Expression Quantification",
  workflow.type = "STAR - Counts"
)
#  Load the data into a SummarizedExperiment object
chol_se <- GDCprepare(query_chol) 

|                                                    |  0%                      
|=                                                   |2.272727% ~4 s remaining  
|==                                                  |4.545455% ~3 s remaining  
|===                                                 |6.818182% ~2 s remaining  
|====                                                |9.090909% ~2 s remaining  
|=====                                               |11.36364% ~2 s remaining  
|=======                                             |13.63636% ~2 s remaining  
|========                                            |15.90909% ~4 s remaining  
|=========                                           |18.18182% ~3 s remaining  
|==========                                          |20.45455% ~3 s remaining  
|===========                                         |22.72727% ~3 s remaining  
|=============                                       | 25% ~3 s remaining       
|==============                                      |27.27273% ~2 s remaining  
|===============                                     |29.54545% ~2 s remaining  
|================                                    |31.81818% ~2 s remaining  
|=================                                   |34.09091% ~2 s remaining  
|==================                                  |36.36364% ~2 s remaining  
|====================                                |38.63636% ~2 s remaining  
|=====================                               |40.90909% ~2 s remaining  
|======================                              |43.18182% ~2 s remaining  
|=======================                             |45.45455% ~2 s remaining  
|========================                            |47.72727% ~2 s remaining  
|==========================                          | 50% ~2 s remaining       
|===========================                         |52.27273% ~2 s remaining  
|============================                        |54.54545% ~1 s remaining  
|=============================                       |56.81818% ~1 s remaining  
|==============================                      |59.09091% ~1 s remaining  
|===============================                     |61.36364% ~1 s remaining  
|=================================                   |63.63636% ~1 s remaining  
|==================================                  |65.90909% ~1 s remaining  
|===================================                 |68.18182% ~1 s remaining  
|====================================                |70.45455% ~1 s remaining  
|=====================================               |72.72727% ~1 s remaining  
|=======================================             | 75% ~1 s remaining       
|========================================            |77.27273% ~1 s remaining  
|=========================================           |79.54545% ~1 s remaining  
|==========================================          |81.81818% ~1 s remaining  
|===========================================         |84.09091% ~0 s remaining  
|============================================        |86.36364% ~0 s remaining  
|==============================================      |88.63636% ~0 s remaining  
|===============================================     |90.90909% ~0 s remaining  
|================================================    |93.18182% ~0 s remaining  
|=================================================   |95.45455% ~0 s remaining  
|==================================================  |97.72727% ~0 s remaining  
|====================================================|100% ~0 s remaining       
|====================================================|100%                      Completed after 3 s 

#######

# 2. Download and prepare the data
# Extract fpkm using 
fpkm_data <- assay(chol_se, "fpkm_unstrand")

# 3. Quick Metadata Link
metadata <- as.data.frame(colData(chol_se)) %>%
  select(barcode, sample_type)
head(metadata)

                                                  barcode         sample_type
TCGA-W5-AA39-01A-11R-A41I-07 TCGA-W5-AA39-01A-11R-A41I-07       Primary Tumor
TCGA-3X-AAVB-01A-31R-A41I-07 TCGA-3X-AAVB-01A-31R-A41I-07       Primary Tumor
TCGA-W5-AA2R-11A-11R-A41I-07 TCGA-W5-AA2R-11A-11R-A41I-07 Solid Tissue Normal
TCGA-W5-AA38-01A-11R-A41I-07 TCGA-W5-AA38-01A-11R-A41I-07       Primary Tumor
TCGA-W5-AA2G-01A-11R-A41I-07 TCGA-W5-AA2G-01A-11R-A41I-07       Primary Tumor
TCGA-W5-AA2Q-11A-11R-A41I-07 TCGA-W5-AA2Q-11A-11R-A41I-07 Solid Tissue Normal

#######

# 4. ML Preparation (Transpose)
# Often in ML, TPM or FPKM values are used directly
# We need the samples as rows and genes as columns
# We filter for low-expression genes to remove noise
# For fpkm, a common threshold is > 5 in a certain % of samples. 
keep <- rowSums(fpkm_data > 5) >= ncol(fpkm_data) * 0.2 # Keep genes expressed in at least 20% of samples
fpkm_filtered <- fpkm_data[keep, ]
dim(fpkm_filtered) # Check dimensions after filtering

[1] 8646   44

## Combine with metadata
df_ml_fpkm <- as.data.frame(t(fpkm_filtered)) %>%
  rownames_to_column("barcode") %>%
  inner_join(metadata, by = "barcode") %>%
  select(-barcode) %>%
  mutate(sample_type = factor(sample_type))
dim(df_ml_fpkm) # Check final dimensions

[1]   44 8647

2.4 The Spending Plan: Data Splitting Strategies

In biology, samples are precious and often limited (\(n < 100\)). This scarcity makes the “Data Spending” phase the most high-stakes part of your project. If you use your data too aggressively during training, you will overfit; if you save too much for testing, your model will be too weak to learn the underlying biology. Finally, if you “peek” at the test set during pre-processing, you will introduce Data Leakage¹, leading to overly optimistic performance estimates.

2.4.1 The Initial Split: Training vs. Testing

To ensure that our model’s performance is a true reflection of its ability to generalize, we must carefully partition our dataset at the very beginning of our analysis. For that, we divide our dataset into two primary components (Figure 2.1):

1. The Training Set: This subset is used to train the model. All pre-processing steps (normalization, feature selection, imputation) must be derived solely from this set. Often, we use 70% to 80% of the data for training.

2. The Test Set: This subset is held out and only used once at the very end to evaluate the model’s performance. It must remain completely unseen during training and pre-processing.

Note: Never perform normalization, feature selection, or imputation on the entire dataset before splitting. If you calculate the mean expression of a gene using all samples, the training set now “knows” something about the distribution of the test set. This is Data Leakage, and it leads to artificially inflated (and ultimately false) performance metrics.

Figure 2.1: Data Spending: Training vs. Testing Sets

The easiest way to implement this split in R is using the initial_split() function from the {rsample} package, which is part of the {tidymodels} ecosystem.

# Split the data into training and testing sets
library(tidymodels)
set.seed(12345) # For reproducibility
data_split <- initial_split(df_ml_fpkm, 
                            prop = 0.8)
train_data_chol <- training(data_split)
test_data_chol  <- testing(data_split)

cat("Training samples:", nrow(train_data_chol), "\nTesting samples:", nrow(test_data_chol))

Training samples: 35 
Testing samples: 9

cat("Proportion of sample types in training set:\n")

Proportion of sample types in training set:

print(prop.table(table(train_data_chol$sample_type)))

      Primary Tumor Solid Tissue Normal 
                0.8                 0.2 

cat("Proportion of sample types in testing set:\n")

Proportion of sample types in testing set:

print(prop.table(table(test_data_chol$sample_type)))

      Primary Tumor Solid Tissue Normal 
          0.7777778           0.2222222 

2.4.2 Stratified Sampling

In several biological contexts, the classes we are trying to predict are imbalanced, which can lead to misleading performance metrics if not handled properly. Menaing that one class (e.g., healthy controls) may vastly outnumber another (e.g., patients with a rare disease). If we randomly split the data, we risk creating training and testing sets that do not accurately represent the overall class distribution. This can lead to models that perform well on the majority class but poorly on the minority class, which is often of greater interest in biological studies.

To address this, we use Stratified Sampling to ensure that the proportion of the outcome (e.g., disease status) is preserved in both the training and testing sets.

In R, we will use the same function as before, but with a new argument.

# Split the data into training and testing sets
library(tidymodels)
set.seed(42) # For reproducibility
data_split <- initial_split(df_ml_fpkm, 
                            prop = 0.8, 
                            strata = sample_type)
train_data_chol <- training(data_split)
test_data_chol  <- testing(data_split)

cat("Training samples:", nrow(train_data_chol), "\nTesting samples:", nrow(test_data_chol))

Training samples: 35 
Testing samples: 9

cat("Proportion of sample types in training set:\n")

Proportion of sample types in training set:

print(prop.table(table(train_data_chol$sample_type)))

      Primary Tumor Solid Tissue Normal 
                0.8                 0.2 

cat("Proportion of sample types in testing set:\n")

Proportion of sample types in testing set:

print(prop.table(table(test_data_chol$sample_type)))

      Primary Tumor Solid Tissue Normal 
          0.7777778           0.2222222 

2.5 Pre-processing for Biological Data

Biological data is rarely “model-ready.” We must address two major issues: Scale and Skewness.

2.5.1 Scaling and Centering

Most ML algorithms (like SVM or Lasso) use distance-based metrics. If Gene A has expression values in the thousands and Gene B in the decimals, the model will unfairly prioritize Gene A.

2.5.2 Near-Zero Variance (NZV)

In RNA-seq, many genes show little to no variation across samples. These are “noise” for a predictive model and increase the “Curse of Dimensionality.”

2.5.3 The recipe Approach

In {tidymodels}, we define a “recipe” – a blueprint of transformations.

library(tidymodels)
# Define a recipe for pre-processing
# We will scale, center, and remove near-zero variance predictors

# Define a recipe for pre-processing
# Focus specifically on numeric predictors to avoid scaling the outcome factor
ml_recipe <- recipe(sample_type ~ ., data = train_data_chol) %>%
  step_zv(all_predictors()) %>%            # Remove zero variance (constant) genes
  step_nzv(all_numeric_predictors()) %>%   # Remove genes with very little variation
  step_normalize(all_numeric_predictors()) # Z-score transformation (mean=0, sd=1)

# Statistical Note for Students: 
# We prep() using ONLY train_data_chol to ensure the mean and SD 
# are not influenced by our test set (avoiding "Data Leakage").
ml_prep <- prep(ml_recipe, training = train_data_chol)

# Apply (bake) to both sets
train_processed <- bake(ml_prep, new_data = NULL) # NULL defaults to the training data
test_processed  <- bake(ml_prep, new_data = test_data_chol)

# Check dimensions: How many "noisy" genes did we drop?
cat("Original Training Dim:", dim(train_data_chol), "\n")

Original Training Dim: 35 8647 

cat("Processed Training Dim:", dim(train_processed), "\n")

Processed Training Dim: 35 8647 

# How about the test set?
cat("Original Testing Dim:", dim(test_data_chol), "\n")

Original Testing Dim: 9 8647 

cat("Processed Testing Dim:", dim(test_processed), "\n")

Processed Testing Dim: 9 8647 

With our data now properly split and pre-processed, we still have an issue: our training set is small. To build a robust model, we need to maximize the utility of our training data without overfitting. This is where Cross-Validation comes into play.

2.6 Cross-Validation: Maximizing Training Data Utility

Cross-Validation (CV) is a resampling technique used to evaluate ML models on a limited data sample. The most common form is k-Fold Cross-Validation. In k-Fold CV, the training data is divided into k subsets (or “folds”). The model is trained on k-1 folds and validated on the remaining fold. This process is repeated k times, with each fold serving as the validation set once (Figure 2.2). The final performance metric is averaged across all folds.

Term	Definition
Fold	A subset of the training data.
Resampling	Repeating the split process.
Hyperparameter	Settings of the model (not learned from data).

Figure 2.2: K-Fold Cross-Validation: The Internal Loop

To implement k-Fold CV in R, we use the vfold_cv() function from {rsample}.

# Create stratified 5-fold cross-validation
set.seed(123)
cv_folds <- vfold_cv(train_data_chol, 
                     v = 5)

cv_folds

#  5-fold cross-validation 
# A tibble: 5 × 2
  splits         id   
  <list>         <chr>
1 <split [28/7]> Fold1
2 <split [28/7]> Fold2
3 <split [28/7]> Fold3
4 <split [28/7]> Fold4
5 <split [28/7]> Fold5

2.6.1 Stratified K-Fold: Ensuring Biological Representation

As we have discussed in Section Section 2.4, biological datasets are often imbalanced. If we use standard K-Fold CV, we risk an iteration where the “Assessment” fold accidentally contains zero cases of a specific disease subtype. Stratified K-Fold ensures that each internal fold preserves the same ratio of classes as the original training set (Figure 2.3).

Figure 2.3: Stratified K-Fold Cross-Validation

In R, implementing this is as simple as adding the strata argument to your resampling function:

# Create stratified 5-fold cross-validation
set.seed(12345)
cv_folds <- vfold_cv(train_data_chol, 
                     v = 5, 
                     strata = sample_type)
cv_folds

#  5-fold cross-validation using stratification 
# A tibble: 5 × 2
  splits         id   
  <list>         <chr>
1 <split [27/8]> Fold1
2 <split [27/8]> Fold2
3 <split [28/7]> Fold3
4 <split [29/6]> Fold4
5 <split [29/6]> Fold5

###### Before start, let us do a quick t-test just to check

# Perform t-tests for each gene
t_test_results <- train_processed %>%
  pivot_longer(cols = -sample_type, 
               names_to = "Gene", 
               values_to = "Expression") %>%
  group_by(Gene) %>%
  summarize(
    p_value = t.test(Expression ~ sample_type)$p.value,
    mean_tumor = mean(Expression[sample_type == "Primary Tumor"]),
    mean_normal = mean(Expression[sample_type == "Solid Tissue Normal"]),
    diff = mean_tumor - mean_normal
  ) %>%
  arrange(p_value)

# Visualize using Esquisse :) 

ggplot(t_test_results) +
  aes(x = diff, y = p_value) +
  geom_point(colour = "#112446") +
  theme_minimal()

2.7 Summary and Concluding Remarks

In this Chapter, we have established the “infrastructure” of a machine learning project. We moved beyond simple data loading to a rigorous pipeline that respects biological complexity and statistical integrity.

2.7.1 Key Takeaways

• The Mindset Shift: We moved from \(p\)-values (inference) to Generalization Error (prediction).

• Data Spending: We learned that the Test Set is a “Locked Vault” that must remain untouched until the very end to prevent Data Leakage.

• Feature Engineering: Using {tidymodels} recipes, we automated the removal of Near-Zero Variance (NZV) genes and normalized high-throughput counts to prevent feature scale bias.

• Validation: We implemented Cross-Validation to maximize our small biological sample size (\(n\)) while maintaining a “mock” testing environment.

2.7.2 Final Remarks

You now have a “processed” dataset (train_processed) and a validation strategy (cv_folds). However, we haven’t actually looked at our data yet. In the next chapter, we will explore Unsupervised Learning. We will use Dimensionality Reduction (PCA, tSNE, UMAP) to see if our biological groups (Cases vs. Controls) naturally separate before we ever try to “force” a model to learn them.

---

1. Data Leakage occurs when information from outside the training dataset is used to create the model. This can lead to overly optimistic performance estimates because the model has effectively “seen” parts of the test data during training.