Contents

suppressPackageStartupMessages({
  library(tidyverse)
  library(survival)
  library(dplyr)
  library(survminer)
  library(Biobase)
  library(ggsci)
  library(ggbeeswarm)
  library(TOP)
  library(curatedOvarianData)
})

theme_set(theme_bw())

1 Installation

# Install the package from Bioconductor
if (!requireNamespace("BiocManager", quietly = TRUE)) {
  install.packages("BiocManager")
}
BiocManager::install("TOP")

2 Overview

The TOP R package provides a transfer learning approach for building predictive models across multiple omics datasets. With the increasing availability of omics data, there is a growing need for methods that can effectively integrate and analyze data from multiple sources. However, merging datasets can be challenging due to batch effects and other sources of variation.

TOP uses transfer learning strategies to build predictive models that can be applied across multiple datasets without the need for extensive batch correction methods. Specifically, TOP employs a lasso regression model to identify important features and construct a transferable predictive model. By leveraging information from multiple datasets, TOP can improve predictive accuracy and identify common biomarkers across different omics datasets.

The package provides several functions for building and evaluating transfer learning models, including options for visualizing results. Additionally, the package includes sample datasets and detailed documentation to help users get started with the package and understand the underlying methods.

Overall, TOP offers a flexible and powerful approach for integrating and analyzing omics data from multiple sources, making it a valuable tool for researchers in a variety of fields.

3 Loading example data

The example data used in this vignette is the curatedOvarianData. Described in the paper from Benjamin Frederick Ganzfried et al (2013) 

data("GSE12418_eset", package = "curatedOvarianData")
data("GSE30009_eset", package = "curatedOvarianData")

4 Building a survival model.

The transferable omics prediction framework is also able to build survival models. In this section, we will demonstrate how to use TOP to build a survival model using example ovarian cancer datasets. First we utilise nine datasets that were preparred by Ganzfield et al. 

data("GSE32062.GPL6480_eset")
japan_a <- GSE32062.GPL6480_eset

data("GSE9891_eset")
tothill <- GSE9891_eset

data("GSE32063_eset")
japan_b <- GSE32063_eset

data("TCGA.RNASeqV2_eset")
selection <- TCGA.RNASeqV2_eset$tumorstage %in% c(3, 4) & TCGA.RNASeqV2_eset$site_of_tumor_first_recurrence == "metastasis"
selection[is.na(selection)] <- FALSE
tcgarnaseq <- TCGA.RNASeqV2_eset[, selection]

data("E.MTAB.386_eset")
bentink <- E.MTAB.386_eset

data("GSE13876_eset")
crijns <- GSE13876_eset
crijns <- crijns[, crijns$grade %in% c(3, 4)]

data("GSE18520_eset")
mok <- GSE18520_eset

data("GSE17260_eset")
yoshihara2010 <- GSE17260_eset

data("GSE26712_eset")
bonome <- GSE26712_eset

list_ovarian_eset <- lst(
  japan_a, tothill, japan_b,
  tcgarnaseq, bonome, mok, yoshihara2010,
  bentink, crijns
)

list_ovarian_eset %>%
  sapply(dim)
#>          japan_a tothill japan_b tcgarnaseq bonome   mok yoshihara2010 bentink
#> Features   20106   19816   20106      20502  13104 19816         20106   10357
#> Samples      260     285      40         51    195    63           110     129
#>          crijns
#> Features  20577
#> Samples      85

5 Common genes between datasets

In order to apply the TOP framework, it is important that the input matrices have identical feature names, such as gene names, across all datasets. In this example, we will identify the common genes present in all the datasets, as these will be the features used for transfer learning.

raw_gene_list <- purrr::map(list_ovarian_eset, rownames)
common_genes <- Reduce(f = intersect, x = raw_gene_list)
length(common_genes)
#> [1] 7517

6 Survival samples

Next, we will prepare the survival data from each of the ovarian cancer datasets. The survival data includes the survival time and the event status (i.e., whether the event of interest, such as death, has occurred)

ov_pdata <- purrr::map(list_ovarian_eset, pData)
list_pdata <- list_ovarian_eset %>%
  purrr::map(pData) %>%
  purrr::map(tibble::rownames_to_column, var = "sample_id")

ov_surv_raw <- purrr::map(
  .x = list_pdata,
  .f = ~ data.frame(
    sample_id = .x$sample_id,
    time = .x$days_to_death %>% as.integer(),
    dead = ifelse(.x$vital_status == "deceased", 1, 0)
  ) %>%
    na.omit() %>%
    dplyr::filter(
      time > 0,
      !is.nan(time),
      !is.nan(dead)
    )
)
ov_surv_raw %>% sapply(nrow)
#>       japan_a       tothill       japan_b    tcgarnaseq        bonome 
#>           260           276            40            51           185 
#>           mok yoshihara2010       bentink        crijns 
#>            53           110           129            85
ov_surv_y <- ov_surv_raw %>%
  purrr::map(~ .x %>%
    dplyr::select(-sample_id)) %>%
  purrr::map(~ Surv(time = .x$time, event = .x$dead))

7 Preparing data for modelling.

In this section, we will prepare the gene expression data and survival data for visualization. We will subset the gene expression data to include only the common genes and samples with survival information. Then, we will create a combined survival data table, with a data_source column to identify the origin of each sample. Finally we plot the distribution of survival times across datasets.

ov_surv_exprs <- purrr::map2(
  .x = list_ovarian_eset,
  .y = ov_surv_raw,
  .f = ~ exprs(.x[common_genes, .y$sample_id])
)

ov_surv_tbl <- ov_surv_raw %>%
  bind_rows(.id = "data_source")
ov_surv_tbl %>%
  ggplot(aes(
    x = time,
    y = ..density..,
    fill = data_source
  )) +
  geom_density(alpha = 0.25) +
  scale_fill_d3()
#> Warning: The dot-dot notation (`..density..`) was deprecated in ggplot2 3.4.0.
#> ℹ Please use `after_stat(density)` instead.
#> This warning is displayed once every 8 hours.
#> Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
#> generated.

# Building a TOP survival model. In this section, we will build the survival model using the TOP framework. To do this, we will first organize the survival data and gene expression data into appropriate data structures. Then, we will apply the TOP_survival function to build the model, specifying the number of features to be selected by the lasso regression.

surv_list <- ov_surv_tbl %>%
  split(ov_surv_tbl$data_source)
surv_list <- lapply(surv_list, function(x) x[, 3:4])

surv_counts <- ov_surv_exprs %>% lapply(t)
surv_list <- surv_list[names(surv_counts)]
surv_model <- TOP_survival(
  x_list = surv_counts, y_list = surv_list, nFeatures = 10
)

8 Visualising performance.

Once we have built the survival model using the TOP framework, it is important to evaluate its performance. In this section, we will visualize the performance of the model by calculating the concordance index (C-index) for each dataset. The C-index measures the agreement between the predicted and observed survival times, with values ranging from 0 to 1,

conf_results <- unlist(lapply(seq_along(surv_counts), function(x) {
  Surv_TOP_CI(
    surv_model,
    newx = surv_counts[[x]], newy = surv_list[[x]]
  )$concordance
}))

conf_results %>%
  tibble::enframe() %>%
  mutate(Metric = "C-index") %>%
  ggplot(aes(y = value, x = Metric)) +
  geom_boxplot(width = 0.5) +
  ylab("C-index") +
  geom_jitter(alpha = 0.7, width = 0.1) +
  theme(axis.text.x = element_blank()) +
  xlab("Survival Model")

9 sessionInfo

sessionInfo()
#> R version 4.4.0 beta (2024-04-15 r86425)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 22.04.4 LTS
#> 
#> Matrix products: default
#> BLAS:   /home/biocbuild/bbs-3.19-bioc/R/lib/libRblas.so 
#> LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.10.0
#> 
#> locale:
#>  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
#>  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
#>  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
#>  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
#>  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
#> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
#> 
#> time zone: America/New_York
#> tzcode source: system (glibc)
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#>  [1] curatedOvarianData_1.41.1 TOP_1.4.0                
#>  [3] ggbeeswarm_0.7.2          ggsci_3.0.3              
#>  [5] Biobase_2.64.0            BiocGenerics_0.50.0      
#>  [7] survminer_0.4.9           ggpubr_0.6.0             
#>  [9] survival_3.6-4            lubridate_1.9.3          
#> [11] forcats_1.0.0             stringr_1.5.1            
#> [13] dplyr_1.1.4               purrr_1.0.2              
#> [15] readr_2.1.5               tidyr_1.3.1              
#> [17] tibble_3.2.1              ggplot2_3.5.1            
#> [19] tidyverse_2.0.0           BiocStyle_2.32.0         
#> 
#> loaded via a namespace (and not attached):
#>   [1] splines_4.4.0               latex2exp_0.9.6            
#>   [3] polyclip_1.10-6             hardhat_1.3.1              
#>   [5] pROC_1.18.5                 rpart_4.1.23               
#>   [7] lifecycle_1.0.4             rstatix_0.7.2              
#>   [9] doParallel_1.0.17           globals_0.16.3             
#>  [11] lattice_0.22-6              MASS_7.3-60.2              
#>  [13] MultiAssayExperiment_1.30.0 backports_1.4.1            
#>  [15] magrittr_2.0.3              limma_3.60.0               
#>  [17] Hmisc_5.1-2                 plotly_4.10.4              
#>  [19] sass_0.4.9                  rmarkdown_2.26             
#>  [21] jquerylib_0.1.4             yaml_2.3.8                 
#>  [23] ClassifyR_3.8.0             abind_1.4-5                
#>  [25] zlibbioc_1.50.0             GenomicRanges_1.56.0       
#>  [27] ggraph_2.2.1                nnet_7.3-19                
#>  [29] tweenr_2.0.3                ipred_0.9-14               
#>  [31] lava_1.8.0                  GenomeInfoDbData_1.2.12    
#>  [33] IRanges_2.38.0              KMsurv_0.1-5               
#>  [35] S4Vectors_0.42.0            ggrepel_0.9.5              
#>  [37] listenv_0.9.1               parallelly_1.37.1          
#>  [39] codetools_0.2-20            DelayedArray_0.30.0        
#>  [41] ggforce_0.4.2               tidyselect_1.2.1           
#>  [43] shape_1.4.6.1               UCSC.utils_1.0.0           
#>  [45] farver_2.1.1                viridis_0.6.5              
#>  [47] matrixStats_1.3.0           stats4_4.4.0               
#>  [49] base64enc_0.1-3             jsonlite_1.8.8             
#>  [51] caret_6.0-94                tidygraph_1.3.1            
#>  [53] Formula_1.2-5               iterators_1.0.14           
#>  [55] foreach_1.5.2               tools_4.4.0                
#>  [57] ggnewscale_0.4.10           Rcpp_1.0.12                
#>  [59] glue_1.7.0                  prodlim_2023.08.28         
#>  [61] gridExtra_2.3               SparseArray_1.4.0          
#>  [63] xfun_0.43                   MatrixGenerics_1.16.0      
#>  [65] ggthemes_5.1.0              GenomeInfoDb_1.40.0        
#>  [67] withr_3.0.0                 BiocManager_1.30.22        
#>  [69] fastmap_1.1.1               fansi_1.0.6                
#>  [71] digest_0.6.35               timechange_0.3.0           
#>  [73] R6_2.5.1                    colorspace_2.1-0           
#>  [75] utf8_1.2.4                  generics_0.1.3             
#>  [77] directPA_1.5.1              calibrate_1.7.7            
#>  [79] data.table_1.15.4           recipes_1.0.10             
#>  [81] class_7.3-22                graphlayouts_1.1.1         
#>  [83] httr_1.4.7                  htmlwidgets_1.6.4          
#>  [85] S4Arrays_1.4.0              ModelMetrics_1.2.2.2       
#>  [87] pkgconfig_2.0.3             gtable_0.3.5               
#>  [89] timeDate_4032.109           XVector_0.44.0             
#>  [91] survMisc_0.5.6              htmltools_0.5.8.1          
#>  [93] carData_3.0-5               bookdown_0.39              
#>  [95] scales_1.3.0                gower_1.0.1                
#>  [97] knitr_1.46                  km.ci_0.5-6                
#>  [99] rstudioapi_0.16.0           tzdb_0.4.0                 
#> [101] reshape2_1.4.4              checkmate_2.3.1            
#> [103] nlme_3.1-164                cachem_1.0.8               
#> [105] zoo_1.8-12                  parallel_4.4.0             
#> [107] vipor_0.4.7                 foreign_0.8-86             
#> [109] pillar_1.9.0                grid_4.4.0                 
#> [111] vctrs_0.6.5                 car_3.1-2                  
#> [113] xtable_1.8-4                cluster_2.1.6              
#> [115] beeswarm_0.4.0              htmlTable_2.4.2            
#> [117] evaluate_0.23               magick_2.8.3               
#> [119] tinytex_0.50                cli_3.6.2                  
#> [121] compiler_4.4.0              rlang_1.1.3                
#> [123] crayon_1.5.2                future.apply_1.11.2        
#> [125] ggsignif_0.6.4              labeling_0.4.3             
#> [127] plyr_1.8.9                  stringi_1.8.3              
#> [129] viridisLite_0.4.2           BiocParallel_1.38.0        
#> [131] assertthat_0.2.1            munsell_0.5.1              
#> [133] lazyeval_0.2.2              glmnet_4.1-8               
#> [135] Matrix_1.7-0                hms_1.1.3                  
#> [137] future_1.33.2               statmod_1.5.0              
#> [139] highr_0.10                  SummarizedExperiment_1.34.0
#> [141] igraph_2.0.3                broom_1.0.5                
#> [143] memoise_2.0.1               bslib_0.7.0