Alon colon-cancer data
Section of High-dimensional analysis, part of Statistical learning
Alon colon-cancer data
| Type | Section |
| Taught in | High-dimensional analysis |
| Part | Statistical learning |
The calculation below uses expression measurements for 2000 genes from 62 colon tissue samples, 40 tumour and 22 normal. Centring and scaling a feature to mean zero and variance one does not make it normally distributed. The code expects an expression matrix x of size 62\times 2000 and labels y. The printed principal-component, testing, and lasso numbers are the results of that analysis.
load("alon.rda") # x: 62 x 2000 expression matrix; y: tumour/normal labels
library(FactoMineR)
numpcs <- 10
res <- PCA(x, scale.unit = FALSE, graph = FALSE, ncp = numpcs)
plot(res$ind$coord[, 1], res$ind$coord[, 2], col = y,
xlab = "PC1", ylab = "PC2")
legend("topleft", legend = c("Normal", "Tumour"),
col = c("black", "red"), pch = c(1, 1))
# Welch t-test of each principal component between the classes
ttestres <- numeric(numpcs)
for (i in 1:10) {
ttestres[i] <- t.test(res$ind$coord[, i] ~ y)$p.value
print(paste(i, signif(ttestres[i], digits = 3)))
}
plot(res$ind$coord[, 3], res$ind$coord[, 6], col = y,
xlab = "PC3", ylab = "PC6")
# top 20 gene contributors, by absolute loading, on PC3 and PC6
lst3 <- sort(abs(res$var$coord[, 3]), index.return = TRUE, decreasing = TRUE)
lst6 <- sort(abs(res$var$coord[, 6]), index.return = TRUE, decreasing = TRUE)
lst3ids <- lapply(lst3, `[`, lst3$x %in% head(unique(lst3$x), 20))$ix
lst6ids <- lapply(lst6, `[`, lst6$x %in% head(unique(lst6$x), 20))$ix
intersect(lst3ids, lst6ids) # none in the reported analysis
# one Welch t-test per gene (40 tumour, 22 normal)
orig_p <- rep(0, 2000)
for (i in 1:2000) {
orig_p[i] <- t.test(x[1:40, i], x[41:62, i])$p.value
}
sorted <- sort(orig_p, index.return = TRUE)
sorted_p <- sorted$x
sorted_p_ix <- sorted$ix
# Benjamini-Hochberg at nominal level 0.01
thres <- 0.01
kk <- rep(0, 2000)
for (i in 1:2000) {
if (sorted_p[i] <= thres * i / 2000) kk[i] <- i
}
cut_off_k <- max(kk) # zero if no inequality succeeds
significant_genes <- sorted_p_ix[seq_len(cut_off_k)]
order_genes_relabel <- 1:2000
significant <- c(rep("red", cut_off_k), rep("blue", 2000 - cut_off_k))
plot(order_genes_relabel, sorted_p, log = "xy",
xlab = "gene number ordered by p-value",
ylab = "unadjusted p-value",
main = "Genes against p-values, BH threshold at FDR 0.01",
col = significant)
abline(v = cut_off_k, lty = 2)
abline(a = 0, b = thres / 2000, untf = TRUE)
intersect(significant_genes, lst3ids)
intersect(significant_genes, lst6ids)
library(glmnet)
cv.lasso <- cv.glmnet(x, y, alpha = 1, family = "binomial")
lambda <- cv.lasso$lambda.min
lasso.model <- glmnet(x, y, alpha = 1, family = "binomial", lambda = lambda)
beta_vals <- as.matrix(lasso.model$beta)
selected_genes <- which(beta_vals != 0)
selected_betas <- beta_vals[selected_genes]
beta_0 <- lasso.model$a0
intersect(significant_genes, selected_genes)Principal components of the matrix
Main article: Principal components of the colon-cancer matrix
The first ten principal components were examined, for three reasons. One or more components might pick up variation associated with the class difference. One can look for outliers, including possibly mislabelled observations. The dimension reduction allows a full suite of classifiers to be constructed. A plot of the first component against the second offered almost no separation between the classes and indicated no outliers. Welch t-tests of each component between the classes gave
| Component | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 |
|---|---|---|---|---|---|---|---|---|---|---|
| p-value | 0.479 | 0.339 | 1.8\times 10^{-7} | 0.0995 | 0.0958 | 0.0166 | 0.346 | 0.84 | 0.915 | 0.431 |
The third component appears very aligned with the class differences, and the sixth is also promising. Their scatter plot separates the classes substantially, with some overlap. A simple classifier such as quadratic discrimination in these two dimensions might separate them fairly accurately. The labels were used in choosing these dimensions, so the analysis would need to be redone with left-out data to estimate classification accuracy honestly. The genes with the largest magnitude in the third and sixth direction vectors were listed, the top 20 of each. There were no genes in common between these two lists.
False discoveries among the genes
Main article: False discoveries among the genes
A Welch t-test was run for each gene, comparing the tumour and normal groups, using Benjamini–Hochberg at nominal level 0.01. Under the required null calibration and dependence assumptions, this targets an expected false-discovery proportion at most 1\%. It does not cap the false fraction in this realised gene list. The procedure declared 173 genes significant. Twelve of the 20 top third-component genes were among them,
1635,\;493,\;625,\;765,\;1511,\;1634,\;1153,\;992,\;1293,\;1002,\;1959,\;609,
and one of the top sixth-component genes, gene 286.
The logistic lasso
Main article: Logistic lasso
The lasso was used to perform model selection at the same time as fitting a logistic regression of the class on all the genes. The binomial deviance was the cost, with \lambda chosen as the value with the lowest binomial deviance on left-out data by cross-validation. This selected 14 gene indices
14,\;249,\;377,\;493,\;576,\;625,\;792,\;1360,\;1473,\;1582,\;1679,\;1772,\;1843,\;1924,
with fitted coefficients, in the same order,
\begin{aligned} &-0.111,\;-0.081,\;-0.588,\;-0.275,\;0.136,\;0.118,\;-0.116,\\ &0.089,\;0.113,\;0.085,\;0.141,\;0.302,\;-0.207,\;-0.141. \end{aligned}
The largest in magnitude are gene 377 (-0.588) and gene 1772 (+0.302), and the smallest is gene 249 (-0.081). Eleven of these were on the false-discovery list (14, 249, 377, 493, 625, 1473, 1582, 1679, 1772, 1843, 1924) and two were on the third-component list (493, 625).
The reported list contains 173 genes. A laboratory able to follow up only ten might prioritise the overlap with the lasso list and then apply an additional ranking. One such prioritisation retains genes 14, 377, 493, 625, 1473, 1582, 1679, 1772, 1843, 1924 after removing gene 249 from the eleven-gene overlap. This is a prioritisation rule, not a new false-discovery theorem: an arbitrary data-dependent subset of the discoveries need not inherit the original bound. Replication, measurement quality, and biological follow-up can test the candidates. The current data do not establish a causal role for them.
A principal-component analysis of the matrix followed by a lasso on the component scores selects combinations of all the genes, since each component is already a combination of all 2000 genes. That does not answer which genes. A lasso on the genes themselves is the fit that matches the question.
Used in
- Taught in High-dimensional analysis
- Also used in Graphical models
- All 9 questions on High-dimensional analysis
See also
Rank bound · Family-wise control implies false-discovery control · Benjamini–Hochberg · Benjamini–Hochberg controls the false discovery rate · Why classical methods fail when p > n · Rank of the sample covariance · Regression when p > n · Where high-dimensional columns come from