Gene-Expression Data Classification: Acute Lymphoblastic Leukemia vs. Acute Myeloid Leukemia
Kassie Shippey
1 Introduction: The Molecular Classification of Cancer: Gene Expression used for Class Predicition and Discovery
This data science project uses gene-expression data, harvested from a
proof-of-concept study published in 1999 by Golub et al. Using methods
such as DNA microarrays, one could identify new cancer classes, and
assign tumors to these classes. In this project, the data is used to
classify leukemia samples: acute lymphoblastic leukemia (ALL) or acute
myeloid leukemia (AML). I use two methods to analyze this data:
k-nearest neighbors (KNN) and Naive
Bayes. The training file contains patients 1–38 and the
independent file contains patients 39–72. The actual.csv
file supplies the known cancer type for all 72 patients. This
information, although part of the public domain, was directly sourced
from Kaggle.
This could be considered a high-dimensional problem; as there are 7,129 gene measurements but only 72 patients. Using every gene directly would make distance-based classification unstable. It would also make Naive Bayes unnecessarily noisy. Because this is true, I perform feature selection using the training data only and retain the 50 genes with the largest absolute two-sample t-statistics between ALL and AML. This is done to avoid using the independent test labels to choose features.
library(class)
set.seed(402)2 Read and organize the data
The gene-expression files are arranged accordingly: genes in rows,
patients in columns. Every patient-expression column is followed by a
qualitative call column. I have chosen to keep only the
numeric expression columns, transpose the data, and then have one row
per patient and one column per gene.
train_raw <- read.csv("data_set_ALL_AML_train.csv", check.names = FALSE)## Warning in scan(file = file, what = what, sep = sep, quote = quote, dec = dec,
## : EOF within quoted stringtest_raw <- read.csv("data_set_ALL_AML_independent.csv", check.names = FALSE)
actual <- read.csv("actual.csv")
# Identify the patient-expression columns.
train_patient_cols <- grep("^[0-9]+$", names(train_raw), value = TRUE)
test_patient_cols <- grep("^[0-9]+$", names(test_raw), value = TRUE)
# Gene names from each data set.
train_genes <- make.unique(as.character(train_raw[[2]]))
test_genes <- make.unique(as.character(test_raw[[2]]))
# Create expression matrices with genes as rows.
train_expr <- as.matrix(
train_raw[, train_patient_cols, drop = FALSE]
)
test_expr <- as.matrix(
test_raw[, test_patient_cols, drop = FALSE]
)
storage.mode(train_expr) <- "numeric"
storage.mode(test_expr) <- "numeric"
# Attach the correct gene names to each data set.
rownames(train_expr) <- train_genes
rownames(test_expr) <- test_genes
# Keep genes that occur in both data sets.
common_genes <- intersect(
rownames(train_expr),
rownames(test_expr)
)
train_expr <- train_expr[common_genes, , drop = FALSE]
test_expr <- test_expr[common_genes, , drop = FALSE]
# Transpose: patients are rows and genes are columns.
X_train <- t(train_expr)
X_test <- t(test_expr)
# Patient IDs.
train_ids <- as.integer(rownames(X_train))
test_ids <- as.integer(rownames(X_test))
# Match each patient to the known cancer diagnosis.
y_train <- factor(
actual$cancer[match(train_ids, actual$patient)],
levels = c("ALL", "AML")
)
y_test <- factor(
actual$cancer[match(test_ids, actual$patient)],
levels = c("ALL", "AML")
)
cat("Training samples:", nrow(X_train), "\n")## Training samples: 38cat("Independent test samples:", nrow(X_test), "\n")## Independent test samples: 34cat("Genes used:", ncol(X_train), "\n")## Genes used: 1647table(y_train)## y_train
## ALL AML
## 27 11table(y_test)## y_test
## ALL AML
## 20 14The training set has 38 patients and the independent test set has 34 patients. There are 7,129 measured genes. Choosing to keep the independent samples completely separate until final evaluation will give a more realistic estimate of how the methods perform on unseen patients.
3 Training-only feature selection
For each gene, I calculated a two-sample t-statistic comparing the ALL and AML patients in the training set. After this, I rank genes by the absolute value of this statistic and retain the top 50. A large absolute t-statistic means the gene has a relatively large difference between the two leukemia groups compared with within-group variation.
all_rows <- y_train == "ALL"
aml_rows <- y_train == "AML"
mean_all <- colMeans(X_train[all_rows, , drop = FALSE])
mean_aml <- colMeans(X_train[aml_rows, , drop = FALSE])
var_all <- apply(X_train[all_rows, , drop = FALSE], 2, var)
var_aml <- apply(X_train[aml_rows, , drop = FALSE], 2, var)
se_diff <- sqrt(var_all / sum(all_rows) + var_aml / sum(aml_rows))
t_score <- (mean_all - mean_aml) / se_diff
t_score[!is.finite(t_score)] <- 0
n_genes_keep <- 50
top_genes <- names(sort(abs(t_score), decreasing = TRUE))[1:n_genes_keep]
train_selected <- X_train[, top_genes, drop = FALSE]
test_selected <- X_test[, top_genes, drop = FALSE]
head(data.frame(Gene = top_genes,
Absolute_t = abs(t_score[top_genes])), 10)## Gene Absolute_t
## L13278_at L13278_at 6.281410
## L47738_at L47738_at 5.733623
## AF009426_at AF009426_at 5.693017
## D49950_at D49950_at 5.617638
## J05243_at J05243_at 5.453125
## HG1612-HT1612_at HG1612-HT1612_at 5.276869
## D63874_at D63874_at 5.187196
## D63880_at D63880_at 5.108403
## J03473_at J03473_at 5.044578
## D86967_at D86967_at 4.982889The test labels are not used in this step. I chose to do this because selecting genes using all 72 labels would leak information from the test set into model building and make the reported accuracy too optimistic and not as accurate as needed.
4 Standardization
KNN depends directly on distances. This means that genes on larger numeric scales could dominate the calculation, which is not what I am looking for. I standardize each selected gene using the training mean and training standard deviation, then apply those same training values to the independent test set. The standardized values are also used for Naive Bayes. This means that both methods receive the same features.
train_means <- colMeans(train_selected)
train_sds <- apply(train_selected, 2, sd)
train_sds[train_sds == 0] <- 1
train_z <- sweep(train_selected, 2, train_means, "-")
train_z <- sweep(train_z, 2, train_sds, "/")
test_z <- sweep(test_selected, 2, train_means, "-")
test_z <- sweep(test_z, 2, train_sds, "/")5 Method 1: k-nearest neighbors
KNN classifies a new patient by finding the k training
patients with the most similar standardized gene-expression profiles and
using their classes to predict the new patient’s class. A small
k can have low bias but high variance, while a larger
k smooths the decision boundary and can increase bias. I
use leave-one-out cross-validation (LOOCV) on the training set to
compare several odd values of k.
k_values <- c(1, 3, 5, 7, 9, 11)
cv_accuracy <- numeric(length(k_values))
for (j in seq_along(k_values)) {
k <- k_values[j]
cv_pred <- character(nrow(train_z))
for (i in seq_len(nrow(train_z))) {
cv_pred[i] <- as.character(
knn(train = train_z[-i, , drop = FALSE],
test = train_z[i, , drop = FALSE],
cl = y_train[-i],
k = k)
)
}
cv_accuracy[j] <- mean(cv_pred == as.character(y_train))
}
cv_results <- data.frame(k = k_values, Accuracy = cv_accuracy)
cv_results## k Accuracy
## 1 1 0.8947368
## 2 3 0.9736842
## 3 5 0.9736842
## 4 7 0.9473684
## 5 9 0.9473684
## 6 11 0.9473684plot(cv_results$k, cv_results$Accuracy,
type = "b", pch = 19,
xlab = "Number of neighbors (k)",
ylab = "LOOCV accuracy",
main = "KNN tuning on the training set")
# If several k values tie, choose the smallest value greater than 1;
# otherwise use the k with the highest CV accuracy.
best_acc <- max(cv_accuracy)
tied_k <- k_values[cv_accuracy == best_acc]
best_k <- if (any(tied_k > 1)) min(tied_k[tied_k > 1]) else min(tied_k)
best_k## [1] 3With this data, k=3 and k=5 are tied for the highest accuracy at about 97.4%. I use k = 3 because it is the smallest tied value above 1, giving some protection against a prediction being determined by a single unusual training observation.
knn_pred <- knn(train = train_z,
test = test_z,
cl = y_train,
k = best_k)
knn_conf <- table(Predicted = knn_pred, Actual = y_test)
knn_conf## Actual
## Predicted ALL AML
## ALL 19 5
## AML 1 9knn_accuracy <- mean(knn_pred == y_test)
knn_accuracy## [1] 0.8235294On the 34 independent patients, KNN correctly classifies 28 of 34 patients, for an accuracy of approximately 82.35%. The confusion matrix shows 19 correctly classified ALL, 9 correctly classified AML, and 6 total errors
6 Method 2: Gaussian Naive Bayes
Naive Bayes estimates the probability of each class using Bayes’ rule. I decided to use a Gaussian version. This means that models each selected standardized gene as normally distributed within each class. The “naive” assumption: genes are conditionally independent given the leukemia class. Gene-expression variables are not truly independent. Although this is true, Naive Bayes can still work well as a simple probabilistic classifier.
I implement the method directly. For each class and gene, the training data provide a mean and standard deviation. For a new patient: I add the log prior probability of the class, and then the Gaussian log-density for each selected gene. The class with the larger total log score is predicted.
fit_gaussian_nb <- function(x, y) {
classes <- levels(y)
priors <- sapply(classes, function(cl) mean(y == cl))
means <- sapply(classes, function(cl) colMeans(x[y == cl, , drop = FALSE]))
sds <- sapply(classes, function(cl) apply(x[y == cl, , drop = FALSE], 2, sd))
# Prevent a zero SD from causing undefined Gaussian densities.
sds[sds < 1e-8] <- 1e-8
list(classes = classes, priors = priors, means = means, sds = sds)
}
predict_gaussian_nb <- function(model, newx) {
scores <- matrix(NA_real_, nrow = nrow(newx), ncol = length(model$classes))
colnames(scores) <- model$classes
for (j in seq_along(model$classes)) {
log_density <- sapply(seq_len(ncol(newx)), function(g) {
dnorm(newx[, g],
mean = model$means[g, j],
sd = model$sds[g, j],
log = TRUE)
})
# rowSums combines evidence from all 50 genes.
scores[, j] <- log(model$priors[j]) + rowSums(log_density)
}
factor(model$classes[max.col(scores, ties.method = "first")],
levels = model$classes)
}nb_model <- fit_gaussian_nb(train_z, y_train)
nb_pred <- predict_gaussian_nb(nb_model, test_z)
nb_conf <- table(Predicted = nb_pred, Actual = y_test)
nb_conf## Actual
## Predicted ALL AML
## ALL 20 5
## AML 0 9nb_accuracy <- mean(nb_pred == y_test)
nb_accuracy## [1] 0.8529412On the independent test set, Gaussian Naive Bayes correctly classifies 29 of 34 patients, for an accuracy of approximately 85.29%. Its confusion matrix shows 20 correctly classified ALL, 9 correctly classified AML, and 5 errors.
7 Compare the two methods
Accuracy alone does not show which class of leukemia is being missed. Because of this, I also calculate class-specific sensitivity, specificity, and balanced accuracy, treating AML as the positive class; as it is often the most mis-diagnosed.
metrics <- function(actual, predicted, positive = "AML") {
actual <- as.character(actual)
predicted <- as.character(predicted)
negative <- setdiff(c("ALL", "AML"), positive)
TP <- sum(actual == positive & predicted == positive)
TN <- sum(actual == negative & predicted == negative)
FP <- sum(actual == negative & predicted == positive)
FN <- sum(actual == positive & predicted == negative)
sensitivity <- TP / (TP + FN)
specificity <- TN / (TN + FP)
c(Accuracy = (TP + TN) / length(actual),
AML_Sensitivity = sensitivity,
ALL_Specificity = specificity,
Balanced_Accuracy = (sensitivity + specificity) / 2)
}
comparison <- rbind(
KNN = metrics(y_test, knn_pred),
Naive_Bayes = metrics(y_test, nb_pred)
)
round(comparison, 3)## Accuracy AML_Sensitivity ALL_Specificity Balanced_Accuracy
## KNN 0.824 0.643 0.95 0.796
## Naive_Bayes 0.853 0.643 1.00 0.821barplot(comparison[, "Accuracy"],
ylim = c(0, 1),
ylab = "Independent test accuracy",
main = "Independent test accuracy by method")
Naive Bayes has slightly higher overall and balanced accuracy, while the two methods have the same AML sensitivity.
8 Examine individual predictions
prediction_table <- data.frame(
Patient = test_ids,
Actual = y_test,
KNN = knn_pred,
Naive_Bayes = nb_pred
)
prediction_table$KNN_Correct <- prediction_table$Actual == prediction_table$KNN
prediction_table$NB_Correct <- prediction_table$Actual == prediction_table$Naive_Bayes
prediction_table## Patient Actual KNN Naive_Bayes KNN_Correct NB_Correct
## 1 39 ALL ALL ALL TRUE TRUE
## 2 40 ALL ALL ALL TRUE TRUE
## 3 42 ALL ALL ALL TRUE TRUE
## 4 47 ALL ALL ALL TRUE TRUE
## 5 48 ALL ALL ALL TRUE TRUE
## 6 49 ALL ALL ALL TRUE TRUE
## 7 41 ALL ALL ALL TRUE TRUE
## 8 43 ALL ALL ALL TRUE TRUE
## 9 44 ALL ALL ALL TRUE TRUE
## 10 45 ALL ALL ALL TRUE TRUE
## 11 46 ALL ALL ALL TRUE TRUE
## 12 70 ALL ALL ALL TRUE TRUE
## 13 71 ALL ALL ALL TRUE TRUE
## 14 72 ALL ALL ALL TRUE TRUE
## 15 68 ALL ALL ALL TRUE TRUE
## 16 69 ALL ALL ALL TRUE TRUE
## 17 67 ALL AML ALL FALSE TRUE
## 18 55 ALL ALL ALL TRUE TRUE
## 19 56 ALL ALL ALL TRUE TRUE
## 20 59 ALL ALL ALL TRUE TRUE
## 21 52 AML ALL AML FALSE TRUE
## 22 53 AML AML ALL TRUE FALSE
## 23 51 AML AML AML TRUE TRUE
## 24 50 AML AML AML TRUE TRUE
## 25 54 AML ALL ALL FALSE FALSE
## 26 57 AML AML AML TRUE TRUE
## 27 58 AML ALL AML FALSE TRUE
## 28 60 AML AML AML TRUE TRUE
## 29 61 AML AML AML TRUE TRUE
## 30 65 AML AML AML TRUE TRUE
## 31 66 AML ALL ALL FALSE FALSE
## 32 63 AML AML ALL TRUE FALSE
## 33 64 AML ALL ALL FALSE FALSE
## 34 62 AML AML AML TRUE TRUEKNN misses 67, 52, 54, 58, 66, and 64, while Naive Bayes misses 53, 54, 66, 63, and 64. Patients 54, 64, and 66 are missed by both. More study and tuning will need to be done.
9 Discussion
Both KNN and Gaussian Naive Bayes performed decently well on this gene-expression dataset. Feature selection reduced the data to 50 genes. The standardization process ensured the models used comparable measurements. KNN achieved an independent-test accuracy of 82.4%. Naive Bayes achieved 85.3% independent-test accuracy. Both methods had the same AML sensitivity of 64.3%, although Naive Bayes correctly classified all ALL patients.
A major limitation that can be seen in this study is the small sample size of only 72 patients. The choice to retain 50 genes (out of 7000+) is also a modeling decision; rather than a biologically established cutoff of any kind. Future work could be used to test different numbers of genes, additional model settings, and larger independent datasets to ensure accuracy.
10 Conclusion
This project used two classification methods on the ALL/AML
gene-expression data provided by Golub and Kaggle. Using training-only
selection of the 50 strongest genes that could be found, KNN with
k = 3 correctly classified 28 of 34 independent patients,
while Naive Bayes correctly classified 29 of 34. The results and
findings show how concepts such as feature selection, standardization,
model tuning, and a genuinely held-out test set can be combined for a
“high-dimensional” classification problem.