Why this activity matters

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?

Learning goals

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

Find the data directory

This activity expects the following files:

# 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.

brca_expr_mat <- readRDS(file.path(data_dir,"brca_expr_mat.rds"))

Inspect the matrix.

dim(brca_expr_mat)
## [1] 18351  1082
brca_expr_mat[1:5, 1:5]
##          TCGA-3C-AAAU TCGA-3C-AALI TCGA-3C-AALJ TCGA-3C-AALK TCGA-4H-AAAK
## TSPAN6      7.5636463     7.705439     9.975045    10.110718     9.881960
## TNMD        0.4272843     1.061776     5.353718     1.164271     2.506120
## DPM1        9.0367945     9.662019     9.919403     8.859174     8.928912
## SCYL3       8.3493520    10.607275     8.395877     9.103957     8.704889
## C1orf112    6.9691735     8.106778     7.922947     7.776269     7.495535

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: Genes are commonly stored as rows because researchers often compare each gene across many patient samples. Patients are stored as columns so their expression profiles can be compared side by side.

Important note: these data are already log-transformed

The values in this expression matrix are mostly between 0 and about 21.

summary(as.vector(brca_expr_mat))
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.000   3.669   8.056   6.852   9.875  20.978

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 gene expression counts can range from 0 to thousands or even millions. A log transformation compresses the large values into smaller ranges. That is why most are between 0-21, which show that the original measurements have been transformed.

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.

Load the clinical data

The clinical data contain patient and tumor information.

brca_clin_df <- read.csv(
  file.path(data_dir, "brca_clin.csv"),
  stringsAsFactors = FALSE
)

dim(brca_clin_df)
## [1] 1082   27
head(brca_clin_df[, 1:6])
##   bcr_patient_barcode gender                      race              ethnicity
## 1        TCGA-3C-AAAU FEMALE                     WHITE NOT HISPANIC OR LATINO
## 2        TCGA-3C-AALI FEMALE BLACK OR AFRICAN AMERICAN NOT HISPANIC OR LATINO
## 3        TCGA-3C-AALJ FEMALE BLACK OR AFRICAN AMERICAN NOT HISPANIC OR LATINO
## 4        TCGA-3C-AALK FEMALE BLACK OR AFRICAN AMERICAN NOT HISPANIC OR LATINO
## 5        TCGA-4H-AAAK FEMALE                     WHITE NOT HISPANIC OR LATINO
## 6        TCGA-5L-AAT0 FEMALE                     WHITE     HISPANIC OR LATINO
##   age_at_diagnosis year_of_initial_pathologic_diagnosis
## 1               55                                 2004
## 2               50                                 2003
## 3               62                                 2011
## 4               52                                 2011
## 5               50                                 2013
## 6               42                                 2010

As we saw previously, the clinical data includes receptor status.

table(brca_clin_df$estrogen_receptor_status)
## 
## [Not Evaluated]   Indeterminate        Negative        Positive 
##              48               2             236             796
table(brca_clin_df$progesterone_receptor_status)
## 
## [Not Evaluated]   Indeterminate        Negative        Positive 
##              49               4             340             689
table(brca_clin_df$her2_receptor_status)
## 
## [Not Available] [Not Evaluated]       Equivocal   Indeterminate        Negative 
##               8             170             177              12             554 
##        Positive 
##             161

Make sure samples are aligned

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.

sample_ids <- colnames(brca_expr_mat)

match_index <- match(sample_ids, brca_clin_df$bcr_patient_barcode)

sum(is.na(match_index))
## [1] 0

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.

clin_matched <- brca_clin_df[match_index, ]

all(clin_matched$bcr_patient_barcode == sample_ids)
## [1] TRUE

Now clin_matched is aligned to the columns of brca_expr_mat.

Average expression across samples

For each gene, we can calculate its average expression across all tumors.

mean_expr <- apply(brca_expr_mat, 1, mean)

summary(mean_expr)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.000   3.759   8.036   6.852   9.771  16.470

What does this look like in boxplot form?

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

# Check out the actual values again:
summary(mean_expr)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.000   3.759   8.036   6.852   9.771  16.470

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 appears larger than the median. This suggests that the distribution is right-skewed because a smaller number of genes with very high average expression 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, a gene with expression below the median is not automatically unimportant. Some genes are only needed in certain cell types or conditions, and even genes expressed at low levels can have important regulatory or signaling roles.

If you want to learn more about boxplots (otherwise known as whisker plots) check out this truly awesome Statquest video.

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: The genes are not all expressed at similar levels. Some appear much more active than others. Cells regulate genes differently because different genes perform different functions.

Variance across samples

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.

var_genes <- apply(brca_expr_mat, 1, var)

summary(var_genes)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.0000  0.3211  0.6607  1.2801  1.5830 26.8820

Let’s look at the distribution of variance values.

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:Genes that vary greatly among patients can help reveal differences between tumors, such as cancer subtype, receptor status, mutation patterns, or disease behavior. Genes that barely change provide less information for separating patients into meaningful groups.

Select the most variable genes

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: I think it would make the heatmap extremely crowded and difficult to read. Important patterns would be hidden by genes that change very little, and the plot would also take more time and computer memory to create.

order_var <- order(var_genes, decreasing = TRUE)

num_genes <- 100

expr_top <- brca_expr_mat[order_var[1:num_genes], ]

dim(expr_top)
## [1]  100 1082

Select a subset of samples

There are many patient samples. For an introductory heatmap, we will plot every fourth sample.

our_samples <- seq(1, ncol(expr_top), by = 4)

expr_sub <- expr_top[, our_samples]

clin_sub <- clin_matched[our_samples, ]

dim(expr_sub)
## [1] 100 271

The most important heatmap correction: scale genes, not samples

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.

expr_sub_scaled <- t(scale(t(expr_sub)))

# Replace any NA values that could occur for genes with zero variance
expr_sub_scaled[is.na(expr_sub_scaled)] <- 0

summary(as.vector(expr_sub_scaled))
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## -3.31947 -0.76921 -0.07953  0.00000  0.71035  3.85797

Now each gene is shown relative to its own average across samples.

**Reflection: Before scaling, some genes naturally have much higher expression than others. After scaling, what does the color represent?

Your answer: After scaling, the color represents how high or low each gene’s expression is compared with that same gene’s average across the selected tumors. Red means above the gene’s average, blue means below its average, and white means close to its average.

Heatmap of variable genes

We will use a blue-white-red palette.

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:es. Selecting the most variable genes reduces unnecessary information and emphasizes genes that differ the most between tumors. This makes clusters and contrasting expression patterns easier to see than they would be in a heatmap containing every gene.

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:I predict that the two tumors share a similar molecular subtype or related biological characteristics. I would also need clinical information such as receptor status, pathologic stage, histology, mutations, treatment response, and patient outcome before concluding that they are biologically similar.

Add estrogen receptor status

Now we will add clinical labels.

We use:

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)
## er_label
##       .   + 
##  13  57 201

Heatmap with ER labels

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: The ER-negative tumors appear somewhat concentrated within particular branches or sections of the heatmap, although the separation is not perfect and some ER-negative samples are mixed with ER-positive tumors.

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. A perfect match would show a strong association, but it would not prove causation. ER status may be connected with a broader breast-cancer subtype, mutations, or other biological factors that also influence gene expression. A controlled experiment and additional evidence would be needed to establish a causal relationship.

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.

Marker gene heatmap

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:

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
## [1] "ESR1"  "PGR"   "ERBB2" "FOXA1" "KRT5"  "KRT14" "MKI67" "EPCAM"
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: ER-positive tumors show higher relative expression of ESR1, while ER-negative tumors show lower relative ESR1 expression. This makes biological sense because ESR1 contains the instructions for the estrogen receptor.

Reflection: Why is this heatmap easier to interpret than the heatmap containing 100 genes?

Your answer: This heatmap is easier to interpret because it contains only eight labeled genes with known roles in breast-cancer biology. Instead of examining 100 unnamed patterns, I can directly compare specific markers such as ESR1, PGR, ERBB2, and the basal-like keratin genes.

This marker-gene heatmap may be easier to explain than the larger unsupervised heatmap.

CHALLENGE 1: PR status

Create labels for progesterone receptor status.

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)
## pr_label
##       .   + 
##  13  81 177

Now plot the heatmap with PR labels.

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: PR status looks similar to ER status, with many PR-positive and ER-positive tumors showing related clustering patterns. However, the PR separation appears less complete, and some PR-positive and PR-negative tumors are mixed together.

CHALLENGE 2: HER2 status

Create labels for HER2 status.

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)
## her2_label
##       .   + 
## 108 125  38

Now plot the heatmap with HER2 labels.

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:HER2 status does not appear to separate as clearly as ER status.

Reflection: Which receptor (ER, PR, or HER2) seems to show the strongest relationship with gene expression patterns? Were you surprised?

Your answer: ER appears to show the strongest relationship with the overall gene-expression patterns, with PR showing a related but weaker pattern and HER2 showing the least clear separation. I was not very surprised because ER status is connected with a broad luminal gene-expression program, not only the expression of a single gene.

CHALLENGE 3: triple-negative breast cancer

Triple-negative breast cancer means:

tnbc <- er_status == "Negative" &
  pr_status == "Negative" &
  her2_status == "Negative"

tnbc_label <- rep("", length(tnbc))
tnbc_label[tnbc] <- "TN"

table(tnbc_label)
## tnbc_label
##      TN 
## 244  27

Now plot the heatmap with TN labels.

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:The triple-negative tumors do not appear as one perfectly clean group. Some are located near one another, but others are mixed with non-triple-negative tumors, showing that triple-negative breast cancer is not a single uniform molecular category.

Reflection: If triple-negative tumors do not all cluster together, what are two possible biological explanations?

Your answer: First, triple-negative breast cancer is biologically diverse and includes multiple molecular subtypes with different mutations and expression programs. Second, each tumor sample contains a mixture of cancer cells and surrounding stromal or immune cells, which can influence the measured gene-expression pattern. Technical noise and differences in tumor purity could also contribute.

Final reflection

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: One conclusion I feel confident making is that gene-expression patterns are related to clinically important breast-cancer features, especially estrogen receptor status, although the groups do not separate perfectly. Next, I would investigate which individual genes or pathways best distinguish ER-negative and triple-negative tumors and whether those patterns are connected with treatment response or patient outcomes.

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.

Knit your report

Click:

Knit → Knit to HTML

Your final report should include: