Rationale

In survivALL: Continuous Biomarker Assessment by Exhaustive Survival Analysis

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)

Rationale

 

Libraries

library(survivALL)
library(survival)
library(survcomp)
library(Biobase)
library(magrittr)
library(ggplot2)
library(GGally)
library(ggthemes)
library(cowplot); theme_set(theme_grey())

 

Rationale

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.

 

Median dichotomisation

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.

 

Hypothesis-driven dichotomisation

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.

 

Data-driven dichotomisation

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

 

 

survivALL documentation built on May 1, 2019, 11:02 p.m.