Why this activity matters

In previous DREAM-High activities, we studied cancer patient data from TCGA. Here we switch to a different kind of cancer model: human cancer cell lines.

Cancer cell lines are cells that can grow in the laboratory. Researchers use them to study cancer biology and to test hypotheses about cancer behavior.

In this activity, we will compare two breast cancer cell lines:

T-47D
MDA-MB-231

We will ask:

Do differences in cell movement correspond to differences in gene expression?

This is a core idea in systems biology: we connect a measurable behavior, or phenotype, to molecular data.

Learning goals

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

load cell line motility and gene expression data
inspect data frames and matrices
subset a data frame using biological criteria
match sample IDs across two datasets
calculate log-transformed expression values
calculate differential gene expression between two cell lines
interpret genes in the tails of a distribution
connect gene expression patterns to cell behavior

Setting the data directory

This activity expects the following files:

cell_speeds.csv
pson_expr_mat.csv

The files should be located in the shared data directory.

# This chunk sets up file path for the activity.

data_dir <- "/shared/dreamhigh/data"

Load the data

We will load two files.

The first file contains measurements of how fast cells move.

cell_speeds_df <- read.csv(
  file.path(data_dir, "cell_speeds.csv"),
  stringsAsFactors = FALSE
)

dim(cell_speeds_df)

## [1] 62  7

head(cell_speeds_df)

##     sample cellLine    diagnosis             experimentalCondition
## 1 mRNA_R17    SW620 Colon Cancer                             Glass
## 2 mRNA_R21    SW620 Colon Cancer           HyaluronicAcid Collagen
## 3 mRNA_R20    SW620 Colon Cancer        HyaluronicAcid Fibronectin
## 4 mRNA_R19    SW620 Colon Cancer    30 kPa polyacrylamide Collagen
## 5 mRNA_R18    SW620 Colon Cancer 30 kPa polyacrylamide Fibronectin
## 6 mRNA_R16    SW620 Colon Cancer    500 Pa polyacrylamide Collagen
##   summary_metric average_value total_number_of_cells_tracked
## 1    speed_um_hr      13.93244                            25
## 2    speed_um_hr      18.47756                            51
## 3    speed_um_hr      53.00602                            36
## 4    speed_um_hr      69.78858                            25
## 5    speed_um_hr      27.07560                            25
## 6    speed_um_hr      35.99601                            51

The second file contains gene expression measurements.

pson_expr_df <- read.csv(
  file.path(data_dir, "pson_expr_mat.csv"),
  check.names = FALSE,
  stringsAsFactors = FALSE
)

dim(pson_expr_df)

## [1] 18682    64

pson_expr_df[1:5, 1:5]

##            mRNA_R17 mRNA_R21 mRNA_R20 mRNA_R19
## 1   TSPAN6    33.56    45.10    39.42    35.61
## 2     TNMD     0.00     0.00     0.00     0.00
## 3     DPM1   169.46   129.88   132.06   113.73
## 4    SCYL3     1.85     1.85     1.77     1.93
## 5 C1orf112     5.73    11.85    10.16     9.28

What is in the cell speed data?

Look at the column names.

colnames(cell_speeds_df)

## [1] "sample"                        "cellLine"                     
## [3] "diagnosis"                     "experimentalCondition"        
## [5] "summary_metric"                "average_value"                
## [7] "total_number_of_cells_tracked"

Each row represents a cell line experiment.

Important columns include:

sample: experiment/sample ID
cellLine: name of the cell line
diagnosis: cancer type
experimentalCondition: surface or substrate on which cells were grown
average_value: average cell speed
total_number_of_cells_tracked: number of tracked cells

What cancer types are represented?

table(cell_speeds_df$diagnosis)

## 
##    Brain Cancer   Breast Cancer    Colon Cancer  Not Applicable Prostate Cancer 
##              14              14              13               7               7 
##     Skin Cancer 
##               7

What cell lines model the cancer types?

unique(cell_speeds_df[, c("cellLine", "diagnosis")])

##      cellLine       diagnosis
## 1       SW620    Colon Cancer
## 8       SW480    Colon Cancer
## 14     RWPE-1  Not Applicable
## 21       A375     Skin Cancer
## 28       T98G    Brain Cancer
## 35      22Rv1 Prostate Cancer
## 42      T-47D   Breast Cancer
## 49       U-87    Brain Cancer
## 56 MDA-MB-231   Breast Cancer

Question: Which cell lines represent breast cancer?

Your answer: MDA MB 231 and T-47D cell lines

Focus on breast cancer cell lines

brca_speeds_df <- cell_speeds_df[
  cell_speeds_df$diagnosis == "Breast Cancer",
]

brca_speeds_df

##      sample   cellLine     diagnosis             experimentalCondition
## 42 mRNA_R52      T-47D Breast Cancer                             Glass
## 43 mRNA_R56      T-47D Breast Cancer           HyaluronicAcid Collagen
## 44 mRNA_R55      T-47D Breast Cancer        HyaluronicAcid Fibronectin
## 45 mRNA_R54      T-47D Breast Cancer    30 kPa polyacrylamide Collagen
## 46 mRNA_R53      T-47D Breast Cancer 30 kPa polyacrylamide Fibronectin
## 47 mRNA_R51      T-47D Breast Cancer    500 Pa polyacrylamide Collagen
## 48 mRNA_R50      T-47D Breast Cancer 500 Pa polyacrylamide Fibronectin
## 56 mRNA_R59 MDA-MB-231 Breast Cancer                             Glass
## 57 mRNA_R63 MDA-MB-231 Breast Cancer           HyaluronicAcid Collagen
## 58 mRNA_R62 MDA-MB-231 Breast Cancer        HyaluronicAcid Fibronectin
## 59 mRNA_R61 MDA-MB-231 Breast Cancer    30 kPa polyacrylamide Collagen
## 60 mRNA_R60 MDA-MB-231 Breast Cancer 30 kPa polyacrylamide Fibronectin
## 61 mRNA_R58 MDA-MB-231 Breast Cancer    500 Pa polyacrylamide Collagen
## 62 mRNA_R57 MDA-MB-231 Breast Cancer 500 Pa polyacrylamide Fibronectin
##    summary_metric average_value total_number_of_cells_tracked
## 42    speed_um_hr     13.770019                            15
## 43    speed_um_hr     16.515478                            81
## 44    speed_um_hr     11.288006                            27
## 45    speed_um_hr      7.475378                            30
## 46    speed_um_hr      8.923515                            56
## 47    speed_um_hr     10.608637                            31
## 48    speed_um_hr     16.734859                            47
## 56    speed_um_hr     22.419595                            34
## 57    speed_um_hr     36.196133                            12
## 58    speed_um_hr     55.548778                            26
## 59    speed_um_hr     31.749973                            27
## 60    speed_um_hr     33.378558                            25
## 61    speed_um_hr     35.848219                            30
## 62    speed_um_hr     34.222731                            29

There are two breast cancer cell lines in this dataset:

T-47D
MDA-MB-231

Both are breast cancer cell lines. MDA-MB-231 is used as a model for triple negative breast cancer. T-47D model ER+/PR+/HER2- breast cancer and is useful as a less motile comparison to MDA-MB-231.

Compare experimental conditions

Cell behavior can depend strongly on the surface or substrate on which cells are grown.

table(cell_speeds_df$experimentalCondition)

## 
##    30 kPa polyacrylamide Collagen 30 kPa polyacrylamide Fibronectin 
##                                 9                                 9 
##    500 Pa polyacrylamide Collagen 500 Pa polyacrylamide Fibronectin 
##                                 9                                 9 
##                             Glass           HyaluronicAcid Collagen 
##                                 8                                 9 
##        HyaluronicAcid Fibronectin 
##                                 9

We will focus on the condition HyaluronicAcid Collagen.

Question Why is a substrate containing Hyaluronic Acid and Collagen a good mimic for the cellular environment?

Your answer: It is a good mimic for the cellular environment because it has several properties similar to the natural environment of a cell in the human body, which allows us to accurately study cell behavior in vitro.

Select breast cancer experiments on hyaluronic acid/collagen

hyal_brca_df <- cell_speeds_df[
  cell_speeds_df$diagnosis == "Breast Cancer" &
    cell_speeds_df$experimentalCondition == "HyaluronicAcid Collagen",
]

hyal_brca_df

##      sample   cellLine     diagnosis   experimentalCondition summary_metric
## 43 mRNA_R56      T-47D Breast Cancer HyaluronicAcid Collagen    speed_um_hr
## 57 mRNA_R63 MDA-MB-231 Breast Cancer HyaluronicAcid Collagen    speed_um_hr
##    average_value total_number_of_cells_tracked
## 43      16.51548                            81
## 57      36.19613                            12

This gives us one T-47D experiment and one MDA-MB-231 experiment under the same substrate condition.

Question Look online: Where do T-47D and MDA-MB-231 cells come from? > Your answer: They come from human breast cancer, like infiltrating ductal cancer and adenocarcinoma that may be metastatic.

Visualize cell speed

barplot(
  hyal_brca_df$average_value,
  names.arg = hyal_brca_df$cellLine,
  ylab = "Average speed (microns per hour)",
  main = "Breast cancer cell speed on hyaluronic acid/collagen",
  las = 1
)

We can extract values from the dataframe as well.

hyal_brca_df[, c("sample", "cellLine", "average_value", "total_number_of_cells_tracked")]

##      sample   cellLine average_value total_number_of_cells_tracked
## 43 mRNA_R56      T-47D      16.51548                            81
## 57 mRNA_R63 MDA-MB-231      36.19613                            12

Question: Which breast cancer cell line moves faster in this experiment?

Your answer:MDA MB-231

Prepare the expression matrix

The expression file contains gene names in the first column and expression values in the remaining columns.

gene_names <- pson_expr_df[, 1]
gene_names <- make.unique(gene_names)

pson_expr_mat <- as.matrix(pson_expr_df[, -1])
storage.mode(pson_expr_mat) <- "numeric"

rownames(pson_expr_mat) <- gene_names

dim(pson_expr_mat)

## [1] 18682    63

pson_expr_mat[1:5, 1:5]

##          mRNA_R17 mRNA_R21 mRNA_R20 mRNA_R19 mRNA_R18
## TSPAN6      33.56    45.10    39.42    35.61    32.10
## TNMD         0.00     0.00     0.00     0.00     0.00
## DPM1       169.46   129.88   132.06   113.73   123.37
## SCYL3        1.85     1.85     1.77     1.93     2.30
## C1orf112     5.73    11.85    10.16     9.28    10.01

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

Your answer: The rows represent gene names and the columns represent the mRNA expression values.

Check expression scale

These expression values are TPM-like abundance values. They include zeros and some very large values.

summary(as.vector(pson_expr_mat))

##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
##     0.00     0.06     4.34    50.36    23.32 48098.71

Question: Do these values seem like they have been log-transformed yet? If not, why?

Your answer: No, because usually the values that are log transformed are under 20, and here, the max is much greater than many of the other dataset values.

We learned in the TCGA activity that a log transformation will make large differences easier to visualize and compare.

pson_log_mat <- log2(pson_expr_mat + 1)

summary(as.vector(pson_log_mat))

##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
##  0.00000  0.08406  2.41684  2.68849  4.60407 15.55374

We add 1 before taking the log so that genes with expression 0 do not become -Inf.

We can visualize the distribution with a barplot as we did previously.

mean_expr <- apply(pson_log_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

Question Is this distribution “right” or “left” skewed? A “right skewed” distribution is one where many large values are out in the right tail, so the mean has a higher value than the median.

Your answer: Yes, this is right skewed because the mean is greater than the median.

Match motility samples to expression samples

The speed table tells us which expression samples correspond to our selected cell line experiments.

hyal_brca_df$sample

## [1] "mRNA_R56" "mRNA_R63"

colnames(pson_log_mat)[1:10]

##  [1] "mRNA_R17" "mRNA_R21" "mRNA_R20" "mRNA_R19" "mRNA_R18" "mRNA_R16"
##  [7] "mRNA_R15" "mRNA_R38" "mRNA_R42" "mRNA_R41"

sample_index <- match(hyal_brca_df$sample, colnames(pson_log_mat))

sample_index

## [1] 44 58

Check that every selected sample was found.

sum(is.na(sample_index))

## [1] 0

If this number is 0, the match worked.

Extract expression for the two breast cancer samples

hyal_brca_log_mat <- pson_log_mat[, sample_index]

colnames(hyal_brca_log_mat) <- hyal_brca_df$cellLine

hyal_brca_log_mat[1:5, ]

##             T-47D MDA-MB-231
## TSPAN6   4.224195   4.060912
## TNMD     0.000000   0.000000
## DPM1     6.959654   6.738903
## SCYL3    3.288359   2.634593
## C1orf112 3.493135   4.095924

For this activity, we will define:

T-47D as the slower cell line in this condition
MDA-MB-231 as the faster cell line in this condition

slow_cell <- hyal_brca_df$cellLine[which.min(hyal_brca_df$average_value)]
fast_cell <- hyal_brca_df$cellLine[which.max(hyal_brca_df$average_value)]

slow_cell

## [1] "T-47D"

fast_cell

## [1] "MDA-MB-231"

Question Based on your examination of the cell lines in a question above, in addition to the speed data, why do you think we call MDA-MD-231 the fast cell line?

Your answer: Because it has higher average speeds in the cell lines than the T-47D cell line.

Differential expression: fast minus slow

We now calculate a simple difference: expression in faster cell line - expression in slower cell line.

Genes with large positive values are higher in the faster cell line.

Genes with large negative values are higher in the slower cell line.

dge <- hyal_brca_log_mat[, fast_cell] - hyal_brca_log_mat[, slow_cell]

DGE_df <- data.frame(
  gene = rownames(hyal_brca_log_mat),
  slow_expression = hyal_brca_log_mat[, slow_cell],
  fast_expression = hyal_brca_log_mat[, fast_cell],
  fast_minus_slow = dge,
  stringsAsFactors = FALSE
)

head(DGE_df)

##              gene slow_expression fast_expression fast_minus_slow
## TSPAN6     TSPAN6        4.224195      4.06091205     -0.16328261
## TNMD         TNMD        0.000000      0.00000000      0.00000000
## DPM1         DPM1        6.959654      6.73890291     -0.22075134
## SCYL3       SCYL3        3.288359      2.63459327     -0.65376529
## C1orf112 C1orf112        3.493135      4.09592442      0.60278950
## FGR           FGR        0.000000      0.07038933      0.07038933

Distribution of differential expression

hist(
  DGE_df$fast_minus_slow,
  breaks = 60,
  main = "Differential gene expression: faster minus slower cell line",
  xlab = "log2 expression difference",
  ylab = "Number of genes"
)

abline(v = 0, lwd = 2)
abline(v = c(-5, 5), lty = 2)

Question: Are most genes very different between the two cell lines, or are most genes closer to zero?

Your answer: Most genes are closer to zero between the two cell lines

Volcano-style effect plot

A standard volcano plot displays the gene-expression difference on the x-axis and statistical significance, usually -log10(p-value), on the y-axis. However, this comparison contains only one expression sample for each cell line. With no biological replicates, we cannot estimate sample-to-sample variability or calculate valid p-values.

Instead, we will make a volcano-style effect plot. The x-axis shows the expression difference between the faster and slower cell lines, while the y-axis shows the average expression of each gene across the two cell lines. This lets us distinguish large differences involving well-expressed genes from large differences involving genes expressed at very low levels.

# Average log2 expression across the two cell lines
DGE_df$mean_expression <- rowMeans(
  DGE_df[, c("slow_expression", "fast_expression")]
)

# Classify genes using the same exploratory cutoff used above
DGE_df$expression_group <- "Similar expression"
DGE_df$expression_group[DGE_df$fast_minus_slow > 5] <- "Higher in faster cell line"
DGE_df$expression_group[DGE_df$fast_minus_slow < -5] <- "Higher in slower cell line"

plot_colors <- c(
  "Similar expression" = "gray70",
  "Higher in faster cell line" = "firebrick",
  "Higher in slower cell line" = "royalblue"
)

plot(
  DGE_df$fast_minus_slow,
  DGE_df$mean_expression,
  pch = 16,
  cex = 0.45,
  col = plot_colors[DGE_df$expression_group],
  xlab = paste("Expression difference:", fast_cell, "minus", slow_cell),
  ylab = "Mean log2(expression + 1)",
  main = "Volcano-style effect plot"
)

abline(v = 0, lwd = 2)
abline(v = c(-5, 5), lty = 2)

legend(
  "topright",
  legend = names(plot_colors),
  col = plot_colors,
  pch = 16,
  cex = 0.8,
  bty = "n"
)

The red points are genes much higher in the faster cell line, and the blue points are genes much higher in the slower cell line. Gray points are closer to equal expression. Points near the top are more highly expressed overall.

Question: Do the genes with the largest expression differences tend to be highly expressed, weakly expressed, or a mixture of both?

Your answer: Genes with the largest expression differences to be highly expressed in one of the cell lines, but this pattern does have some variation.

Important: This is not a statistical volcano plot because it does not contain p-values. A true volcano plot requires biological replicates for each group.

Genes higher in the faster cell line

DGE_ordered <- DGE_df[order(DGE_df$fast_minus_slow, decreasing = TRUE), ]

head(DGE_ordered, 15)

##              gene slow_expression fast_expression fast_minus_slow
## VIM           VIM      1.23878686       10.959234        9.720447
## LDHB         LDHB      0.28688115        9.076067        8.789186
## SERPINE1 SERPINE1      0.01435529        8.718704        8.704349
## MSN           MSN      0.01435529        8.709945        8.695590
## GPX1         GPX1      0.04264434        8.720347        8.677703
## CAV1         CAV1      0.00000000        8.457955        8.457955
## GSTP1       GSTP1      0.59454855        8.961392        8.366843
## FOSL1       FOSL1      0.68706069        8.750841        8.063780
## AXL           AXL      0.18903382        8.201781        8.012747
## F3             F3      0.69599381        8.651411        7.955417
## CST1         CST1      0.00000000        7.920532        7.920532
## PLAT         PLAT      0.09761080        7.883315        7.785705
## MMP14       MMP14      0.61353165        8.398744        7.785212
## AKR1B1     AKR1B1      0.08406426        7.820881        7.736817
## PLAU         PLAU      0.51601515        8.167368        7.651353
##          mean_expression           expression_group
## VIM             6.099010 Higher in faster cell line
## LDHB            4.681474 Higher in faster cell line
## SERPINE1        4.366530 Higher in faster cell line
## MSN             4.362150 Higher in faster cell line
## GPX1            4.381496 Higher in faster cell line
## CAV1            4.228978 Higher in faster cell line
## GSTP1           4.777970 Higher in faster cell line
## FOSL1           4.718951 Higher in faster cell line
## AXL             4.195407 Higher in faster cell line
## F3              4.673702 Higher in faster cell line
## CST1            3.960266 Higher in faster cell line
## PLAT            3.990463 Higher in faster cell line
## MMP14           4.506138 Higher in faster cell line
## AKR1B1          3.952473 Higher in faster cell line
## PLAU            4.341692 Higher in faster cell line

These are genes with higher expression in MDA-MB-231 than in T-47D.

Genes higher in the slower cell line

tail(DGE_ordered, 15)

##              gene slow_expression fast_expression fast_minus_slow
## ST14         ST14        7.778734      0.87184365       -6.906891
## AZGP1       AZGP1        8.022312      0.57531233       -7.447000
## KRT23       KRT23        7.504700      0.00000000       -7.504700
## FOXA1       FOXA1        7.829723      0.04264434       -7.787078
## OLFM1       OLFM1        7.868699      0.00000000       -7.868699
## RAB25       RAB25        8.054957      0.09761080       -7.957346
## AGR3         AGR3        8.373170      0.11103131       -8.262138
## IGFBP5     IGFBP5        8.417684      0.02856915       -8.389115
## CDH1         CDH1        8.760354      0.17632277       -8.584031
## STC2         STC2        9.032679      0.41142625       -8.621253
## CRABP2     CRABP2       11.248586      2.49057013       -8.758016
## SERPINA6 SERPINA6        9.490229      0.04264434       -9.447585
## CRIP1       CRIP1        9.928415      0.31034012       -9.618075
## PIP           PIP       11.931051      0.00000000      -11.931051
## MGP           MGP       12.594747      0.27500705      -12.319740
##          mean_expression           expression_group
## ST14            4.325289 Higher in slower cell line
## AZGP1           4.298812 Higher in slower cell line
## KRT23           3.752350 Higher in slower cell line
## FOXA1           3.936184 Higher in slower cell line
## OLFM1           3.934350 Higher in slower cell line
## RAB25           4.076284 Higher in slower cell line
## AGR3            4.242100 Higher in slower cell line
## IGFBP5          4.223126 Higher in slower cell line
## CDH1            4.468338 Higher in slower cell line
## STC2            4.722053 Higher in slower cell line
## CRABP2          6.869578 Higher in slower cell line
## SERPINA6        4.766437 Higher in slower cell line
## CRIP1           5.119377 Higher in slower cell line
## PIP             5.965526 Higher in slower cell line
## MGP             6.434877 Higher in slower cell line

These are genes with higher expression in T-47D than in MDA-MB-231.

Count genes in the tails

num_fast_tail <- sum(DGE_df$fast_minus_slow > 5)
num_slow_tail <- sum(DGE_df$fast_minus_slow < -5)

num_fast_tail

## [1] 134

num_slow_tail

## [1] 94

Summarize your data:

paste(
  num_fast_tail,
  "genes are much higher in the faster cell line and",
  num_slow_tail,
  "genes are much higher in the slower cell line using a cutoff of 5."
)

## [1] "134 genes are much higher in the faster cell line and 94 genes are much higher in the slower cell line using a cutoff of 5."

The cutoff of 5 is arbitrary but useful for exploration.

Biological marker genes

Some genes are already known to relate to cell state and cancer biology.

Here are a few examples:

VIM: vimentin, often associated with mesenchymal/invasive cell states
CDH1: E-cadherin, often associated with epithelial cell identity
EPCAM: epithelial cell adhesion molecule
KRT8, KRT18, KRT23: keratin genes
MKI67: proliferation marker

marker_genes <- c("VIM", "CDH1", "APOE", "LDHB", "KRT18", "KRT23", "MKI67", "MSN", "PRLR", "KRAS")

marker_genes <- marker_genes[marker_genes %in% DGE_df$gene]

DGE_df[match(marker_genes, DGE_df$gene), ]

##        gene slow_expression fast_expression fast_minus_slow mean_expression
## VIM     VIM      1.23878686      10.9592340      9.72044717        6.099010
## CDH1   CDH1      8.76035403       0.1763228     -8.58403126        4.468338
## APOE   APOE      8.35416121       1.8519988     -6.50216237        5.103080
## LDHB   LDHB      0.28688115       9.0760673      8.78918619        4.681474
## KRT18 KRT18     11.35680957       8.7070830     -2.64972658       10.031946
## KRT23 KRT23      7.50469983       0.0000000     -7.50469983        3.752350
## MKI67 MKI67      3.91552090       5.6452982      1.72977726        4.780410
## MSN     MSN      0.01435529       8.7099454      8.69559009        4.362150
## PRLR   PRLR      7.70494111       1.3785116     -6.32642949        4.541726
## KRAS   KRAS      3.61588707       3.6170633      0.00117627        3.616475
##                 expression_group
## VIM   Higher in faster cell line
## CDH1  Higher in slower cell line
## APOE  Higher in slower cell line
## LDHB  Higher in faster cell line
## KRT18         Similar expression
## KRT23 Higher in slower cell line
## MKI67         Similar expression
## MSN   Higher in faster cell line
## PRLR  Higher in slower cell line
## KRAS          Similar expression

Question: Which marker genes are higher in the faster cell line? Which are higher in the slower cell line?

Your answer: VIM, MSN, and KRT18 have higher expression in the faster cell line, whereas CDH1, APOE and again KRT18 have higher expression in the slower cell line.

Scatter plot: slow vs fast expression

A scatter plot helps us compare expression in the two cell lines.

Each point is one gene.

plot(
  DGE_df$slow_expression,
  DGE_df$fast_expression,
  pch = 16,
  cex = 0.35,
  xlab = paste("Expression in", slow_cell),
  ylab = paste("Expression in", fast_cell),
  main = "Gene expression comparison"
)

abline(0, 1, lwd = 2)

Genes near the diagonal have similar expression in both cell lines.

Genes far above the diagonal are higher in the faster cell line.

Genes far below the diagonal are higher in the slower cell line.

Highlight marker genes

plot(
  DGE_df$slow_expression,
  DGE_df$fast_expression,
  pch = 16,
  cex = 0.35,
  xlab = paste("Expression in", slow_cell),
  ylab = paste("Expression in", fast_cell),
  main = "Marker genes in expression space"
)

abline(0, 1, lwd = 2)

marker_rows <- match(marker_genes, DGE_df$gene)

points(
  DGE_df$slow_expression[marker_rows],
  DGE_df$fast_expression[marker_rows],
  pch = 19,
  cex = 1.2,
  col="red"
)

text(
  DGE_df$slow_expression[marker_rows],
  DGE_df$fast_expression[marker_rows],
  labels = marker_genes,
  pos = 4,
  cex = 0.8,
  col="red"
)

Question: Where is VIM on this plot? Where are keratin genes such as KRT8, KRT18, or KRT23?

Your answer: VIM is an outlier on this plot that does not exhibit a strong pattern of expression between both cell lines. It is expressed more in the faster cell lines. Keratin genes have the opposite effect of being expressed more in slower cell lines, following the data.

Make gene lists for enrichment analysis

We can write the top genes to files and use them in a gene set enrichment tool.

The top genes are higher in the faster cell line.

N <- 50

fast_genes <- DGE_ordered$gene[1:N]

write.table(
  fast_genes,
  "fast_cell_line_genes.csv",
  row.names = FALSE,
  col.names = FALSE,
  quote = FALSE
)

head(fast_genes)

## [1] "VIM"      "LDHB"     "SERPINE1" "MSN"      "GPX1"     "CAV1"

The bottom genes are higher in the slower cell line.

slow_genes <- tail(DGE_ordered$gene, N)

write.table(
  slow_genes,
  "slow_cell_line_genes.csv",
  row.names = FALSE,
  col.names = FALSE,
  quote = FALSE
)

head(slow_genes)

## [1] "SLC7A8"   "TGFB3"    "RERG"     "PRSS8"    "GREB1"    "KIAA1324"

You can copy these gene lists into a gene set enrichment tool such as (Gene Set AI)[https://idekerlab.ucsd.edu/gsai/] or another enrichment analysis server, e.g. (g:Profiler)[biit.cs.ut.ee/gprofiler/gost].

Question: What biological processes are enriched among the genes higher in the faster cell line?

Your answer: Pathways for cell motility and collagen production.

Question: What biological processes are enriched among the genes higher in the slower cell line?

Your answer: Pathways for keratinization and intracellular adhesion.

CHALLENGE 1: Compare all breast cancer conditions

The dataset contains additional measurements for both breast cancer cell lines on other substrates.

brca_speeds_df[, c("sample", "cellLine", "experimentalCondition", "average_value")]

##      sample   cellLine             experimentalCondition average_value
## 42 mRNA_R52      T-47D                             Glass     13.770019
## 43 mRNA_R56      T-47D           HyaluronicAcid Collagen     16.515478
## 44 mRNA_R55      T-47D        HyaluronicAcid Fibronectin     11.288006
## 45 mRNA_R54      T-47D    30 kPa polyacrylamide Collagen      7.475378
## 46 mRNA_R53      T-47D 30 kPa polyacrylamide Fibronectin      8.923515
## 47 mRNA_R51      T-47D    500 Pa polyacrylamide Collagen     10.608637
## 48 mRNA_R50      T-47D 500 Pa polyacrylamide Fibronectin     16.734859
## 56 mRNA_R59 MDA-MB-231                             Glass     22.419595
## 57 mRNA_R63 MDA-MB-231           HyaluronicAcid Collagen     36.196133
## 58 mRNA_R62 MDA-MB-231        HyaluronicAcid Fibronectin     55.548778
## 59 mRNA_R61 MDA-MB-231    30 kPa polyacrylamide Collagen     31.749973
## 60 mRNA_R60 MDA-MB-231 30 kPa polyacrylamide Fibronectin     33.378558
## 61 mRNA_R58 MDA-MB-231    500 Pa polyacrylamide Collagen     35.848219
## 62 mRNA_R57 MDA-MB-231 500 Pa polyacrylamide Fibronectin     34.222731

Make a barplot comparing all breast cancer speed experiments.

barplot(
  brca_speeds_df$average_value,
  names.arg = paste(brca_speeds_df$cellLine, brca_speeds_df$experimentalCondition, sep = "\n"),
  las = 2,
  cex.names = 0.55,
  ylab = "Average speed (microns per hour)",
  main = "Breast cancer cell speed across substrates"
)

Question: Is MDA-MB-231 always faster than T-47D, or does substrate matter?

Your answer: It is usually faster, but the changes in substrate do result in difference between the range of which the speeds vary between the cell lines.

CHALLENGE 2: TPM column sums

The expression values are TPM-like values. TPM stands for transcripts per million.

If a complete TPM matrix is used, the values in each sample column often sum to about 1,000,000.

expr_col_sums <- colSums(pson_expr_mat)

summary(expr_col_sums)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  922843  934158  940822  940788  949631  953272

hist(
  expr_col_sums,
  breaks = 30,
  main = "Column sums of expression matrix",
  xlab = "Column sum"
)

Question: Do the columns sum to 1,000,000? Why might they not?

Your answer: Not exactly, this may be due to processing the matrix which may not reflect the raw nature of the data, which is more likely to sum to a million.

Possible explanations include filtering, preprocessing, rounding, or use of a transformed/processed expression matrix rather than raw full TPM values.

OPTIONAL CHALLENGE 3: Try a different cancer type

Choose another cancer type from the dataset.

table(cell_speeds_df$diagnosis)

## 
##    Brain Cancer   Breast Cancer    Colon Cancer  Not Applicable Prostate Cancer 
##              14              14              13               7               7 
##     Skin Cancer 
##               7

Try repeating the same analysis for another pair of cell lines from one cancer type.

This is a more open-ended challenge. Colon cell lines are a nice comparison here:

Colon cancer cell lines

SW620 cells are considered more aggressive than SW480 cells. Both cell lines were derived from the same patient, but at different times and locations:

SW480 cells are from the primary tumor, and
SW620 cells are from a lymph node metastasis.

This makes them a valuable pair for studying colon cancer progression and metastasis in vitro.

We looked at the speed and data for breast cancer cell lines on the HyaluronicAcid Collagen substrate. Use the colon cancer data for the experimental condition HyaluronicAcid Fibronectin.

Final reflection

In this activity, you connected two kinds of data:

Cell behavior: how fast cells move
Gene expression: which genes are more or less active

The faster breast cancer cell line, MDA-MB-231, shows higher expression of genes such as VIM, while the slower T-47D cell line shows higher expression of several epithelial/keratin-associated genes.

This does not prove that any one gene causes the motility difference. But it gives us a biological hypothesis:

differences in cell state and cytoskeletal/epithelial markers may help explain differences in cancer cell movement.

That is how exploratory computational biology often works: we use data to generate careful, testable hypotheses.

DREAM-High: Breast Cancer Cell Lines, Cell Motility, and Gene Expression

Cancer Systems Biology Consortium DREAM-High Team

Sinchana Nirguna Joshi

Summer 2026