In the previous heatmap activity, we used the small built-in
mtcars dataset. That was useful because the dataset was
small enough to see clearly.
Now we will use real gene expression data from TCGA breast cancer samples. Gene expression data which tells us how many messenger RNAs (mRNAs) per gene are present in a patient sample. The amount of a gene’s mRNA corresponds (roughly) to the amount of protein in the sample.
This is more realistic, but also more challenging:
That is normal in real computational biology.
Our goal is to use heatmaps to ask:
Do breast tumors with similar gene expression patterns also share clinical features, such as estrogen receptor status?
By the end of this activity, you should be able to:
This activity expects the following files:
brca_expr_matbrca_clin.csv```{r setup, message=FALSE, warning=FALSE} # This chunk sets up file path for the activity.
data_dir <- “/shared/dreamhigh/data”
---
## Load the expression data
The expression file is an RDS file.
An RDS file (which ends in .rds) is a special file format used by R to save exactly one specific piece of data (like a single data table, a list, or a machine learning model) from your computer's memory onto your hard drive.
Reading and writing RDS files is significantly faster than processing text-based files.
Rows are genes.
Columns are patient tumor samples.
```{r load-expression}
brca_expr_mat <- readRDS(file.path(data_dir,"brca_expr_mat.rds"))
Inspect the matrix.
```{r inspect-expression} dim(brca_expr_mat)
brca_expr_mat[1:5, 1:5]
**Question:** What do the rows represent? What do the columns represent?
> Your answer: Rows are genes; columns are patient tumor samples.
**Reflection:** Why do you think genes are stored as rows and patients as columns? Could the data have been organized the other way around?
> Your answer: Storing genes (~20,000) as rows and patients (hundreds) as columns is standard convention in bioinformatics. Data could be transposed, but standard tools expect this layout.
---
### Important note: these data are already log-transformed
The values in this expression matrix are mostly between 0 and about 21.
```{r expression-range}
summary(as.vector(brca_expr_mat))Reflection: The largest expression values are only around 20 instead of thousands or millions. Why does this suggest the data have already been log-transformed?
Your answer: Raw counts reach millions. Log transformation compresses numbers into a 0–20 scale, preventing massive raw values from distorting analysis.
The distribution of values is a strong clue that these values are already on a transformed scale, likely a log-like expression scale.
We log-transform gene expression data to make highly skewed numbers more symmetrical. This fixes a common problem where a few highly active genes distort your and data dominate statistical analyses purely due to their massive raw numerical values, rather than their actual biological relevance.
The clinical data contain patient and tumor information.
```{r load-clinical} brca_clin_df <- read.csv( file.path(data_dir, “brca_clin.csv”), stringsAsFactors = FALSE )
dim(brca_clin_df) head(brca_clin_df[, 1:6])
As we saw previously, the clinical data includes receptor status.
```{r receptor-tables}
table(brca_clin_df$estrogen_receptor_status)
table(brca_clin_df$progesterone_receptor_status)
table(brca_clin_df$her2_receptor_status)This is a very important step.
The expression matrix columns are sample IDs.
The clinical data rows are patient/sample IDs.
We should match them by name, not just assume they are in the same order.
```{r match-samples} sample_ids <- colnames(brca_expr_mat)
match_index <- match(sample_ids, brca_clin_df$bcr_patient_barcode)
sum(is.na(match_index))
If the result is 0, which should be the case here, every expression sample matched a clinical row.
If the result wasn't zero, we can use `match_index` to sort the rows of the clinical data to match the columns of the expression data.
```{r align-clinical}
clin_matched <- brca_clin_df[match_index, ]
all(clin_matched$bcr_patient_barcode == sample_ids)Now clin_matched is aligned to the columns of
brca_expr_mat.
For each gene, we can calculate its average expression across all tumors.
```{r mean-expression} mean_expr <- apply(brca_expr_mat, 1, mean)
summary(mean_expr)
What does this look like in boxplot form?
```{r}
mean_expr <- apply(brca_expr_mat, 1, mean)
# Draw boxplot but don't plot outliers
bp <- boxplot(mean_expr,
boxwex = 0.35,
horizontal = TRUE,
col = "lightblue",
outline = FALSE,
main = "Distribution of Mean Gene Expression",
xlab = "Mean Expression")
# Add the mean as a red diamond
mean.val <- mean(mean_expr)
points(mean.val, 1, pch = 23, bg = "red", cex = 1.5)
# Label the mean
text(mean.val, 1.15,
labels = paste0("Mean = ", sprintf("%.2f", mean.val)),
col = "red")
# Label the five-number summary
stats <- bp$stats
text(stats[1], 0.82, sprintf("%.1f", stats[1])) # Min
text(stats[2] - 0.10, 0.82, sprintf("%.1f", stats[2])) # Q1
text(stats[3], 0.82, sprintf("%.1f", stats[3])) # Median
text(stats[4] + 0.10, 0.82, sprintf("%.1f", stats[4])) # Q3
text(stats[5], 0.82, sprintf("%.1f", stats[5])) # Max
{r} # Check out the actual values again: summary(mean_expr)
Reflection: Is the mean (red diamond) larger or smaller than the median (black line)? What does this tell you about the distribution of gene expression?
(Clue: Skewness measures the asymmetry of data. In a symmetrical distribution, the mean and median are identical. In an asymmetrical (skewed) distribution, extreme values or a long tail “pull” the mean toward the direction of the tail, while the median remains closer to the center of the data.)
Your answer: The mean is larger than the median, showing a right-skewed distribution where a few highly expressed genes pull the mean upward.
Reflection: About half of the genes have an average expression below the median. Does that mean half of the genes are “unimportant”? Why or why not?
Your answer: No. Crucial regulatory proteins (like transcription factors) work at low expression levels, while high-expression genes are often basic housekeeping genes.
If you want to learn more about boxplots (otherwise known as whisker plots) check out this truly awesome Statquest video.
{r mean-expression-hist} hist( mean_expr, breaks = 50, main = "Mean gene expression across breast tumors", xlab = "Mean expression" )
Reflection: Do all genes appear to be expressed at similar levels, or do some genes appear much more active than others? Why might cells regulate genes differently?
Your answer: Cells adjust gene expression based on tissue type, function, metabolic needs, and external signals.
A gene can have a high average expression but not vary much between patients.
For heatmaps, genes that vary across patients are often more informative.
```{r gene-variance} var_genes <- apply(brca_expr_mat, 1, var)
summary(var_genes)
Let's look at the distribution of variance values.
```{r variance-hist}
hist(
var_genes,
breaks = 50,
main = "Variance of gene expression across breast tumors",
xlab = "Variance"
)Most genes have relatively low variance. *
Reflection: Why might genes that change a lot from patient to patient be more useful for studying cancer than genes whose expression hardly changes?
Your answer: High-variance genes highlight differences between tumor subtypes and patient outcomes. Invariant genes look identical across all patients and offer no useful contrast.
We will begin with a small number of highly variable genes.
This is easier to interpret than trying to plot all genes at once.
Prediction: What do you think would happen if we plotted all 20,000 genes instead of only the 100 most variable genes?
Your answer:It would create visual clutter and introduce background noise, hiding important clustering patterns.
```{r select-variable-genes} order_var <- order(var_genes, decreasing = TRUE)
num_genes <- 100
expr_top <- brca_expr_mat[order_var[1:num_genes], ]
dim(expr_top)
---
### Select a subset of samples
There are many patient samples. For an introductory heatmap, we will plot every fourth sample.
```{r select-samples}
our_samples <- seq(1, ncol(expr_top), by = 4)
expr_sub <- expr_top[, our_samples]
clin_sub <- clin_matched[our_samples, ]
dim(expr_sub)This is the key idea.
For a gene expression heatmap, we usually want to ask:
Is each gene higher or lower than its own average across tumors?
That means we should scale each row of the matrix, because rows are genes.
Base R’s scale() function scales columns by default.
Since our columns are patients, this would scale patients, not
genes.
So we use t(scale(t(matrix))).
This transposes the matrix, scales the genes, and transposes it back.
```{r scale-genes} expr_sub_scaled <- t(scale(t(expr_sub)))
expr_sub_scaled[is.na(expr_sub_scaled)] <- 0
summary(as.vector(expr_sub_scaled))
Now each gene is shown relative to its own average across samples.
- red = expression is higher than that gene's average
- blue = expression is lower than that gene's average
- white = expression near that gene's average
**Reflection^^: Before scaling, some genes naturally have much higher expression than others. After scaling, what does the color represent?
> Your answer: Color shows expression relative to each gene's own average: Red = above average, Blue = below average, White = near average.
---
## Heatmap of variable genes
We will use a blue-white-red palette.
```{r heatmap-variable-genes, fig.width=9, fig.height=8}
heat_colors <- colorRampPalette(c("blue", "white", "red"))(100)
heatmap(
expr_sub_scaled,
labRow = "",
labCol = "",
margins = c(3, 3),
xlab = "Tumor samples",
ylab = "Variable genes",
col = heat_colors,
zlim = c(-2, 2),
main = "Most variable genes in TCGA breast cancer samples"
)Reflection: Do you think selecting only the most variable genes helped make patterns easier to see? Explain your reasoning.
Your answer: Yes. Filtering out non-variable genes removes background noise and brings main tumor patterns into focus.
Reflection: If you saw two tumors with nearly identical expression patterns, what might you predict about those tumors? What additional information would you need before concluding they are biologically similar?
Your answer: You would predict they share a subtype or clinical outcome. You would need DNA mutation data, protein levels, and clinical treatment responses to confirm.
Now we will add clinical labels.
We use:
+ for ER-positive. for ER-negative```{r er-labels} er_status <- clin_sub$estrogen_receptor_status
er_label <- rep(““, length(er_status)) er_label[er_status == ”Positive”] <-”+” er_label[er_status == “Negative”] <- “.”
table(er_label)
---
### Heatmap with ER labels
```{r heatmap-er-labels, fig.width=9, fig.height=8}
heatmap(
expr_sub_scaled,
labRow = "",
labCol = er_label,
cexCol = 0.5,
margins = c(3, 3),
xlab = "Tumor samples labeled by ER status",
ylab = "Variable genes",
col = heat_colors,
zlim = c(-2, 2),
main = "Gene expression heatmap with ER status labels"
)Interpretation question: Do ER-negative tumors appear concentrated in any part of the heatmap?
Your answer: Yes. ER-negative tumors group together with lower expression of luminal genes and higher expression of basal genes.
Reflection: Suppose the ER labels matched the heatmap perfectly. Would that prove that estrogen receptor status causes the expression patterns? Why or why not?
Your answer: No. Correlation does not equal causation. Observed patterns could stem from downstream signaling, cell-of-origin traits, or shared mutations.
Careful science note:
It is okay if the separation is not perfect. Real tumor data are
complex. We are looking for patterns, not expecting every sample to
behave perfectly.
Sometimes a small set of biologically meaningful genes is easier to interpret than the top 100 variable genes.
Here are several genes related to breast cancer subtype or tumor biology:
ESR1: estrogen receptorPGR: progesterone receptorERBB2: HER2FOXA1: luminal breast cancer biologyKRT5, KRT14: basal-like featuresMKI67: proliferationEPCAM: epithelial marker```{r marker-genes} marker_genes <- c(“ESR1”, “PGR”, “ERBB2”, “FOXA1”, “KRT5”, “KRT14”, “MKI67”, “EPCAM”)
marker_genes <- marker_genes[marker_genes %in% rownames(brca_expr_mat)]
marker_mat <- brca_expr_mat[marker_genes, our_samples]
marker_scaled <- t(scale(t(marker_mat))) marker_scaled[is.na(marker_scaled)] <- 0
marker_genes
```{r marker-heatmap, fig.width=9, fig.height=5}
heatmap(
marker_scaled,
labRow = rownames(marker_scaled),
labCol = er_label,
cexRow = 0.9,
cexCol = 0.5,
margins = c(4, 8),
xlab = "Tumor samples labeled by ER status",
ylab = "Marker genes",
col = heat_colors,
zlim = c(-2, 2),
main = "Breast cancer marker genes"
)Question: How does ESR1 expression
relate to ER status?
Your answer: ESR1 encodes Estrogen Receptor Alpha. ER-positive tumors show high ESR1 expression (red); ER-negative tumors show low expression (blue).
Reflection: Why is this heatmap easier to interpret than the heatmap containing 100 genes?
Your answer: It uses a small group of well-studied genes, making patterns easy to match directly with known breast cancer biology.
This marker-gene heatmap may be easier to explain than the larger unsupervised heatmap.
Create labels for progesterone receptor status.
```{r pr-labels} pr_status <- clin_sub$progesterone_receptor_status
pr_label <- rep(““, length(pr_status)) pr_label[pr_status == ”Positive”] <-”+” pr_label[pr_status == “Negative”] <- “.”
table(pr_label)
Now plot the heatmap with PR labels.
```{r heatmap-pr-labels, fig.width=9, fig.height=8}
heatmap(
expr_sub_scaled,
labRow = "",
labCol = pr_label,
cexCol = 0.5,
margins = c(3, 3),
xlab = "Tumor samples labeled by PR status",
ylab = "Variable genes",
col = heat_colors,
zlim = c(-2, 2),
main = "Gene expression heatmap with PR status labels"
)Question: Does PR status look similar to ER status?
Your answer: Yes. PGR is driven by ER signaling, so PR status closely follows ER status.
Create labels for HER2 status.
```{r her2-labels} her2_status <- clin_sub$her2_receptor_status
her2_label <- rep(““, length(her2_status)) her2_label[her2_status == ”Positive”] <-”+” her2_label[her2_status == “Negative”] <- “.”
table(her2_label)
Now plot the heatmap with HER2 labels.
```{r heatmap-her2-labels, fig.width=9, fig.height=8}
heatmap(
expr_sub_scaled,
labRow = "",
labCol = her2_label,
cexCol = 0.5,
margins = c(3, 3),
xlab = "Tumor samples labeled by HER2 status",
ylab = "Variable genes",
col = heat_colors,
zlim = c(-2, 2),
main = "Gene expression heatmap with HER2 status labels"
)Question: Does HER2 status separate as clearly as ER status?
Your answer: No. HER2 overexpression stems from localized gene amplification rather than widespread transcriptomic shifts, so it separates less cleanly on broad heatmaps.
Reflection: Which receptor (ER, PR, or HER2) seems to show the strongest relationship with gene expression patterns? Were you surprised?
Your answer:ER status. It marks the core biological split between luminal and non-luminal breast cancer lineages.
Triple-negative breast cancer means:
```{r triple-negative} tnbc <- er_status == “Negative” & pr_status == “Negative” & her2_status == “Negative”
tnbc_label <- rep(““, length(tnbc)) tnbc_label[tnbc] <-”TN”
table(tnbc_label)
Now plot the heatmap with TN labels.
```{r heatmap-tnbc-labels, fig.width=9, fig.height=8}
heatmap(
expr_sub_scaled,
labRow = "",
labCol = tnbc_label,
cexCol = 0.5,
margins = c(3, 3),
xlab = "Tumor samples labeled by triple-negative status",
ylab = "Variable genes",
col = heat_colors,
zlim = c(-2, 2),
main = "Gene expression heatmap with triple-negative labels"
)Question: Do triple-negative tumors appear as one clean group, or are they mixed with other tumors?
Your answer: They mostly cluster together in one main group, but display sub-clustering due to biological variation.
Reflection: If triple-negative tumors do not all cluster together, what are two possible biological explanations?
Your answer: TNBC consists of multiple sub-types, and variations in sample purity or immune cell presence affect overall expression.
Final Reflection: Imagine you are a cancer researcher who has never seen these data before. Based on today’s analyses, what is one conclusion you feel confident making, and what is one question you would want to investigate next?
Your answer: Conclusion: ER status is a primary driver of gene expression patterns in breast cancer. Next Question: What molecular sub-types exist within triple-negative tumors, and how do they impact treatment response?
The most important biological lesson is:
Gene expression patterns can reflect important tumor features, but real cancer data are complex and must be interpreted carefully.
Click:
Knit → Knit to HTML
Your final report should include: