Home About Blog pRojects

Determine optimal cutpoints for numerical variables in survival plots

Nov 21, 2016 • Marcin Kosiński

The often demand in the biostatistical research is to group patients depending on explanatory variables that are continuous. In some cases the requirement is to test overall survival of the subjects that suffer on a mutation in specific gene and have high expression (over expression) in other given gene. To visualize differences in the Kaplan-Meier estimates of survival curves between groups, first the discretization of continuous variable is performed. Problems caused by categorization of continuous variables are known and widely spread (Harrel, 2015), but in this case there appear a simplification requirement for the discretization. In this post I present the maxstat(maximally selected rank statistics) statistic to determine the optimal cutpoint for continuous variables, which was provided it in the survminer package by Alboukadel Kassambara kassambara.

In this post I will use data from TCGA study, that are provided in the RTCGA package Star and survminer Star package to determine the optimal cutpoint for continuous variable.

Data preparation

I wrote about TCGA datasets and their preprocessing in my earlier posts: RTCGA factory of R packages - Quick Guide and BioC 2016 Conference Overview and Few Ways of Downloading TCGA Data. If your are not familiar with RTCGA family of data packages, you can visit the RTCGA website. Below I join survival information with ABCD4|5826 gene expression for patients suffering from BRCA (breast cancer) and HNSC (head and neck cancer). It can be done due to bcr_patient_barcode column which identifies each patient.

library(tidyverse) # pipes (%>%) and dplyr data munging
# installation
# source("https://bioconductor.org/biocLite.R")
# biocLite("RTCGA.clinical")
# biocLite("RTCGA.rnaseq")
library(RTCGA.clinical) # survival times
library(RTCGA.rnaseq) # genes' expression


survivalTCGA(BRCA.clinical, 
             HNSC.clinical) -> BRCA_HNSC.surv

expressionsTCGA(
   BRCA.rnaseq, HNSC.rnaseq,
   extract.cols = c("ABCD4|5826")) %>%
   rename(cohort = dataset,
          ABCD4 = `ABCD4|5826`) %>%
   filter(substr(bcr_patient_barcode, 14, 15) == "01") %>% 
   # only cancer samples
   mutate(bcr_patient_barcode = 
             substr(bcr_patient_barcode, 1, 12)) -> BRCA_HNSC.rnaseq

head(BRCA_HNSC.surv)
  times bcr_patient_barcode patient.vital_status
1  3767        TCGA-3C-AAAU                    0
2  3801        TCGA-3C-AALI                    0
3  1228        TCGA-3C-AALJ                    0
4  1217        TCGA-3C-AALK                    0
5   158        TCGA-4H-AAAK                    0
6  1477        TCGA-5L-AAT0                    0
head(BRCA_HNSC.rnaseq)
# A tibble: 6 × 3
  bcr_patient_barcode      cohort    ABCD4
                <chr>       <chr>    <dbl>
1        TCGA-3C-AAAU BRCA.rnaseq 322.2560
2        TCGA-3C-AALI BRCA.rnaseq 486.6775
3        TCGA-3C-AALJ BRCA.rnaseq 308.2502
4        TCGA-3C-AALK BRCA.rnaseq 589.1601
5        TCGA-4H-AAAK BRCA.rnaseq 615.7447
6        TCGA-5L-AAT0 BRCA.rnaseq 896.3191

Joining survival times and ABCD4 gene’ expression.

BRCA_HNSC.surv %>%
   left_join(BRCA_HNSC.rnaseq,
             by = "bcr_patient_barcode") ->
   BRCA_HNSC.surv_rnaseq

head(BRCA_HNSC.surv_rnaseq)
  times bcr_patient_barcode patient.vital_status      cohort    ABCD4
1  3767        TCGA-3C-AAAU                    0 BRCA.rnaseq 322.2560
2  3801        TCGA-3C-AALI                    0 BRCA.rnaseq 486.6775
3  1228        TCGA-3C-AALJ                    0 BRCA.rnaseq 308.2502
4  1217        TCGA-3C-AALK                    0 BRCA.rnaseq 589.1601
5   158        TCGA-4H-AAAK                    0 BRCA.rnaseq 615.7447
6  1477        TCGA-5L-AAT0                    0 BRCA.rnaseq 896.3191

13 patients have clinical info but they do not have expression information so I remove them from the analysis.

table(BRCA_HNSC.surv_rnaseq$cohort, useNA = "always")
BRCA.rnaseq HNSC.rnaseq        <NA> 
       1093         519          13 
BRCA_HNSC.surv_rnaseq <- BRCA_HNSC.surv_rnaseq %>%
   filter(!is.na(cohort))

The complete data used for further analysis is printed below

dim(BRCA_HNSC.surv_rnaseq)
[1] 1612    5
head(BRCA_HNSC.surv_rnaseq)
  times bcr_patient_barcode patient.vital_status      cohort    ABCD4
1  3767        TCGA-3C-AAAU                    0 BRCA.rnaseq 322.2560
2  3801        TCGA-3C-AALI                    0 BRCA.rnaseq 486.6775
3  1228        TCGA-3C-AALJ                    0 BRCA.rnaseq 308.2502
4  1217        TCGA-3C-AALK                    0 BRCA.rnaseq 589.1601
5   158        TCGA-4H-AAAK                    0 BRCA.rnaseq 615.7447
6  1477        TCGA-5L-AAT0                    0 BRCA.rnaseq 896.3191

maxstat - maximally selected rank statistics

kassambara prepared a functionality that uses the maxstat(maximally selected rank statistics) statistic to determine the optimal cutpoint for continuous variables and provided it in the survminer package. The development process is described on the survminer issues track

#devtools::install_github('kassambara/survminer')
library(survminer)

Articles in which the maxstat statistics was used:

Explanation of the selection process for the cutpoint is described in this vignette, chapter 2

Determining the optimal cutpoint for ABCD4 gene’s expression

BRCA_HNSC.surv_rnaseq.cut <- surv_cutpoint(
   BRCA_HNSC.surv_rnaseq,
   time = "times",
   event = "patient.vital_status",
   variables = c("ABCD4", "cohort")
)
summary(BRCA_HNSC.surv_rnaseq.cut)
      cutpoint statistic
ABCD4 468.6417  1.944138

Plot the cutpoint for gene ABCD4

# palette = "npg" (nature publishing group), see ?ggpubr::ggpar
plot(BRCA_HNSC.surv_rnaseq.cut, "ABCD4", palette = "npg")
$ABCD4

plot of chunk unnamed-chunk-8

Categorize ABCD4 variable

BRCA_HNSC.surv_rnaseq.cat <- surv_categorize(BRCA_HNSC.surv_rnaseq.cut)

Fit and visualize Kaplan-Meier estimates of survival curves

Below I divided patients into 4 groups, denoting the membership to cancer type and patient’s ABCD4 gene’s expression level.

RTCGA way

RTCGA::kmTCGA(
   BRCA_HNSC.surv_rnaseq.cat, 
   explanatory.names = c("ABCD4", "cohort"),
   pval = TRUE,
   conf.int = TRUE,
   xlim = c(0,3000)
)

plot of chunk unnamed-chunk-10

survminer way

library(survival)
fit <- survfit(Surv(times, patient.vital_status) ~ ABCD4 + cohort,
               data = BRCA_HNSC.surv_rnaseq.cat)
ggsurvplot(
   fit,                     # survfit object with calculated statistics.
   risk.table = TRUE,       # show risk table.
   pval = TRUE,             # show p-value of log-rank test.
   conf.int = TRUE,         # show confidence intervals for 
                            # point estimaes of survival curves.
   xlim = c(0,3000),        # present narrower X axis, but not affect
                            # survival estimates.
   break.time.by = 1000,    # break X axis in time intervals by 500.
   ggtheme = theme_RTCGA(), # customize plot and risk table with a theme.
 risk.table.y.text.col = T, # colour risk table text annotations.
  risk.table.y.text = FALSE # show bars instead of names in text annotations
                            # in legend of risk table
)

plot of chunk unnamed-chunk-11