About this activity

This module will introduce you to loading and manipulating the data from the breast cancer METABRIC dataset, visualizing and working with gene expression measurements, and building predictive models based on the expression of many different genes and the associated clinical data.

The METABRIC dataset contains more patient survivial information than we had with the TCGA data.

Loading & inspecting data

load("/shared/dreamhigh/data/metabric_disc_clin_df.RData")
load("/shared/dreamhigh/data/metabric_disc_expr_mat.RData")
ls()
## [1] "metabric_disc_clin_df"  "metabric_disc_expr_mat"

Check out the clinical data frame. Many of the features are similar to those in the TCGA breast cancer clinical data frame.

# What's in the clinical data frame?

View(metabric_disc_clin_df)

Note that the expression matrix has already been log-transformed. How many rows and columns does metabric_disc_expr_mat have?

Reminder: The dimensions of a matrix can be found with the dim() function.

# How many rows and columns in the expression matrix?

 dim(metabric_disc_expr_mat)
## [1] 15832   997

What’s in the expression matrix?

Reminder: Indexing can be used to get specific rows and some columns of a matrix.

# The expression matrix is already log-transformed

metabric_disc_expr_mat[1:5 ,1:5]
##         MB_0362 MB_0346 MB_0386 MB_0574 MB_0185
## KIR2DL4   5.715   5.894   5.510   5.515   6.239
## NBPF10    7.954   7.872   8.684   8.535   8.407
## TSPAN6    7.579   7.819   7.366   8.110   5.993
## TNMD      6.211   5.557   8.108   5.789   5.484
## DPM1      9.150   9.794   9.229   9.672  10.159

Each row of the clinical data frame is one patient, and each column of the expression data frame is for one patient.

Check whether the patient IDs in the expression and clinical data frames are in the same order:

# Are the patient IDs the same and in the same order?

identical(metabric_disc_clin_df$metabric_id,colnames(metabric_disc_expr_mat))
## [1] TRUE

Now we are ready to do some actual analysis!

Exploring AURKA expression & survival

We happen to know that high AURKA expression in a tumor is bad for the patient’s outcome.

Question: What exactly does AURKA do in the cell cycle of healthy cells?

Your answer: It helps healthy cells divide by supporting centrosomes and spindle maturation.

Since AURKA (Aurora kinase A) encodes for a protein that is involved in regulation of the cell cycle, various studies have associated it with tumour progression and metastasis.

# Row number for AURKA in the expression data

aurka_row <- which(rownames(metabric_disc_expr_mat) == "AURKA")
aurka_row
## [1] 1597

Make an object that contains the expression of AURKA across patients:

# as.vector turns a row into a one-dimensional array
aurka <- metabric_disc_expr_mat[aurka_row, ]
colnames(aurka)<-NULL
aurka<-unlist(aurka)
head(aurka)
## [1] 7.918 9.626 7.689 7.738 8.612 7.083

Patient survival status

Since high AURKA expression is associated with tumor progression, we are interested in clinical features for which AURKA expression may be predictive.

The feature last_follow_up_status tells whether patients being followed in this study were alive or had passed at that time.

Remember we have used the function table to get summaries of TCGA clinical features.

# What status can patients have at the last time they were checked?

table(metabric_disc_clin_df$last_follow_up_status)
## 
##      a      d d-d.s. d-o.c. 
##    548     46    259    126

What do each of the categories mean?

At time of last follow up, a patient can have a status:

259 patients were specifically indicated to have died from breast cancer (d-d.s.). There are 46 patients labelled as d. We think they died of cancer, but the cause for these patients could be unknown.

Even for the patients who died of unknown cause d-o.c., the cause may have been indirectly caused by having cancer, so we usually want to be careful about how we treat these cases.

Censored records

Another perspective is to consider only those patients who died of disease and examine how long since diagnosis that death occurred.

We have to exclude the rows and columns for patients who are alive and those who died of other causes. These cases are known as censored events.

Rather than filtering out patients with last_follow_up_status equal to a or d-o.c., we can use the censored column in the data frame.

# How many patients are alive or have died?

table(metabric_disc_clin_df$censored)
## 
## FALSE  TRUE 
##   323   674

Question Compare these number to the numbers for a, d, d-d.s. and d-o.c. What do TRUE and FALSE correspond to?

Your answer: False would be diead from disease, while true would be alive or dead from different causes

Why do we want to remove patients with censored equal to TRUE?

Your answer: Because these patients died of causes not related to the disease or are alive, which is not relevant to this property.

# Remove the censored patients from the data
# The function which() gives the row number for censored = TRUE

drop_pats <- which(metabric_disc_clin_df$censored)

# Create clincial dataframe and expression matrix for patients who have died
dead_clin <- metabric_disc_clin_df[-drop_pats, ]
dead_expr <- metabric_disc_expr_mat[, -drop_pats]

dim(dead_clin)
## [1] 323  32
dim(dead_expr)
## [1] 15832   323

To test whether the code did what you expected, check out the new clinical data frame:

View(dead_clin)
# What is the status of patients in dead_clin?

table(dead_clin$last_follow_up_status)
## 
##      d d-d.s. 
##     46    259

Question Is this what you expect?

Your answer: Yes, as this matches the removing of the censored values.

The time-to-death (the feature T) is measured in days since cancer diagnosis.

Modeling survival and AURKA expression

How is time-to-death related to the expression level of AURKA?

We’ll create a plot that includes a trend line representing the linear model from lm().

# make sure to remove censored patients from the `aurka` vector
aurka_dead = aurka[-drop_pats]


# model time to death as a function of AURKA expression
aurka_model <- lm(dead_clin$T ~ aurka_dead)

# create a scatterplot to check your result
plot(dead_clin$T ~ aurka_dead,
     main = "Time to death as a function of AURKA expression",
     ylab = "Time to death, days",
     xlab = "AURKA expression")

abline(aurka_model, col = "red")

Time-to-death is inversely correlation with AURKA expression.

As we saw in the mtcars_linear_regression.Rmd activity, there are several things in the output of lm() that are interesting. We’ll store the summary information as an object and extract the values we want.

# Store the summary of the model you already made
aurka_model_summary <- summary(aurka_model)

coef_aurka_model = coefficients(aurka_model_summary)

# Get slope and p-value
slope = coef_aurka_model[2, 1]
pval = coef_aurka_model[2, 4]

r_squared <- aurka_model_summary$r.squared

slope
## [1] -711.7943
r_squared
## [1] 0.1723047
pval
## [1] 8.426028e-15

The time to death decreases as AURKA expression increases.

Add clinical variables to your model

In the code chunk above, you should have found that AURKA expression explains 17% of the variance (r.squares) in time to death (dead_clin$T).

Let’s include more factors in the model.

Continuous variables

AURKA expression in aurka_dead is a continuous variable. You can include other continuous variables in your model by adding them. For example:

Number of positive lymph nodes

The number of lymph nodes found to be positive for cancer cells at the time of diagnosis, lymph_nodes_positive, can be a strong predictor of time until the person dies. We would probably want to allow for variables like this in a survival model.

# Add lymph_nodes_positive as a feature to the model

summary(lm(dead_clin$T ~ dead_clin$lymph_nodes_positive + aurka_dead))
## 
## Call:
## lm(formula = dead_clin$T ~ dead_clin$lymph_nodes_positive + aurka_dead)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2400.1  -857.5  -179.4   679.9  4308.6 
## 
## Coefficients:
##                                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                     7599.93     676.37  11.236  < 2e-16 ***
## dead_clin$lymph_nodes_positive   -62.72      14.07  -4.459 1.14e-05 ***
## aurka_dead                      -678.83      85.19  -7.968 2.90e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1214 on 318 degrees of freedom
##   (2 observations deleted due to missingness)
## Multiple R-squared:  0.221,  Adjusted R-squared:  0.2161 
## F-statistic: 45.11 on 2 and 318 DF,  p-value: < 2.2e-16

Question Does time-to-death increase or decrease with lymph_nodes_positive? How much variance is explained with the new model?

Your answer: Time to death increases with the lymph nodes being positive, which may explain some variance, as the model still may be showing a relationship that is based on many more factors.

NPI

Try the same with NPI. The Nottingham Prognostic Index (NPI) is a tool used to predict the prognosis of breast cancer patients after surgery. It combines three key factors: tumor size, lymph node involvement, and tumor grade. A higher score corresponds to poorer prognosis.

summary(lm(dead_clin$T ~ dead_clin$NPI + aurka_dead))
## 
## Call:
## lm(formula = dead_clin$T ~ dead_clin$NPI + aurka_dead)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2830.8  -782.4  -260.5   605.3  4298.4 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    7971.96     677.71  11.763  < 2e-16 ***
## dead_clin$NPI  -294.20      63.74  -4.616 5.69e-06 ***
## aurka_dead     -579.69      89.40  -6.484 3.40e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1212 on 318 degrees of freedom
##   (2 observations deleted due to missingness)
## Multiple R-squared:  0.2243, Adjusted R-squared:  0.2194 
## F-statistic: 45.97 on 2 and 318 DF,  p-value: < 2.2e-16

Question Are your results consistent with the definition of the NPI score? Why?

Your answer: No, because there is a large amount of error as seen above.

Modeling survival for all genes

Putting all the bits together, we can write a function, gene_lm() that takes the time-to-death (y) and fits it to some variable (x), such as the expression of AURKA, and extracts the slope (direction of the association) and statistical significance (pval).

More generally, yes, we can write our own functions in R!

# Do a linear fit of y ~ x and return the slope and p-value
gene_lm <- function(y, x) {
    # Keep only paired, finite observations.
    ok <- is.finite(y) & is.finite(x)

    # A slope cannot be estimated without at least two complete observations
    # and two distinct x values.
    if (sum(ok) < 2 || length(unique(x[ok])) < 2) {
        return(c(pval = NA_real_, slope = NA_real_))
    }

    result <- tryCatch({
        gene_model_summary <- summary(lm(y[ok] ~ x[ok]))
        coef_gene_model <- coefficients(gene_model_summary)
        c(pval = coef_gene_model[2, "Pr(>|t|)"],
          slope = coef_gene_model[2, "Estimate"])
    }, error = function(e) {
        c(pval = NA_real_, slope = NA_real_)
    })

    return(result)
}

Check that you get the same result as before.

# Run your new function

gene_lm(dead_clin$T, aurka_dead)
##          pval         slope 
##  8.426028e-15 -7.117943e+02

We could check every gene in the expression matrix to see if it is significantly correlated with time-to-death of the METABRIC patients. Fortunately with R, we can easily loop over all the genes. We’ll create an empty matrix of results for storing the p-value and slope for each gene. Altogether, 15,834 genes were measured and included in the expression matrix.

# Try the prediction first for 5000 genes

ngenes <- min(5000, nrow(dead_expr))
lm_results <- matrix(nrow = ngenes, ncol = 2, data = NA_real_)
colnames(lm_results) <- c("pval", "slope")
rownames(lm_results) <- rownames(dead_expr)[1:ngenes]
for (i in seq_len(ngenes)) {
    lm_results[i, ] <- gene_lm(dead_clin$T, as.numeric(dead_expr[i, ]))
}

And we check on AURKA again.

lm_results[rownames(lm_results) == "AURKA", ]
##          pval         slope 
##  8.426028e-15 -7.117943e+02

That’s it. We now have the regression of every gene against time-to-death. There are some other things we need to do in practice, like correcting the p-values for the number of tests we did. If you check 1000 genes, then for p = 0.01 you would expect roughly 10 genes to be that significant by chance (0.01*1000).

From the results table lm_results, we can order by p-value to see which genes are most significantly associated with death (either favourably or unfavourably).

lm_results = lm_results[order(lm_results[, 1]), ]

head(lm_results,20)
##                  pval      slope
## AVP      5.637819e-18  1678.9654
## CENPA    1.130745e-16  -855.0454
## ANLN     7.649051e-16  -995.2773
## RABEP1   2.573654e-15   627.6337
## AURKA    8.426028e-15  -711.7943
## FAM83D   4.566357e-14  -733.2090
## MCM10    7.459059e-14  -751.7670
## TRIP13   9.654274e-14  -710.3487
## STMN1    4.681205e-13  -807.2726
## SEC14L2  7.382927e-13   501.6181
## KIAA0141 8.697885e-13  1449.0480
## MXD4     9.772895e-13   934.6374
## KARS     1.139240e-12 -1169.7592
## ESR1     1.412970e-12   236.0715
## TPX2     1.464180e-12  -732.2328
## MZF1     1.700396e-12   822.2763
## CYBRD1   2.220941e-12   499.1818
## CDC20    2.256791e-12  -485.7147
## HPN      2.285301e-12   482.4860
## CPEB3    2.874603e-12  1068.1870

AURKA is near the top. The ANLN gene codes for the protein Anillin, an actin-binding protein required for cytokinesis, is a prognostic marker panel in breast cancer..

Several of the genes seem to be favorable for survival, such as MXD4 and RABEP1. You can learn more about their functions at the UniProt Knowledgebase.

Prediction

Let’s use your model for prediction. You know already that your model accounts for only ~20% of the variance in time-to-death, so we don’t expect stellar performance.

If you look at dead_clin, you’ll see there are a couple of cases where dead_clin$T has a value of NA:

which(is.na(dead_clin$T))
## [1]  87 316

Predictions, which are numeric values for time-to-death, will be made for all patients with a recorded time-to-death. To compare predicted values with observed values, remove patients whose T value is missing or non-finite. Use a logical index rather than hard-coded row numbers so the clinical data and expression matrix remain aligned if the data change.

keep_patients <- is.finite(dead_clin$T)
new_dead_clin <- dead_clin[keep_patients, , drop = FALSE]
new_dead_expr <- dead_expr[, keep_patients, drop = FALSE]

# Confirm that patients are still aligned in both objects.
stopifnot(
    nrow(new_dead_clin) == ncol(new_dead_expr),
    identical(new_dead_clin$metabric_id, colnames(new_dead_expr))
)

Let’s create a model based on expression of a gene that may be important for survival. We’ve been focussing on AURKA, but you can choose another- perhaps one of the most significant genes from lm_result when you modeled survival for all genes.

For prediction, it is easiest if we add the column for gene expression to the clinical data frame.

# Pick the gene to use in your model here.
# AURKA is known to be present in this activity's expression matrix.
my_model_gene <- "AURKA"

if (!my_model_gene %in% rownames(new_dead_expr)) {
    stop("The selected gene is not present in the expression matrix: ",
         my_model_gene)
}

# add gene expression to the training clinical table above
new_dead_clin$gene_exp <- as.numeric(new_dead_expr[my_model_gene, ])

if (sum(is.finite(new_dead_clin$gene_exp)) < 2 ||
    length(unique(new_dead_clin$gene_exp[is.finite(new_dead_clin$gene_exp)])) < 2) {
    stop("The selected gene does not have enough usable variation for linear regression: ",
         my_model_gene)
}

Now train your model for survival as a function of gene expression…

gene_model <- lm(T ~ gene_exp, data = new_dead_clin, na.action = na.exclude)
summary(gene_model)
## 
## Call:
## lm(formula = T ~ gene_exp, data = new_dead_clin, na.action = na.exclude)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3003.3  -844.4  -256.2   773.7  4508.6 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  7678.72     695.86  11.035  < 2e-16 ***
## gene_exp     -711.79      87.35  -8.149 8.43e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1250 on 319 degrees of freedom
## Multiple R-squared:  0.1723, Adjusted R-squared:  0.1697 
## F-statistic: 66.41 on 1 and 319 DF,  p-value: 8.426e-15

Look at the coefficient named gene_exp. Based on this coefficient, is the relationship between the expression of your gene and survival time positive or negative?

You see that your gene has some explanatory value (Multiple R-squared), let’s see if it has predictive value.

The function predict.lm returns predicted values for the dependent feature, in our case time-to-death T.

# The object my_prediction will hold your predicted values for T

my_prediction <- data.frame(
    metabric_id = new_dead_clin$metabric_id,
    T = predict(gene_model, newdata = new_dead_clin)
)
head(my_prediction)
##    metabric_id         T
## 1      MB_0362 2042.7343
## 2      MB_0346  826.9897
## 5      MB_0185 1548.7491
## 11     MB_0189 2528.1780
## 19     MB_0223 2548.1082
## 26     MB_0906 2123.1670

To get a sense of how well you model did, let’s compare the predicted values with the real values:

# correlation
cor(my_prediction$T, new_dead_clin$T, use = "complete.obs")
## [1] 0.415096

Question Given the correlation you just calculated, what do you think of your model?

Your answer: My model is not sophisticated enough, as the R value does not show strong positive correlation, but a moderate value. So, my model needs more work to reach a higher R value, such as a number near 0.7/0.8, to better draw conclusions from.

Wrap-up

Congratulations on accomplishing another activity!! You’ve done a tremendous amount of work this summmer and learned a wide range of data analysis methods with many different data types. Great work.

Don’t forget to knit and publish to Rpubs!