Principal Component Analysis (PCA)

Author

Daniel Sun, Amanda Ng

Introduction

Welcome to the Principal Component Analysis (PCA) module. By the end of this workshop, you will be able to:

  • Explain what PCA does and why it is useful in biological research
  • Recognize what PCA output looks like (scores, loadings, a score plot)
  • Apply three standard criteria to decide how many principal components to retain:
    • Kaiser criterion (eigenvalue > 1)
    • Variance explained threshold (>=80-90%)
    • Scree plot visual inspection
  • Interpret loadings to understand what each principal component represents biologically
  • Carry out a basic PCA in both R and Python

Remember: if you get stuck on a coding exercise, click the Hint button for a nudge, or the Solution button if you need the full answer.

Motivation

Why do we need PCA?

In human biology, datasets are almost always high-dimensional. You might measure dozens of clinical traits, hundreds of gene expression values, or thousands of metabolite concentrations on each individual. Analyzing every variable separately is slow, and many of those variables are correlated with each other.

PCA is a technique that takes a set of correlated variables and replaces them with a smaller set of uncorrelated summary variables called principal components (PCs). Each PC is a weighted combination of the original variables, constructed so that:

  1. PC1 captures as much variation in the data as possible.
  2. PC2 captures as much of the remaining variation as possible (and is uncorrelated with PC1).
  3. And so on…
Note

PCA does not require you to have a response variable. It is an unsupervised method used to explore and summarize your data.

A visual intuition

The best way to understand PCA is to see what is happening. The interactive 3-D plot and static 2-D plot below show the same simulated dataset of 200 people (height, weight, BMI) before and after PCA. Drag to rotate the 3D plot on the left to see how the cloud of points is oriented in three dimensions. The plot on the right shows the same data after PCA collapses all three variables into two axes, together preserving 99.9% of the variance.

The 3D plot shows all three variables simultaneously as a cloud of points. No single viewing angle captures all the structure at once. The PCA plot on the right finds the best possible 2D view automatically. Note that being able to reduce to just 2 components and still make a useful scatterplot is not always possible it depends on whether most variance is concentrated in the first two PCs.

NoteThe challenge of high-dimensional data

While 3D data is manageable, the real issue arises with high-dimensional biological datasets. If you collect 4, 10, or 10,000 variables (like gene expression networks or extensive metabolic panels) on each patient, it becomes completely impossible to visualize the data geometrically. You cannot plot a 4D or 50D scatterplot. PCA is powerful because it compresses these hyper-dimensional spaces down to fewer principal components, allowing us to visualize, cluster, and analyze complex biological structures that would otherwise remain hidden.

Section 1: The Data

Study background

This module uses data from by a published RNA-seq study of human folliculogenesis (GSE107746, Zhang et al. 2018). That study profiled oocytes and granulosa cells from 148 samples across five follicular development stages using RNA-seq, generating log2 FPKM expression values for thousands of genes. PCA was a central tool in that paper for showing how samples from different stages separate in gene expression space.

For this module we use four marker genes that show strong stage-specific expression in folliculogenesis, based on the biology reported in the paper:

Gene Role
NOBOX Oocyte-specific transcription factor; high in early stages
GDF9 Oocyte growth factor; rises through primary and secondary stages
FSHR FSH receptor in granulosa cells; increases toward antral stage
LHCGR LH/CG receptor; peaks at preovulatory stage

The 148 samples span five stages: Primordial (n=25), Primary (n=40), Secondary (n=18), Antral (n=47), and Preovulatory (n=18).

Note

The expression values below are simulated to match the group structure and approximate expression patterns described in the paper. They are not the raw GEO data, but they capture the biological signal that PCA was used to reveal.

1.1 Explore the data

Run the code below to load and inspect the dataset. You do not need to change anything. Just click Run Code and examine the output.

1.2 Visualize pairwise relationships

Before running PCA, it helps to see how correlated the variables are. Highly correlated variables are a sign that PCA may be able to compress the data effectively.

Section 2: Running PCA

2.1 Why scale first?

PCA is sensitive to the scale of variables. A variable with a large range will dominate one with a small range, even if they are conceptually equally important. The standard fix is to standardize each variable to have mean 0 and standard deviation 1 before running PCA (scale. = TRUE in R; StandardScaler in Python).

View code for plot 
raw    <- as.matrix(follicle_data[, 1:4])
scaled <- scale(raw)
par(mfrow = c(1, 2), mar = c(6, 4, 3, 1))
boxplot(raw,    main = "Before scaling",          ylab = "log2 FPKM",         col = "steelblue",   las = 2)
boxplot(scaled, main = "After scaling (mean=0, sd=1)", ylab = "Standardized value", col = "forestgreen", las = 2)
par(mfrow = c(1, 1))
View code for plot 
from sklearn.preprocessing import StandardScaler
raw    = follicle_data[["NOBOX","GDF9","FSHR","LHCGR"]].values
scaled = StandardScaler().fit_transform(raw)
fig, axes = plt.subplots(1, 2, figsize=(9, 4))
axes[0].boxplot(raw,    labels=["NOBOX","GDF9","FSHR","LHCGR"]); axes[0].set_title("Before scaling"); axes[0].set_ylabel("log2 FPKM")
axes[1].boxplot(scaled, labels=["NOBOX","GDF9","FSHR","LHCGR"]); axes[1].set_title("After scaling (mean=0, sd=1)"); axes[1].set_ylabel("Standardized value")
plt.tight_layout(); plt.show()
Tip

As a rule of thumb: always scale your variables unless you have a strong reason to believe they are already on the same scale.

2.2 Perform PCA: fill in the blanks

Complete the code to run PCA on all four gene expression columns, with scaling. The pca_result object you create here will carry forward into Sections 3, 4, and 5.

NoteHint

The four gene columns are columns 1:4. Use scale. = TRUE to standardize before PCA.

TipSolution
pca_result <- prcomp(follicle_data[, 1:4], scale. = TRUE) summary(pca_result) 
pca_result <- prcomp(follicle_data[, 1:4], scale. = TRUE)
summary(pca_result)
NoteHint

Use scaler.fit_transform(X) to scale. Then PCA(n_components=4) retains all 4 PCs.

TipSolution 
scaler   = StandardScaler()
X_scaled = scaler.fit_transform(X)
pca      = PCA(n_components=4)
pca.fit(X_scaled)
print("Variance explained by each PC:")
print(np.round(pca.explained_variance_ratio_, 4))
print("\nCumulative variance explained:")
print(np.round(np.cumsum(pca.explained_variance_ratio_), 4))

Sections 3-5: How Many Components to Keep?

Running PCA gives you as many PCs as input variables (4 here). The goal is to reduce dimensionality, so we need principled rules for deciding how many PCs are sufficient.

There are three widely used approaches, each covered in its own section:

Method Rule
Section 3 — Kaiser criterion Keep PCs with eigenvalue > 1
Section 4 — Variance threshold Keep PCs until cumulative variance >= 80-90%
Section 5 — Scree plot Keep PCs before the “elbow” in the eigenvalue plot

Section 3: The Kaiser Criterion

The eigenvalue of a PC equals the variance it explains (in standardized units, where each original variable has variance = 1). A PC with eigenvalue > 1 captures more information than any single original variable, so it is worth retaining.

Note

In R: eigenvalues = pca_result$sdev^2
In Python: eigenvalues = pca.explained_variance_

Exercise: Extract eigenvalues and apply the Kaiser criterion

NoteHint

pca_result$sdev gives standard deviations. Square them with ^2 to get eigenvalues. The Kaiser threshold is > 1.

TipSolution
eigenvalues <- pca_result$sdev ^ 2 print(round(eigenvalues, 3)) kaiser_keep <- eigenvalues > 1 print(paste("Keep PCs:", which(kaiser_keep))) 
eigenvalues <- pca_result$sdev ^ 2
print(round(eigenvalues, 3))
kaiser_keep <- eigenvalues > 1
print(paste("Keep PCs:", which(kaiser_keep)))
NoteHint

Use pca.explained_variance_ (not explained_variance_ratio_). The threshold is > 1.

TipSolution 
eigenvalues = pca.explained_variance_
print("Eigenvalues:", np.round(eigenvalues, 3))
kaiser_keep = eigenvalues > 1
print("PCs to keep (0-based index):", np.where(kaiser_keep)[0])
print("Number of PCs to keep:", np.sum(kaiser_keep))
WarningKaiser caveat

The Kaiser criterion is a rule of thumb, not a law. Always use it alongside the variance threshold and scree plot before deciding how many PCs to retain.

Section 4: Variance Explained Threshold

Keep enough PCs to explain at least 80-90% of the total variance. The right threshold depends on your field and how much information loss you can tolerate.

Exercise: Compute cumulative variance explained

NoteHint

cumsum() builds a running cumulative total. The threshold values are 0.80 and 0.90.

TipSolution
prop_var <- summary(pca_result)$importance["Proportion of Variance", ] cum_var <- cumsum(prop_var) print(round(prop_var, 3)) print(round(cum_var, 3)) n_80 <- which(cum_var >= 0.80)[1] cat("PCs needed for >= 80% variance:", n_80, "\n") n_90 <- which(cum_var >= 0.90)[1] cat("PCs needed for >= 90% variance:", n_90, "\n") 
prop_var <- summary(pca_result)$importance["Proportion of Variance", ]
cum_var  <- cumsum(prop_var)
print(round(prop_var, 3))
print(round(cum_var, 3))
n_80 <- which(cum_var >= 0.80)[1]
cat("PCs needed for >= 80% variance:", n_80, "\n")
n_90 <- which(cum_var >= 0.90)[1]
cat("PCs needed for >= 90% variance:", n_90, "\n")
NoteHint

np.cumsum() gives a running total. Use >= 0.80 and >= 0.90 as thresholds. np.argmax() on a boolean array returns the first True index.

TipSolution 
prop_var = pca.explained_variance_ratio_
cum_var  = np.cumsum(prop_var)
print("Variance per PC:      ", np.round(prop_var, 3))
print("Cumulative variance:  ", np.round(cum_var, 3))
n_80 = int(np.argmax(cum_var >= 0.80)) + 1
print(f"PCs needed for >= 80% variance: {n_80}")
n_90 = int(np.argmax(cum_var >= 0.90)) + 1
print(f"PCs needed for >= 90% variance: {n_90}")

Now let’s visualize the cumulative variance as a bar chart.

View code for plot 
prop_var <- summary(pca_result)$importance["Proportion of Variance", ]
cum_var  <- cumsum(prop_var)
barplot(cum_var, names.arg = paste0("PC", 1:4), ylim = c(0, 1.05),
        col = ifelse(cum_var >= 0.9, "forestgreen", "steelblue"),
        ylab = "Cumulative proportion of variance",
        main = "Cumulative variance explained (folliculogenesis PCA)")
abline(h = 0.80, lty = "dashed", col = "orange", lwd = 2)
abline(h = 0.90, lty = "dashed", col = "red", lwd = 2)
legend("bottomright", legend = c("80% threshold","90% threshold"),
       lty = "dashed", col = c("orange","red"), lwd = 2, bty = "n")
View code for plot 
cum_var_plot = np.cumsum(pca_cv.explained_variance_ratio_)
pc_labels    = [f"PC{i+1}" for i in range(4)]
colors       = ["forestgreen" if v >= 0.9 else "steelblue" for v in cum_var_plot]
plt.figure(figsize=(6, 4))
plt.bar(pc_labels, cum_var_plot, color=colors)
plt.axhline(0.80, color="orange", linestyle="--", lw=2, label="80% threshold")
plt.axhline(0.90, color="red",    linestyle="--", lw=2, label="90% threshold")
plt.ylim(0, 1.05); plt.ylabel("Cumulative proportion of variance")
plt.title("Cumulative variance explained (folliculogenesis PCA)")
plt.legend(loc="lower right"); plt.tight_layout(); plt.show()

Section 5: Scree Plot

A scree plot shows eigenvalues in descending order. Look for the “elbow” where the slope changes sharply. PCs before the elbow are retained.

Exercise: Draw a scree plot

NoteHint

Pass eigenvalues as the first argument to plot(). The Kaiser line sits at h = 1.

TipSolution
eigenvalues <- pca_result$sdev ^ 2 plot(eigenvalues, type = "b", xlab = "Principal Component", ylab = "Eigenvalue", main = "Scree Plot: Folliculogenesis PCA", pch = 19, col = "steelblue", lwd = 2) abline(h = 1, lty = "dashed", col = "red") legend("topright", legend = "Kaiser threshold (eigenvalue = 1)", lty = "dashed", col = "red", bty = "n") 
eigenvalues <- pca_result$sdev ^ 2
plot(eigenvalues, type = "b",
     xlab = "Principal Component", ylab = "Eigenvalue",
     main = "Scree Plot: Folliculogenesis PCA",
     pch = 19, col = "steelblue", lwd = 2)
abline(h = 1, lty = "dashed", col = "red")
legend("topright", legend = "Kaiser threshold (eigenvalue = 1)",
       lty = "dashed", col = "red", bty = "n")
NoteHint

Use pca.explained_variance_ for eigenvalues. The x-axis values are pc_numbers and the Kaiser line is at y = 1.

TipSolution 
eigenvalues_scree = pca.explained_variance_
pc_numbers = list(range(1, len(eigenvalues_scree) + 1))
plt.figure(figsize=(5, 4))
plt.plot(pc_numbers, eigenvalues_scree, "o-", color="steelblue", lw=2, markersize=8)
plt.axhline(y=1, color="red", linestyle="--", label="Kaiser threshold (eigenvalue = 1)")
plt.xlabel("Principal Component"); plt.ylabel("Eigenvalue")
plt.title("Scree Plot: Folliculogenesis PCA")
plt.xticks(pc_numbers); plt.legend(); plt.tight_layout(); plt.show()

Section 6: What Do the Components Mean? Loadings

We have decided to retain 2 PCs. But what do PC1 and PC2 actually mean biologically? The answer lies in the loadings.

Loadings are the weights assigned to each gene in the construction of each PC. A large absolute loading means that gene strongly drives that PC. Genes with loadings near zero have little influence.

TipInterpreting loadings biologically

Look for genes with large absolute loadings. If a group of functionally related genes all load heavily on the same PC, that PC likely represents an underlying biological process, such as the maturation trajectory from primordial to preovulatory follicle.

6.1 Extract and view loadings

NoteHint

Loadings live in pca_result$rotation. Pass that to round() to print neatly.

TipSolution
loadings <- pca_result$rotation print(round(loadings, 3)) 
loadings <- pca_result$rotation
print(round(loadings, 3))
NoteHint

Loadings are in pca.components_. Transpose with .T so rows are genes.

TipSolution 
loadings = pca.components_.T
loadings_df = pd.DataFrame(
    np.round(loadings, 3),
    index=feature_names,
    columns=[f"PC{i+1}" for i in range(4)]
)
print(loadings_df)

6.2 Visualize loadings as a bar chart

View code for plot 
loadings  <- pca_result$rotation
var_names <- rownames(loadings)
par(mfrow = c(1, 2), mar = c(6, 4, 3, 1))
barplot(loadings[,1], names.arg=var_names, main="PC1 Loadings", ylab="Loading",
        col=ifelse(loadings[,1]>0,"steelblue","tomato"), las=2, ylim=c(-1,1))
abline(h=0)
barplot(loadings[,2], names.arg=var_names, main="PC2 Loadings", ylab="Loading",
        col=ifelse(loadings[,2]>0,"steelblue","tomato"), las=2, ylim=c(-1,1))
abline(h=0)
par(mfrow=c(1,1))
View code for plot 
loadings = pca_plot.components_.T
fig, axes = plt.subplots(1, 2, figsize=(9, 4), sharey=True)
for ax, pc_idx in zip(axes, [0, 1]):
    vals = loadings[:, pc_idx]
    cols = ["steelblue" if v > 0 else "tomato" for v in vals]
    ax.bar(feature_names, vals, color=cols)
    ax.axhline(0, color="black", lw=0.8)
    ax.set_title(f"PC{pc_idx+1} Loadings"); ax.set_ylabel("Loading"); ax.set_ylim(-1, 1)
    ax.set_xticks(range(len(feature_names)))
    ax.set_xticklabels(feature_names, rotation=25, ha="right", fontsize=9)
plt.suptitle("Loadings for PC1 and PC2: Folliculogenesis PCA")
plt.tight_layout(); plt.show()

Section 7: Biological Interpretation of the Loadings

Based on the loading bar charts in Section 6:

  • PC1 is dominated by NOBOX (large negative loading) and FSHR (large positive loading), with GDF9 also contributing. Since NOBOX is an early-stage oocyte factor and FSHR is a late-stage granulosa cell receptor, PC1 represents the follicular maturation gradient: samples progress from negative to positive PC1 scores as follicles develop from primordial to preovulatory.

  • PC2 is driven primarily by LHCGR and partly by GDF9. LHCGR spikes specifically at the preovulatory stage in response to the LH surge. PC2 therefore captures a preovulatory transition signal that is orthogonal to the general maturation gradient on PC1.

Note

The sign of a PC is arbitrary. PCA can flip any axis and the result is mathematically equivalent. What matters is the magnitude of loadings and their relative signs across genes.

Section 8: Score Plot

The score plot places each sample in the new PC coordinate system, coloured by follicle stage. This is where the biological interpretation pays off: we can see whether samples from the same stage cluster together, and whether the stages separate along axes that reflect the biology encoded in the loadings.

This is analogous to Figure 1 in Zhang et al. (2018), where PCA of the full transcriptome showed clear separation of oocyte samples by developmental stage.

View code for plot 
library(ggplot2)
scores        <- as.data.frame(pca_result$x)
scores$Stage  <- follicle_data$Stage
pct           <- round(summary(pca_result)$importance[2, 1:2] * 100, 1)
stage_cols    <- c("Primordial"="#4e79a7","Primary"="#f28e2b",
                   "Secondary"="#59a14f","Antral"="#e15759","Preovulatory"="#b07aa1")

ggplot(scores, aes(x=PC1, y=PC2, colour=Stage)) +
  geom_point(size=2.2, alpha=0.85) +
  scale_colour_manual(values=stage_cols) +
  labs(title = "PCA Score Plot: Folliculogenesis (GSE107746-inspired)",
       x = paste0("PC1 (",pct[1],"% variance): maturation gradient"),
       y = paste0("PC2 (",pct[2],"% variance): preovulatory transition")) +
  theme_bw(base_size=12)
View code for plot 
df = pd.DataFrame(scores_arr, columns=["PC1","PC2"])
df["Stage"] = pd.Categorical(stage_labels, categories=stage_order)
stage_cols = {"Primordial":"#4e79a7","Primary":"#f28e2b",
              "Secondary":"#59a14f","Antral":"#e15759","Preovulatory":"#b07aa1"}
var_pct = pca2.explained_variance_ratio_ * 100

fig, ax = plt.subplots(figsize=(7,5))
for sp, grp in df.groupby("Stage", observed=True):
    ax.scatter(grp["PC1"], grp["PC2"], label=sp, color=stage_cols[sp], alpha=0.85, s=35)
ax.set_xlabel(f"PC1 ({var_pct[0]:.1f}% variance): maturation gradient")
ax.set_ylabel(f"PC2 ({var_pct[1]:.1f}% variance): preovulatory transition")
ax.set_title("PCA Score Plot: Folliculogenesis (GSE107746-inspired)")
ax.legend(title="Stage", bbox_to_anchor=(1.02,1), loc="upper left", fontsize=8)
plt.tight_layout(); plt.show()

Notice how the five follicular stages form a roughly ordered sequence along PC1, consistent with the maturation gradient interpretation from the loadings. The Preovulatory samples separate upward on PC2, capturing the LHCGR-driven preovulatory transition signal. This mirrors the kind of stage-level separation that Zhang et al. observed in the full transcriptome data (GSE107746).

Summary

In this module you have learned:

What PCA is and why it is essential for high-dimensional biological data like RNA-seq.

How to build and explore a gene expression dataset based on a published folliculogenesis study (GSE107746).

How to run PCA in R (prcomp) and Python (sklearn.decomposition.PCA) with scaling.

How to apply the Kaiser criterion (eigenvalue > 1) to select components.

How to use cumulative variance explained (>=80-90%) to select components.

How to draw and interpret a scree plot to identify the elbow.

How to extract and interpret loadings to give PCs a biological meaning.

How to produce and interpret a score plot coloured by biological group.

Key takeaway

For this folliculogenesis dataset, all three criteria agree: retain 2 principal components. Together they explain approximately 90% of the total variance in four marker genes. The loadings show that PC1 represents the follicular maturation gradient (NOBOX declining, FSHR rising) and PC2 captures the preovulatory LHCGR spike. This biologically interpretable structure in just two dimensions mirrors what Zhang et al. observed in thousands of genes, demonstrating why PCA is a first-line tool in transcriptomic analysis.

In practice, the three criteria sometimes disagree. Always consider them together and use your biological judgement.

Return to Learning Hub Homepage

Learning Hub Homepage