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 `survivALL`

package.

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 `survivALL`

package.

```
data(nki_subset)
```

pData(nki_subset)[1:3, ]

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()`

function.

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

Embedding an R snippet on your website

Add the following code to your website.

For more information on customizing the embed code, read Embedding Snippets.