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)
###########
# 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)
#######
# 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)
#######
# 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
## 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
#######
# 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))
cat("Proportion of sample types in training set:\n")
print(prop.table(table(train_data_chol$sample_type)))
cat("Proportion of sample types in testing set:\n")
print(prop.table(table(test_data_chol$sample_type)))
#######
# 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))
cat("Proportion of sample types in training set:\n")
print(prop.table(table(train_data_chol$sample_type)))
cat("Proportion of sample types in testing set:\n")
print(prop.table(table(test_data_chol$sample_type)))
### Just to check the mean and sd of
### the expression values for each gene
### in the training set
train_data_chol %>%
pivot_longer(cols = -sample_type,
names_to = "Gene",
values_to = "Expression") %>%
group_by(Gene, sample_type) %>%
summarize(mean_expression = mean(Expression),
sd = sd(Expression))
#######
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")
cat("Processed Training Dim:", dim(train_processed), "\n")
# How about the test set?
cat("Original Testing Dim:", dim(test_data_chol), "\n")
cat("Processed Testing Dim:", dim(test_processed), "\n")
# Create stratified 5-fold cross-validation
set.seed(123)
cv_folds <- vfold_cv(train_data_chol,
v = 5)
cv_folds
# Create stratified 5-fold cross-validation
set.seed(12345)
cv_folds <- vfold_cv(train_data_chol,
v = 5,
strata = sample_type)
cv_folds
###### 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()
####### Sec 3: Unsupervised
library(embed)
# Define the PCA Recipe
pca_rec <- recipe(sample_type ~ .,
data = train_processed) %>%
step_pca(all_predictors(), num_comp = 5)
# Train the recipe on the CHOL training set
pca_estimates <- prep(pca_rec)
# Extract the coordinates for visualization
pca_plot_data <- juice(pca_estimates)
# Visualization of PC1 vs PC2
pca_plot_data %>%
ggplot() +
aes(x = PC1, y = PC2, color = sample_type) +
geom_point(alpha = 0.7, size = 2) +
scale_color_manual(values = c("Primary Tumor" = "#2c7bb6", "Solid Tissue Normal" = "#d7191c")) +
theme_minimal() +
labs(title = "Principal Component Analysis: TCGA-CHOL",
x = paste0("PC1 (", round(pca_estimates$steps[[1]]$res$sdev[1]^2 / sum(pca_estimates$steps[[1]]$res$sdev^2)*100, 1), "%)"),
y = paste0("PC2 (", round(pca_estimates$steps[[1]]$res$sdev[2]^2 / sum(pca_estimates$steps[[1]]$res$sdev^2)*100, 1), "%)"),
color = "Sample Type") +
theme(legend.position = "bottom")
library(Rtsne)
# Prepare data for t-SNE
# Assuming 'train_data_chol' is a data frame with the sample type in the first column and gene expression in the remaining columns
# We will use the gene expression data for t-SNE
tsne_data <- train_processed %>%
select(-sample_type) %>%
as.matrix()
# Run t-SNE
# Set a random seed for reproducibility
set.seed(123)
tsne_result <- Rtsne(tsne_data,
perplexity = 1, # Play with the perplexity parameter to see how it affects the visualization (try 0.1, 1, 5, 10, 20)
verbose = TRUE,
max_iter = 1000)
# Create a data frame for plotting
tsne_plot_data <- data.frame(
X = tsne_result$Y[, 1],
Y = tsne_result$Y[, 2],
sample_type = train_data_chol$sample_type
)
# Visualization of t-SNE results
# Using ggplot2 for visualization
tsne_plot_data %>%
ggplot() +
aes(x = X, y = Y, color = sample_type) +
geom_point(alpha = 0.7, size = 2) +
scale_color_manual(values = c("Primary Tumor" = "#2c7bb6",
"Solid Tissue Normal" = "#d7191c")) +
theme_minimal() +
labs(title = "t-SNE: TCGA-CHOL",
x = "t-SNE Dimension 1",
y = "t-SNE Dimension 2",
color = "Sample Type") +
theme(legend.position = "bottom")
library(tidymodels)
library(embed)
# Define UMAP recipe
umap_rec <- recipe(sample_type ~ ., data = train_processed) %>%
step_umap(
all_predictors(),
num_comp = 2, # Plotting in 2D, but try 3D as well
neighbors = 5, # Try different values (5, 10, 15) to see how it affects the visualization
min_dist = 0.1 # Try different values (0.1, 0.5, 0.9) to see how it affects the clustering of points
)
# Train recipe
umap_estimates <- prep(umap_rec)
# Extract coordinates
umap_plot_data <- juice(umap_estimates)
# Visualization
umap_plot_data %>%
ggplot() +
aes(x = UMAP1, y = UMAP2, color = sample_type) +
geom_point(alpha = 0.7, size = 2) +
scale_color_manual(values = c("Primary Tumor" = "#2c7bb6",
"Solid Tissue Normal" = "#d7191c")) +
theme_minimal() +
labs(title = "UMAP Projection: TCGA-CHOL",
x = "UMAP 1",
y = "UMAP 2",
color = "Sample Type") +
theme(legend.position = "bottom")
# Compute distance matrix
# A common choice for gene expression data is the Euclidean distance
distance_matrix <- dist(
train_processed %>%
select(-sample_type),
method = "minkowski") # different distance metrics can be tried (euclidean, manhattan, minkowski, etc.)
## Euclidean is used for continuous data, while Manhattan can be more robust to outliers.
## Minkowski is a generalization of both
# Perform hierarchical clustering
# We can use the complete linkage method, which considers the maximum distance between points in different clusters
hc <- hclust(distance_matrix,
method = "complete")
# Plot the dendrogram
plot(hc,
labels = train_processed$sample_type,
main = "Hierarchical Clustering Dendrogram",
xlab = "",
sub = "")
# Determine the optimal number of clusters using the Elbow Method
# We will compute the total within-cluster sum of squares for a range of K values
wss <- sapply(1:10, function(k) {
kmeans(train_processed %>%
select(-sample_type),
centers = k,
nstart = 25)$tot.withinss
})
# Plot the Elbow Method
plot(1:10, wss,
type = "b",
las = 1,
pch = 19,
xlab = "Number of clusters K",
ylab = "Total within-clusters sum of squares")
# Perform K-means clustering with K=2
set.seed(123)
kmeans_result <- kmeans(train_processed %>%
select(-sample_type),
centers = 2,
nstart = 25)
cluster <- as.factor(kmeans_result$cluster)
# Visualize the clusters
# We can use PCA to visualize the clusters in a 2D space
pca_plot_data %>%
mutate(cluster = cluster) %>%
ggplot() +
aes(x = PC1, y = PC2, color = cluster, shape = sample_type) +
geom_point(alpha = 0.7, size = 2) +
scale_color_manual(values = c("1" = "#2c7bb6",
"2" = "#d7191c")) +
theme_minimal() +
labs(title = "K-Means Clustering (K=2) on PCA Projection",
x = paste0("PC1 (", round(pca_estimates$steps[[1]]$res$sdev[1]^2 / sum(pca_estimates$steps[[1]]$res$sdev^2)*100, 1), "%)"),
y = paste0("PC2 (", round(pca_estimates$steps[[1]]$res$sdev[2]^2 / sum(pca_estimates$steps[[1]]$res$sdev^2)*100, 1), "%)"),
color = "Cluster") +
theme(legend.position = "bottom")
library(dbscan)
# Prepare data for DBSCAN
# We will use the PCA-reduced data for DBSCAN to reduce computational complexity
dbscan_data <- pca_plot_data %>%
select(PC1, PC2) %>%
as.matrix()
# Run DBSCAN
# Set eps to a value that captures the local density of points and minPts to a value that reflects the expected cluster size
dbscan_result <- dbscan(dbscan_data,
eps = 20,
minPts = 5,
borderPoints = T
)
dbscan_result
# Add cluster assignments to the original data
pca_plot_data$cluster <- as.factor(dbscan_result$cluster)
# Visualize the DBSCAN clusters
pca_plot_data %>%
ggplot() +
aes(x = PC1, y = PC2, color = cluster) +
geom_point(alpha = 0.7, size = 2) +
scale_color_manual(values = c("0" = "#2c7bb6", # Cluster 0 (noise)
"1" = "#d7191c", # Cluster 1
"2" = "#fdae61", # Cluster 2
"3" = "#abdda4", # Cluster 3
"4" = "#2b83ba")) + # Cluster 4
theme_minimal() +
labs(title = "DBSCAN Clustering on PCA Projection",
x = paste0("PC1 (", round(pca_estimates$steps[[1]]$res$sdev[1]^2 / sum(pca_estimates$steps[[1]]$res$sdev^2)*100, 1), "%)"),
y = paste0("PC2 (", round(pca_estimates$steps[[1]]$res$sdev[2]^2 / sum(pca_estimates$steps[[1]]$res$sdev^2)*100, 1), "%)"),
color = "Cluster") +
theme(legend.position = "bottom")
####### Supervised
library(rpart)
library(rpart.plot)
### Create a mini Train, so it does not take a long time to run
main_genes <- t_test_results %>%
mutate(p_adj = p.adjust(p_value, method = "BH")) %>%
filter(p_adj > 0.8) %>%
head(50) %>%
pull(Gene)
train_processed_mini <- train_processed %>%
select(sample_type, main_genes)
# 1. Prepare the Data.
# We have already done it ;)
# 2. Factorize
train_processed_mini$sample_type <- as.factor(train_processed_mini$sample_type)
train_processed_mini$sample_type <- as.factor(train_processed_mini$sample_type)
# 3. Fit Model with Complexity Control
# cp (complexity parameter) helps prevent the tree from getting too "wild"
tree_model <- rpart(sample_type ~ .,
data = train_processed_mini,
method = "class",
control = rpart.control(cp = 0.001,
minsplit = 2))
# 4. Plot
# 'type = 5' shows the split labels clearly; 'extra = 104' shows probabilities and percentages
rpart.plot(tree_model,
type = 5,
extra = 104,
box.palette = "RdYlGn",
shadow.col = "gray",
main = "Decision Tree (CHOL)")
#######
# 1. Define the Recipe (Pre-processing)
# We select our outcome and predictors, then filter for the top most variable genes
tree_recipe <- recipe(sample_type ~ .,
data = train_processed_mini)
full_recipe <- recipe(sample_type ~ .,
data = train_processed )
# 2. Define the Model Specification
# Here we specify 'rpart' as the engine for a classification task
tree_spec <- decision_tree(cost_complexity = 0.01,
tree_depth = 5,
min_n = 2
) %>%
set_engine("rpart") %>%
set_mode("classification")
# 3. Create a Workflow
# This bundles the recipe and the model together
tree_workflow <- workflow() %>%
add_recipe(tree_recipe) %>%
add_model(tree_spec)
# 4. Fit the model
tree_fit <- tree_workflow %>%
fit(data = train_processed)
# 5. Extract the fitted model to plot the tree
extract_fit_engine(tree_fit) %>%
rpart.plot(type = 5,
extra = 104,
box.palette = "RdYlGn",
main = "Decision Tree: CHOL Dataset")
# Make predictions on the test dataset
test_processed$sample_type <- as.factor(test_processed$sample_type)
predictions <- predict(tree_model, newdata = test_processed, type = "class")
# Evaluate the accuracy of the model
# Create a confusion matrix
require(caret)
confusionMatrix(predictions, test_processed$sample_type)
#######
library(randomForest)
set.seed(12345) # For reproducibility
# Fit a random forest model
#
rf_model <- randomForest(sample_type ~ .,
mtry = 10, # Number of variables randomly sampled as candidates at each split,
data = train_processed,
maxnodes = 5,
nodesize = 2,
ntree = 100) # Number of trees to grow
# Print the random forest model
rf_model
#######
# Plot the random forest model
plot(rf_model)
#######
# Get feature importance
importance(rf_model) %>%
as.data.frame() %>%
arrange(desc(MeanDecreaseGini)) %>%
head(10) %>%
knitr::kable(caption = "Top 10 Most Important Features in the Random Forest Model")
#######
library(ggplot2)
imp_df <- importance(rf_model) %>%
as.data.frame() %>%
tibble::rownames_to_column("Gene") %>%
arrange(desc(MeanDecreaseGini)) %>%
head(20)
ggplot(imp_df) +
aes(x = reorder(Gene, MeanDecreaseGini),
y = MeanDecreaseGini) +
geom_point(size = 3,
color = "steelblue") +
geom_segment(aes(x = Gene,
xend = Gene,
y = 0,
yend = MeanDecreaseGini),
color = "skyblue") +
coord_flip() +
labs(title = "Variable Importance (Random Forest)",
x = "Genes",
y = "Mean Decrease in Gini Index") +
theme_minimal()
#######
# Look at the structure of the 1st tree
# k = 1 is the tree number
# labelVar = TRUE ensures it uses gene names instead of index numbers
randomForest::getTree(rf_model,
k = 1,
labelVar = TRUE) %>%
head(10) # Showing the first 10 nodes
#######
library(tidymodels)
set.seed(123) # For reproducibility
## The data step is the same as previously, we need only to update our model specification and workflow
# 1. Define the Random Forest Specification
# We'll grow 100 trees (trees = 100)
# mtry is the number of genes randomly sampled at each split
rf_spec <- rand_forest(
mtry = tune(), # We can tune this later
trees = 100,
min_n = 2
) %>%
set_engine("ranger",
importance = "impurity") %>%
set_mode("classification")
# 2. Update the Workflow
# We reuse the 'tree_recipe' from the previous section
rf_workflow <- workflow() %>%
add_recipe(full_recipe) %>% # Reusing the recipe
add_model(rf_spec %>% finalize_model(list(mtry = 2))) # Setting a starting mtry
# 3. Fit the Model
rf_fit <- rf_workflow %>%
fit(data = train_processed)
# 4. Extract the fitted model to get feature importance
rf_fit %>%
extract_fit_parsnip() %>%
vip::vip(num_features = 20,
geom = "point") +
theme_minimal() +
labs(title = "Random Forest: Top 20 Gene Biomarkers",
x = "Importance (Mean Decrease in Gini)",
y = "Genes")
#######
# Make predictions on the test dataset
# We will use our test data here
predictions <- predict(rf_model, newdata = test_processed)
# Evaluate the accuracy of the model
# Create a confusion matrix
require(caret)
confusionMatrix(predictions, test_processed$sample_type)
#######
# Make predictions on the test dataset
# We will use our test data here
predictions <- predict(rf_model, newdata = test_processed)
# Evaluate the accuracy of the model
# Create a confusion matrix
require(caret)
confusionMatrix(predictions, test_processed$sample_type)
#######
## AUC and AUC-PR
require(ROCR)
# Get predicted probabilities for the positive class
pred_prob <- predict(rf_model,
newdata = test_processed,
type = "prob")[, 1] # Tumour Prob is the First Column
sample_type <- ifelse(test_processed$sample_type == "Primary Tumor", 1, 0) # Convert to binary labels (1 for Tumor, 0 for Normal)
plot(pred_prob, sample_type,
xlab = "Predicted Probability of Tumor",
ylab = "Actual Sample Type (1=Tumor, 0=Normal)",
main = "Predicted Probabilities vs Actual Labels")
# Create a prediction object
pred <- ROCR::prediction(pred_prob, sample_type)
# Calculate AUC
auc <- ROCR::performance(pred, measure = "auc")@y.values[[1]]
auc
#######
# Calculate AUC-PR
auc_pr <- performance(pred, measure = "aucpr")@y.values[[1]]
auc_pr
#######
# plot the Precision-Recall curve
pr_perf <- performance(pred, measure = "prec", x.measure = "rec")
plot(pr_perf,
col = "red",
lwd = 2,
main = "Precision-Recall Curve",
xlab = "Recall",
ylab = "Precision")
#######
library(xgboost)
set.seed(12345) # For reproducibility
# Prepare the data for xgboost
# We need to convert the data into a matrix format and the labels into a numeric format
train_matrix <- train_processed_mini %>%
select(-sample_type) %>% # Exclude the label
as.matrix()
train_labels <- as.numeric(train_processed_mini$sample_type) - 1 # Convert to numeric (0 and 1)
# Fit a gradient boosting model
gb_model <- xgboost(data = train_matrix,
label = as.factor(train_labels),
nrounds = 100, # Number of boosting rounds
objective = "binary:logistic") # Binary classification
# Print the gradient boosting model
print(gb_model)
#######
# Get feature importance
importance_matrix <- xgb.importance(feature_names = colnames(train_matrix), model = gb_model)
importance_matrix %>%
head(10) %>%
knitr::kable(caption = "Top 10 Most Important Features in the Gradient Boosting Model")
#######
library(tidymodels)
library(xgboost)
library(vip)
set.seed(123)
# 1. Define the Gradient Boosting Specification
# Note: XGBoost usually needs many more trees than a Random Forest
# because it learns slowly (boosting vs bagging).
gb_spec <- boost_tree(
trees = 100,
tree_depth = 3,
learn_rate = 0.1
) %>%
set_engine("xgboost") %>%
set_mode("classification")
# 2. Update the Workflow (reusing your tree_recipe)
gb_workflow <- workflow() %>%
add_recipe(full_recipe) %>%
add_model(gb_spec)
# 3. Fit the Model
gb_fit <- gb_workflow %>%
fit(data = train_processed)
# 4. Extract and Plot Importance
# In tidymodels, you simply pass the parsnip object directly to vip().
# It will automatically interface with xgboost and retrieve the feature names.
gb_fit %>%
extract_fit_parsnip() %>%
vip(num_features = 20,
geom = "point",
aesthetics = list(color = "midnightblue", size = 3)) +
theme_minimal() +
labs(title = "Gradient Boosting: Top 20 Gene Biomarkers",
subtitle = "Importance calculated via Gain",
x = "Importance",
y = "Genes")
#######
# Prepare the test data for xgboost
test_matrix <- as.matrix(test_processed %>% select (-sample_type)) # Exclude the label
test_labels <- as.numeric(test_processed$sample_type) - 1 # Convert to numeric (0 and 1)
# Make predictions on the test dataset
predictions <- predict(gb_model, newdata = test_matrix)
# Convert predicted probabilities to class labels (0 or 1)
predicted_labels <- ifelse(predictions > 0.5, 1, 0 ) %>%
as.factor()
# Evaluate the accuracy of the model
caret::confusionMatrix(predicted_labels, as.factor(test_labels))
#######
## AUC and AUC-PR
## Get predicted probabilities for the positive class
pred_prob <- predict(gb_model,
newdata = test_matrix) # Predicted probabilities for the positive class
sample_type <- test_labels # Binary labels (1 for Tumor, 0 for Normal)
# Create a prediction object
pred <- ROCR::prediction(pred_prob, sample_type)
# Calculate AUC
auc <- ROCR::performance(pred,
measure = "auc")@y.values[[1]]
auc
#######
# plot the ROC curve
roc_perf <- ROCR::performance(pred, measure = "tpr", x.measure = "fpr")
plot(roc_perf, col = "blue",
lwd = 2, main = "ROC Curve",
xlab = "False Positive Rate",
ylab = "True Positive Rate")
#######
# Calculate AUC-PR
auc_pr <- ROCR::performance(pred,
measure = "aucpr")@y.values[[1]]
auc_pr
# plot the Precision-Recall curve
pr_perf <- ROCR::performance(pred,
measure = "prec",
x.measure = "rec")
plot(pr_perf,
col = "red",
lwd = 2,
main = "Precision-Recall Curve",
xlab = "Recall",
ylab = "Precision")
#######
# Fit a linear regression model
linear_model <- lm(ENSG00000001629.10 ~ sample_type, data = train_processed)
# Print the linear regression model
summary(linear_model)
### Ok, but I don't want to do one gene at a time, or put into a for loop.
require(broom)
fit = train_processed %>%
pivot_longer(cols = -sample_type,
names_to = "Gene",
values_to = "Expression") %>%
group_by(Gene) %>%
do(tidy(lm(Expression ~ sample_type, data = .))) %>%
filter(term != "(Intercept)") %>%
mutate(p_adj = p.adjust(p.value, method = "BH"))
fit %>%
ggplot() +
aes(x = estimate, y = -log10(p_adj)) +
geom_point(color = "steelblue") +
theme_minimal() +
labs(title = "Linear Regression: Effect of Sample Type on Gene Expression",
x = "Estimated Effect (Coefficient)",
y = "log10 Adjusted P-value (BH)") +
theme(legend.position = "bottom")
#######
library(tidymodels)
# 1. Define the Recipe
regression_recipe <- recipe(ENSG00000001629.10 ~ sample_type, data = train_processed)
regression_spec <- linear_reg() %>%
set_engine("lm") %>%
set_mode("regression")
# 3. Create a Workflow
regression_workflow <- workflow() %>%
add_recipe(regression_recipe) %>%
add_model(regression_spec)
# 4. Fit the model
regression_fit <- regression_workflow %>%
fit(data = train_processed)
# 5. Print the model summary
extract_fit_parsnip(regression_fit)
#######
## We start by making predictions on the test dataset
predictions <- predict(linear_model,
newdata = test_processed)
# Calculate Mean Squared Error (MSE)
mse <- mean((predictions - test_processed$ENSG00000001629.10)^2)
mse
# Calculate Mean Absolute Error (MAE)
mae <- mean(abs(predictions - test_processed$ENSG00000001629.10))
mae
# Calculate R-squared
ss_total <- sum((test_processed$ENSG00000001629.10 - mean(test_processed$ENSG00000001629.10))^2)
ss_residual <- sum((test_processed$ENSG00000001629.10 - predictions)^2)
r_squared <- 1 - (ss_residual / ss_total)
r_squared
#######
# Fit a logistic regression model
logistic_model <- glm(sample_type ~ ENSG00000001629.10,
data = train_processed,
family = binomial)
# Print the logistic regression model
summary(logistic_model)
## Stepwise selection
## # To reduce the amount of genes we test, we will use the genes that were found in the decision tree model
keep_imp = t_test_results %>%
mutate(p_adj = p.adjust(p_value, method = "BH")) %>%
filter(p_adj < 0.00001) %>%
head(50) %>%
pull(Gene)
train_processed_cl = train_processed %>%
select(sample_type, all_of(keep_imp))
null_model <- glm(sample_type ~ 1,
data = train_processed_cl,
family = binomial) # Model with only the intercept
full_model <- glm(sample_type ~ .,
data = train_processed_cl,
family = binomial) # Model with all features
stepwise_model <- stats::step(null_model,
scope = list(lower = null_model,
upper = full_model),
direction = "both",
trace = 0)
summary(stepwise_model)
#######
library(tidymodels)
# 1. Define the Recipe (Pre-processing)
log_recipe <- recipe(sample_type ~ .,
data = train_processed_cl %>%
select(sample_type, 2:5))
# 2. Define the Model Specification
# Here we specify 'glm' as the engine for a classification task
logistic_spec <- logistic_reg() %>%
set_engine("glm") %>%
set_mode("classification")
# 3. Create a Workflow
# adding the stepwise selection is a bit tricky in tidymodels, as it does not have a built-in function for stepwise selection.
logistic_workflow <- workflow() %>%
add_recipe(log_recipe) %>% # Reusing the recipe with Variance Filtering
add_model(logistic_spec)
# 4. Fit the model
logistic_fit <- logistic_workflow %>%
fit(data = train_processed_cl)
#######
# Make predictions on the test dataset
predictions <- predict(logistic_model,
newdata = test_processed,
type = "response")
predicted_labels <- ifelse(predictions < 0.5, "Primary Tumor", "Solid Tissue Normal") %>%
as.factor()
# Evaluate the accuracy of the model
caret::confusionMatrix(predicted_labels, test_processed$sample_type)
#######
library(glmnet)
set.seed(123) # For reproducibility
# Prepare the data for glmnet
# Data must be in matrix format and labels must be numeric
train_matrix <- train_processed %>%
select(-sample_type) %>%
as.matrix()
train_labels <- as.numeric(train_processed$sample_type) - 1 # Convert
# Fit a Lasso regression model
# Before fitting a LASSO we have to find the best lambda (regularization parameter) using cross-validation
cv_lasso <- cv.glmnet(train_matrix,
train_labels,
alpha = 1, # Lasso
family = "binomial")
plot(cv_lasso)
best_lambda <- cv_lasso$lambda.min
best_lambda
lasso_model <- glmnet(train_matrix,
train_labels,
alpha = 1,
lambda = best_lambda,
family = "binomial")
# Print the Lasso regression model
# The coefficients of the model can be extracted using the coef function
coef(lasso_model) %>%
as.matrix() %>%
as.data.frame() %>%
tibble::rownames_to_column("Gene") %>%
filter(s0 != 0) %>%
arrange(desc(abs(s0))) %>%
filter(Gene != "(Intercept)") %>%
knitr::kable(caption = "Most Important Features in the Lasso Regression Model")
#######
library(tidymodels)
# 1. Define Lasso Specification
lasso_spec <- logistic_reg(
penalty = tune(), # We will find this via CV
mixture = 1 # 1 = Lasso
) %>%
set_engine("glmnet") %>%
set_mode("classification")
# 2. Define Cross-Validation Folds (10-fold)
set.seed(123)
folds <- vfold_cv(train_processed, v = 10)
# 3. Create a Grid of Lambda values to test
lambda_grid <- grid_regular(penalty(range = c(-5, 5)), levels = 10)
# 4. Run the Tuning
lasso_grid <- tune_grid(
workflow() %>%
add_recipe(full_recipe) %>%
add_model(lasso_spec),
resamples = folds,
grid = lambda_grid
)
# 5. Visualize the Tuning Results
# This replaces the plot(cv_lasso) call
autoplot(lasso_grid) +
theme_minimal() +
labs(title = "Lasso Tuning Results: Impact of Penalty on Accuracy")
#######
# 6. Select the best penalty (highest ROC AUC)
best_penalty <- lasso_grid %>%
select_best(metric = "roc_auc")
# 7. Finalize and Fit
final_lasso <- workflow() %>%
add_recipe(full_recipe) %>%
add_model(lasso_spec) %>%
finalize_workflow(best_penalty) %>%
fit(data = train_processed)
# 8. Extract Non-Zero Coefficients (Biomarkers)
final_lasso %>%
extract_fit_parsnip() %>%
tidy() %>%
filter(estimate != 0 & term != "(Intercept)") %>%
arrange(desc(abs(estimate))) %>%
knitr::kable(caption = "Biomarkers Selected by Lasso Regression")
#######
# Prepare the test data for prediction
# We need to ensure that the test data has the same features as the training data, and that the labels are in the correct format
test_matrix <- test_processed %>%
select(-sample_type) %>% # Exclude the label
as.matrix()
test_labels <- as.numeric(test_processed$sample_type) - 1 # Convert
# Make predictions on the test dataset
predictions <- predict(lasso_model,
newx = test_matrix,
type = "response")
predicted_labels <- ifelse(predictions > 0.5, 1, 0
) %>%
as.factor()
# Evaluate the accuracy of the model
confusionMatrix(predicted_labels, as.factor(test_labels))
## AUC and AUC-PR
# Get predicted probabilities for the positive class
pred_prob <- predict(lasso_model, newx = test_matrix, type = "response")
sample_type <- test_labels # Binary labels (1 for Tumor, 0 for Normal)
# Create a prediction object
pred <- ROCR::prediction(pred_prob, sample_type)
# Calculate AUC
auc <- ROCR::performance(pred,
measure = "auc")@y.values[[1]]
auc
# plot the ROC curve
roc_perf <- ROCR::performance(pred, measure = "tpr", x.measure = "fpr")
plot(roc_perf, col = "blue",
lwd = 2, main = "ROC Curve",
xlab = "False Positive Rate",
ylab = "True Positive Rate")
# Calculate AUC-PR
auc_pr <- ROCR::performance(pred, measure = "aucpr")@y.values[[1]]
auc_pr
# plot the Precision-Recall curve
pr_perf <- ROCR::performance(pred, measure = "prec", x.measure = "rec")
plot(pr_perf,
col = "red",
lwd = 2,
main = "Precision-Recall Curve",
xlab = "Recall",
ylab = "Precision")
#######
library(glmnet)
set.seed(123) # For reproducibility
# We use the same data as prepared for the Lasso regression, but we will set alpha = 0 for Ridge regression
cv_ridge <- cv.glmnet(train_matrix,
train_labels,
alpha = 0, # alpha = 0 for Ridge
family = "binomial")
plot(cv_ridge)
best_lambda_ridge <- cv_ridge$lambda.min
best_lambda_ridge
ridge_model <- glmnet(train_matrix,
train_labels,
alpha = 0, # alpha = 0 for Ridge
lambda = best_lambda_ridge,
family = "binomial")
# Print the Ridge regression model
coef(ridge_model) %>%
as.matrix() %>%
as.data.frame() %>%
tibble::rownames_to_column("Gene") %>%
arrange(desc(abs(s0))) %>%
filter(abs(s0) > 0.001) %>% # Filter for coefficients with a magnitude greater than a certain value
filter(Gene != "(Intercept)") %>%
knitr::kable(caption = "Most Important Features in the Ridge Regression Model")
#######
library(tidymodels)
set.seed(123)
# 1. Define Ridge Specification
ridge_spec <- logistic_reg(
penalty = tune(), # We will find this via CV
mixture = 0 # 0 = Ridge
) %>%
set_engine("glmnet") %>%
set_mode("classification")
# 2. Define Cross-Validation Folds (10-fold)
folds <- vfold_cv(train_processed, v = 10)
# 3. Create a Grid of Lambda values to test
lambda_grid <- grid_regular(penalty(range = c(-5, 0)), levels = 50)
# 4. Run the Tuning
ridge_grid <- tune_grid(
workflow() %>% add_recipe(full_recipe) %>% add_model(ridge_spec),
resamples = folds,
grid = lambda_grid
)
# 5. Visualize the Tuning Results
autoplot(ridge_grid) +
theme_minimal() +
labs(title = "Ridge Tuning Results: Impact of Penalty on Accuracy")
#######
# Prepare the test data for prediction
test_matrix <- as.matrix(test_processed %>%
select (-sample_type)) %>% # Exclude the label
as.matrix()
test_labels <- as.numeric(test_processed$sample_type) - 1 # Convert
# Make predictions on the test dataset
predictions <- stats::predict(ridge_model,
newx = test_matrix,
type = "response")
predicted_labels <- ifelse(predictions > 0.5, 1, 0)
predicted_labels <- as.factor(predicted_labels)
# Evaluate the accuracy of the model
# Evaluate the accuracy of the model
confusionMatrix(predicted_labels, as.factor(test_labels))
#######
## AUC and AUC-PR
# Get predicted probabilities for the positive class
pred_prob <- predict(lasso_model,
newx = test_matrix,
type = "response")
sample_type <- test_labels # Binary labels (1 for Tumor, 0 for Normal)
# Create a prediction object
pred <- ROCR::prediction(pred_prob, sample_type)
# Calculate AUC
auc <- ROCR::performance(pred, measure = "auc")@y.values[[1]]
auc
# plot the ROC curve
roc_perf <- ROCR::performance(pred, measure = "tpr", x.measure = "fpr")
plot(roc_perf, col = "blue",
lwd = 2, main = "ROC Curve",
xlab = "False Positive Rate",
ylab = "True Positive Rate")
# Calculate AUC-PR
auc_pr <- ROCR::performance(pred, measure = "aucpr")@y.values[[1]]
auc_pr
# plot the Precision-Recall curve
pr_perf <- ROCR::performance(pred, measure = "prec", x.measure = "rec")
plot(pr_perf,
col = "red",
lwd = 2,
main = "Precision-Recall Curve",
xlab = "Recall",
ylab = "Precision")