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: Rows represent cancer samples and columns represent genes.

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: As it can be observed on the plots pc7 pc8 is where the line starts to flatten, the “elbow” point. PCs after that point do not have as much variance and can be discarded to simplify the data.

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: Melanoma seems to be clustered on the bottom mid-left side. Colon cancer samples are also clustered in the top mid-right section. CNS and renal cancer types also seem not to be too far from each other.

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: Each cancer type is much seperated than before. The most noticable changes in proximity can be observed with Leukemia, CNS and ovarian cancer types.

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 PC2 captures the highest percentage of variance. It clearly seperates melanoma, leukemia and colon cancer from other cancer types, while also showing clear clustering. Therefore, it is the most informative PCA plot.

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: On the top, where there are already many points, adding labels makes it more crowded. It was already easy to differentiate cancer types due to coloring, so in this particular case, labels are counterproductive, in my opinion.

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(2, 4, legend_position = "topright")

Which pair gives the clearest separation?

Your answer: Out of all the combinations includinf PC4, the one with the clearest seperation is PC1 vs PC4, since PC1 has the highest percentage of variance and gives more room for seperation.

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, it is much easier. With only 3 cancer types to pay attention too, we get much clearer picture without all the additional noise. It can be useful when we need information only about these 3 types.

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.