Why this activity matters

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.

Learning goals

By the end of this activity, you should be able to:

  • describe why PCA is useful for large biological datasets
  • run PCA in R using prcomp()
  • interpret a scree plot
  • make PC1 vs PC2 and related PCA plots
  • color samples by cancer type
  • add a legend to a base R plot
  • think carefully about when a visualization is or is not informative

Load the NCI-60 data

[The NCI-60 dataset is included in the ISLR package. It contains gene expression data for cancer cell lines representing several cancer types.]

# Load the package containing the NCI-60 data.
# If this fails, ask an instructor to install the ISLR package.

library(ISLR)

# Cancer type labels for the cell lines
nci_labs <- NCI60$labs

# Gene expression data
# Rows = cancer cell lines
# Columns = genes
nci_data <- NCI60$data

First look at the data

Before running PCA, we should understand the objects we are working with.

# What kind of object is nci_labs?
class(nci_labs)
## [1] "character"
# How many labels are there?
length(nci_labs)
## [1] 64
# What cancer types are represented?
table(nci_labs)
## nci_labs
##      BREAST         CNS       COLON K562A-repro K562B-repro    LEUKEMIA 
##           7           5           7           1           1           6 
## MCF7A-repro MCF7D-repro    MELANOMA       NSCLC     OVARIAN    PROSTATE 
##           1           1           8           9           6           2 
##       RENAL     UNKNOWN 
##           9           1
# What kind of object is nci_data?
class(nci_data)
## [1] "matrix" "array"
# What are its dimensions?
dim(nci_data)
## [1]   64 6830
# Look at a small piece of the matrix
nci_data[1:5, 1:5]
##           1         2         3         4         5
## V1 0.300000  1.180000  0.550000  1.140000 -0.265000
## V2 0.679961  1.289961  0.169961  0.379961  0.464961
## V3 0.940000 -0.040000 -0.170000 -0.040000 -0.605000
## V4 0.280000 -0.310000  0.680000 -0.810000  0.625000
## V5 0.485000 -0.465000  0.395000  0.905000  0.200000

Question: What do the rows represent? What do the columns represent?

Your answer: The rows are individual cancer cell lines. The columns are gene expression levels.

Why PCA is needed

This dataset has many more genes than samples.

n_cell_lines <- nrow(nci_data)
n_genes <- ncol(nci_data)

paste("There are", n_cell_lines, "cell lines and", n_genes, "genes.")
## [1] "There are 64 cell lines and 6830 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.

Run PCA

We use the base R function prcomp().

The argument scale. = TRUE scales each gene to have standard deviation 1 before PCA is performed.

pca_out <- prcomp(nci_data, scale. = TRUE)

The PCA results contain several components.

names(pca_out)
## [1] "sdev"     "rotation" "center"   "scale"    "x"

The most important part for plotting samples is:

head(pca_out$x[, 1:5])
##          PC1       PC2         PC3         PC4        PC5
## V1 -19.68245  3.527748  -9.7354382   0.8177816 -12.511081
## V2 -22.90812  6.390938 -13.3725378  -5.5911088  -7.972471
## V3 -27.24077  2.445809  -3.5053437   1.3311502 -12.466296
## V4 -42.48098 -9.691742  -0.8830921  -3.4180227 -41.938370
## V5 -54.98387 -5.158121 -20.9291076 -15.7253986 -10.361364
## V6 -26.96488  6.727122 -21.6422924 -13.7323153   7.934827

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.

pve <- 100 * pca_out$sdev^2 / sum(pca_out$sdev^2)

head(pve)
## [1] 11.358942  6.756203  5.751842  4.247554  3.734972  3.618630

Scree plot and cumulative variance plot

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.

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: Around 5–10 PCs offer a streamlined yet faithful representation of the dataset, since variance gains taper off beyond that range.

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.

# 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
##      BREAST         CNS       COLON K562A-repro K562B-repro    LEUKEMIA 
##   "#FF0000"   "#FF6D00"   "#FFDB00"   "#B6FF00"   "#49FF00"   "#00FF24" 
## MCF7A-repro MCF7D-repro    MELANOMA       NSCLC     OVARIAN    PROSTATE 
##   "#00FF92"   "#00FFFF"   "#0092FF"   "#0024FF"   "#4900FF"   "#B600FF" 
##       RENAL     UNKNOWN 
##   "#FF00DB"   "#FF006D"

PCA scatter plots

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.

Helper function for PCA plots

We can write our own functions in R! This function makes PCA plots and automatically adds:

  • axis labels with percent variance explained
  • cancer-type colors
  • a legend
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

plot_pca_pair(1, 2, legend_position = "topright")

Question: Do any cancer types cluster together clearly?

Your answer: Yes, melanoma, renal, and leukemia samples form tight, visually distinct clusters, while breast and colon show moderate grouping and NSCLC appears more scattered across the space.

PC1 vs PC3

Sometimes PC1 vs PC2 is not the most informative view.

plot_pca_pair(1, 3, legend_position = "topright")

Question: Does this view separate any cancer types better than PC1 vs PC2?

Your answer: Yes, PC1 vs PC3 gives a cleaner separation for a few groups. Leukemia stands out more distinctly along PC3, and melanoma also shows a clearer vertical spread compared with the PC1 vs PC2 view. Other types such as NSCLC and ovarian remain fairly diffuse in both projections.

PC2 vs PC3

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 PC3 is the most informative view. It gives clearer separation for leukemia and melanoma, and it highlights structure that is less visible in the PC1 vs PC2 and PC2 vs PC3 projections.

A cleaner legend outside the plot

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.

# Save old graphics settings
old_par <- par(no.readonly = TRUE)

# Make extra room on the right side
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"
)

# Restore old graphics settings
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.

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 labels add useful context because they reveal which clusters correspond to specific cancer types, but the plot becomes somewhat crowded as overlapping points and similar colors make fine distinctions harder to see.

Why the PCA plot may not look dramatic

A PCA plot does not always produce clean, separated clusters.

That does not mean PCA failed.

Possible reasons include:

  • cancer biology is complex
  • some cancer types are more similar than others
  • PC1 and PC2 may capture technical or broad expression differences
  • important subtype differences may appear in later PCs
  • gene expression varies within each cancer type

A good scientific visualization does not have to show perfect separation. It should help us ask better questions.

Challenge 1: Try other PCs

Modify this code to try other PC pairs.

plot_pca_pair(3, 4, legend_position = "topright")

Which pair gives the clearest separation?

Your answer: PC1 vs PC3 gives the clearest separation because it highlights distinct spacing for leukemia and melanoma, while the other pairs show more overlap among major groups.

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.

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. Showing only a few cancer types makes the structure much easier to see because the clusters no longer compete for space and the separation along the principal components becomes more visually distinct.

Final reflection

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:

  • cancer type
  • gene expression heatmaps
  • clinical variables
  • tumor subtype labels

This is one way computational biology connects large molecular datasets to biological interpretation.