Cancer datasets often contain thousands of measurements for each sample. For example, a single cell line or patient tumor may have expression measurements for thousands of genes.
That creates an important problem:
How can we see patterns in data with thousands of dimensions?
Principal Component Analysis, or PCA, is one way to reduce a large dataset to a smaller number of summary variables while keeping as much of the original information as possible.
In this activity, we will use PCA to explore the NCI-60 cancer cell line dataset.
By the end of this activity, you should be able to:
prcomp()The NCI-60 dataset is included in the ISLR package. It
contains gene expression data for cancer cell lines representing several
cancer types.
```{r load-data, message=FALSE, warning=FALSE} # Load the package containing the NCI-60 data. # If this fails, ask an instructor to install the ISLR package.
library(ISLR)
nci_labs <- NCI60$labs
nci_data <- NCI60$data
---
## First look at the data
Before running PCA, we should understand the objects we are working with.
```{r inspect-data}
# What kind of object is nci_labs?
class(nci_labs)
# How many labels are there?
length(nci_labs)
# What cancer types are represented?
table(nci_labs)
```{r inspect-expression-data} # What kind of object is nci_data? class(nci_data)
dim(nci_data)
nci_data[1:5, 1:5]
**Question:** What do the rows represent? What do the columns represent?
> Your answer: Rows represent the individual cancer cell lines (samples). Columns represent the gene expression measurements for thousands of genes (features).
---
## Why PCA is needed
This dataset has many more genes than samples.
```{r sample-gene-counts}
n_cell_lines <- nrow(nci_data)
n_genes <- ncol(nci_data)
paste("There are", n_cell_lines, "cell lines and", n_genes, "genes.")
It would be impossible to directly visualize all genes at once. PCA gives us a smaller number of new variables called principal components.
Each principal component is a weighted combination of the original gene expression measurements.
We use the base R function prcomp().
The argument scale. = TRUE scales each gene to have
standard deviation 1 before PCA is performed.
```{r run-pca} pca_out <- prcomp(nci_data, scale. = TRUE)
The PCA results contain several components.
```{r pca-structure}
names(pca_out)
The most important part for plotting samples is:
```{r pca-scores-preview} head(pca_out$x[, 1:5])
The object `pca_out$x` contains the coordinates of each cell line in PCA space.
---
## How much variation does each PC explain?
We calculate the percent of variance explained by each principal component.
```{r pve}
pve <- 100 * pca_out$sdev^2 / sum(pca_out$sdev^2)
head(pve)
The scree plot shows how much variance each PC explains.
The cumulative plot shows how much total variance is explained as we include more PCs.
```{r pve-plots, fig.width=10, fig.height=4} par(mfrow = c(1, 2))
plot( pve, type = “b”, pch = 19, xlab = “Principal component”, ylab = “Percent variance explained”, main = “Scree plot” )
plot( cumsum(pve), type = “b”, pch = 19, xlab = “Number of principal components”, ylab = “Cumulative percent variance explained”, main = “Cumulative variance” )
par(mfrow = c(1, 1))
**Question:** About how many PCs would you keep if you wanted a simplified view of the data?
> Your answer: Typically around 5 to 7 PCs. You look for the "elbow" in the scree plot where the curve flattens out, indicating later PCs mostly capture background noise.
---
## Prepare colors for cancer types
We want each cancer type to have its own color.
This is useful because PCA itself does not use the cancer labels. PCA is unsupervised. We add the labels afterward to see whether the PCA structure lines up with known biology.
```{r colors}
# Get the unique cancer types
cancer_types <- sort(unique(nci_labs))
# Choose colors
cancer_colors <- rainbow(length(cancer_types))
# Name the colors by cancer type
names(cancer_colors) <- cancer_types
# Assign one color to each cell line
point_colors <- cancer_colors[nci_labs]
# Check the mapping
cancer_colors
Plotting PC1 (Principal Component 1) versus PC2 (Principal Component 2) is the standard method for visualizing high-dimensional data. It allows you to see the maximum amount of information (variance) in a dataset on a single 2D graph.
We can write our own functions in R! This function makes PCA plots and automatically adds:
```{r pca-plot-function} plot_pca_pair <- function(pc_x, pc_y, legend_position = “topright”) {
x_values <- pca_out\(x[, pc_x] y_values <- pca_out\)x[, pc_y]
x_label <- paste0(“PC”, pc_x, ” (“, round(pve[pc_x], 1),”%)“) y_label <- paste0(”PC”, pc_y, ” (“, round(pve[pc_y], 1),”%)“)
plot( x_values, y_values, col = point_colors, pch = 19, xlab = x_label, ylab = y_label, main = paste(“NCI-60 PCA: PC”, pc_x, “vs PC”, pc_y) )
legend( legend_position, legend = cancer_types, col = cancer_colors, pch = 19, cex = 0.75, bty = “n” ) }
---
### PC1 vs PC2
```{r pc1-pc2, fig.width=8, fig.height=6}
plot_pca_pair(1, 2, legend_position = "topright")
Question: Do any cancer types cluster together clearly?
Your answer: Yes, Leukemia and Melanoma typically form tight, distinct clusters. Others, like Breast cancer, are more scattered, indicating higher biological variation within those types.
Sometimes PC1 vs PC2 is not the most informative view.
```{r pc1-pc3, fig.width=8, fig.height=6} plot_pca_pair(1, 3, legend_position = “topright”)
**Question:** Does this view separate any cancer types better than PC1 vs PC2?
> Your answer: Yes, different principal components capture different sources of variance. Clusters that overlap in PC2 might separate cleanly along PC3.
---
### PC2 vs PC3
```{r pc2-pc3, fig.width=8, fig.height=6}
plot_pca_pair(2, 3, legend_position = "topright")
Question: Which PCA plot is most informative: PC1 vs PC2, PC1 vs PC3, or PC2 vs PC3?
Your answer: PC1 vs PC2 explains the most variance mathematically and is usually the best starting point, but the most “informative” plot depends on the specific biological variation or subtype differences you are trying to study.
Sometimes the legend covers points. In base R, we can make extra space on the right side and place the legend outside the plotting region.
```{r legend-outside, fig.width=9, fig.height=6} # Save old graphics settings old_par <- par(no.readonly = TRUE)
par(mar = c(5, 5, 4, 10), xpd = TRUE)
plot( pca_out\(x[, 1], pca_out\)x[, 2], col = point_colors, pch = 19, xlab = paste0(“PC1 (”, round(pve[1], 1), “%)”), ylab = paste0(“PC2 (”, round(pve[2], 1), “%)”), main = “NCI-60 PCA with legend outside plot” )
legend( “topright”, inset = c(-0.35, 0), legend = cancer_types, col = cancer_colors, pch = 19, cex = 0.8, bty = “n” )
par(old_par)
This is often the best base R option when there are many groups.
---
## Add cell line labels
Colors show cancer type, but we may also want to know which point is which cell line.
This plot adds labels. It can become crowded, so use it only when helpful.
```{r text-labels, fig.width=9, fig.height=7}
plot(
pca_out$x[, 1],
pca_out$x[, 2],
col = point_colors,
pch = 19,
xlab = paste0("PC1 (", round(pve[1], 1), "%)"),
ylab = paste0("PC2 (", round(pve[2], 1), "%)"),
main = "NCI-60 PCA with cancer-type labels"
)
text(
pca_out$x[, 1],
pca_out$x[, 2],
labels = nci_labs,
pos = 3,
cex = 0.55
)
Question: Are the labels helpful, or does the plot become too crowded?
Your answer: The plot becomes far too crowded (overplotting). Text labels are only helpful if you filter the data down to just a few samples first to identify specific outliers.
A PCA plot does not always produce clean, separated clusters.
That does not mean PCA failed.
Possible reasons include:
A good scientific visualization does not have to show perfect separation. It should help us ask better questions.
Modify this code to try other PC pairs.
```{r challenge-pc-pairs} plot_pca_pair(3, 4, legend_position = “topright”)
Which pair gives the clearest separation?
> Your answer: PC1 vs PC2 generally provides the best overall separation of major groups, though looking at deeper PCs can occasionally separate closely related subtypes that are lumped together earlier.
---
## Challenge 2: Focus on a few cancer types
Sometimes plots are easier to interpret if we subset the data.
Here we keep only three cancer types.
```{r subset-cancer-types}
keep_types <- c("BREAST", "LEUKEMIA", "MELANOMA")
keep <- nci_labs %in% keep_types
plot(
pca_out$x[keep, 1],
pca_out$x[keep, 2],
col = point_colors[keep],
pch = 19,
xlab = paste0("PC1 (", round(pve[1], 1), "%)"),
ylab = paste0("PC2 (", round(pve[2], 1), "%)"),
main = "Subset of cancer types"
)
legend(
"topright",
legend = keep_types,
col = cancer_colors[keep_types],
pch = 19,
bty = "n"
)
Question: Is it easier to see structure when fewer cancer types are shown?
Your answer: Yes, drastically. Removing the visual noise from overlapping clusters makes the relationships between the selected cancer types much easier to see.
In this activity, you used PCA to reduce thousands of gene expression measurements to a few principal components.
You also learned an important scientific lesson:
A visualization may be useful even if it does not produce perfect clusters.
PCA helps us explore structure, generate hypotheses, and decide what questions to ask next.
Later, we can compare PCA patterns with:
This is one way computational biology connects large molecular datasets to biological interpretation.