knitr::opts_chunk$set(collapse = TRUE, comment = "#>", warning = FALSE, message = FALSE, fig.align = 'center', results = 'asis', fig.show = 'hold', fig.width = 7, fig.height = 5)
library(survivALL) library(survival) library(survcomp) library(Biobase) library(magrittr) library(ggplot2) library(GGally) library(ggthemes) library(cowplot); theme_set(theme_grey())
Survival analysis typically stratifies a patient cohort by some measure and compares the relative survival between groups. In practice, for a continuous measure, the decision of where to draw this division relies on previous knowledge or, more problematically, it may simply be a researchers best guess - for example separating around the median. Below are outlined 3 alternative survival analysis approaches --- median, hypothesis-driven and our own data-driven approach as part of the
We begin with an expressionSet
nki, which details both gene expression as well as survival information (i.e. event and time-to-event data) for 337 invasive breast cancer patients of no specific subtype. Here the event we will be measuring is distant metastasis free survival or dmfs.
It is important to ensure that there are no
NA values in our survival information. In our example this has already been accounted for and the data is included as part of the
library(pander) library(magrittr) pData(nki_subset)[1:3, ] %>% pandoc.table(caption = "Survival Information")
We can see that our complete survival data includes the three critical components of our analysis, namely sample names, events (e.dmfs) and times-to-event (t.dmfs).
To determine a variables prognostic capacity we can apply the Kaplan-Meier estimator and plot. Here we will split our cohort using the expression of ERBB2 a gene whose increased expression is known to be associated with poor prognosis in a mixed population of invasive breast cancers.
We will stratify the cohort into high and low ERBB2 expression using median ERBB2 expression...
#ERBB2 expression vector erbb2_xpr <- exprs(nki_subset)["NM_004448",] #convert to binary classifier erbb2_med <- ifelse(erbb2_xpr >= median(erbb2_xpr), "high", "low")
...and produce our Kaplan-Meier plot and statistics
srv_obj <- survival::Surv(nki_subset$t.dmfs, nki_subset$e.dmfs) median_fit <- survival::survfit(srv_obj ~ erbb2_med) p_med <- GGally::ggsurv(median_fit, surv.col = c("#525252", "#bdbdbd")) + ylim(0.4, 1) + labs(title = "Figure S1", subtitle = "Median approach", x = "Years") + theme_pander() + theme(axis.line = element_line(size = 0.1), legend.position = 'bottom') p_med <- ggdraw(add_sub(p_med, "p = 0.9", vpadding=grid::unit(0, "lines"), y = 17, x = 0.55, hjust = 0, size = 10)) p_med
broom::tidy(survival::coxph(srv_obj ~ erbb2_med)) %>% pandoc.table()
Surprisingly and despite ERBB2 known to be highly prognostic in invasive breast cancer, there is no association between ERBB2 expression and prognosis.
However, knowing that ERBB2 is overexpressed in ~25% of invasive breast cancer cases, we can modify our cohort stratification to reflect this.
erbb2_hypothesis <- ifelse(erbb2_xpr >= quantile(erbb2_xpr, probs = 0.75), "high", "low") hypothesis_fit <- survival::survfit(srv_obj ~ erbb2_hypothesis) p_hyp <- GGally::ggsurv(hypothesis_fit, surv.col = c("#525252", "#bdbdbd")) + ylim(0.4, 1) + labs(title = "Figure S2", subtitle = "Hypothesis-driven approach", x = "Years") + theme_pander() + theme(axis.line = element_line(size = 0.1), legend.position = 'none') p_hyp <- ggdraw(add_sub(p_hyp, "p = 0.043", vpadding=grid::unit(0, "lines"), y = 17, x = 0.55, hjust = 0, size = 10)) p_hyp
broom::tidy(survival::coxph(srv_obj ~ erbb2_hypothesis)) %>% pandoc.table()
And consequently demonstrate a significant association between prognosis and ERBB2 expression.
By comparing multiple patient stratifications we have successful revealed ERBB2's relationship with prognosis based on a reasoned estimate of where we expect optimal stratification to occur. It may still be, however, that our cohort's composition in terms of ERBB2 expression does not exactly mirror population levels - i.e. that the best point-of-separation is not found with a 75-25% split. This will become even more problematic for a novel biomarker for which population levels may be unknown. A solution to this unknown quality is, instead of estimating where we anticipate our point-of-separation to occur, we calculate survival statistics for all separations possible.
For our 319 patients and having calculated our median and hypothesis divisions, there remain an additional 316 points at which we could otherwise stratify our cohort, any of which may demonstrate a more significant association with prognosis. By additionally calculating survival statistics for each of these possible point-of-separation we guarantee that we reveal the true optimal stratification.
We can visually inspect the best point of separation by using the
plotALL(measure = erbb2_xpr, srv = pData(nki_subset), time = "t.dmfs", event = "e.dmfs") + labs(title = "Figure S3", subtitle = "survivALL plot output") + theme(plot.title = element_text(hjust = 0))
We can now see that the best point at which we could separate our cohort based on ERBB2 expression is closer to a 85-15% split, rather than 75-25%, indicating a sampling bias towards HER2 negative patients in the NKI cohort.
We can further investigate the most significant separation using
survivALL(), which follows the same procedure as
plotALL() but returns a dataframe of calculations rather than plotting.
Using this we can select the point of most significant separation...
srvall <- survivALL(measure = erbb2_xpr, srv = pData(nki_subset), time = "t.dmfs", event = "e.dmfs", measure_name = "ERBB2") srvall[which.min(srvall$p),] %>% pandoc.table()
...and re-draw our Kaplan-Meier
erbb2_data <- ifelse(srvall$clsf == 0, "low", "high") srv_obj <- survival::Surv(as.numeric(srvall$event_time), srvall$event) data_fit <- survival::survfit(srv_obj ~ erbb2_data) p_data <- GGally::ggsurv(data_fit, surv.col = c("#525252", "#bdbdbd")) + ylim(0.4, 1) + labs(title = "Figure S4", subtitle = "Data-drive approach", x = "Years") + theme_pander() + theme(axis.line = element_line(size = 0.1), legend.position = 'none') p_data <- ggdraw(add_sub(p_data, "p = 0.001", vpadding=grid::unit(0, "lines"), y = 17, x = 0.55, hjust = 0, size = 10)) p_data
broom::tidy(survival::coxph(srv_obj ~ erbb2_data)) %>% pandoc.table()
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.