0.1Instalation

if (!require('BiocManager')) install.packages('BiocManager') BiocManager::install('glmSparseNet')

1Required Packages

library(dplyr) library(ggplot2) library(survival) library(futile.logger) library(curatedTCGAData) library(TCGAutils) # library(glmSparseNet) # # Some general options for futile.logger the debugging package .Last.value <- flog.layout(layout.format('[~l] ~m')) .Last.value <- glmSparseNet:::show.message(FALSE) # Setting ggplot2 default theme as minimal theme_set(ggplot2::theme_minimal())

2Load data

The data is loaded from an online curated dataset downloaded from TCGA usingcuratedTCGADatabioconductor package and processed.

To accelerate the process we use a very reduced dataset down to 107 variables only(genes), which is stored as a data object in this package. However, the procedure to obtain the data manually is described in the following chunk.

skcm <- curatedTCGAData(diseaseCode = 'SKCM', assays = 'RNASeq2GeneNorm', version = '1.1.38', dry.run = FALSE)

构建survival data from the clinical columns.

• Merge survival times for patients, which have different columns in case they are alive or dead.
• Build two matrix objects that fit the dataxdataandydata

skcm。转移性< - TCGAutils: TCGAsplitAssays (skcm, '06') xdata.raw <- t(assay(skcm.metastatic[[1]])) # Get survival information ydata.raw <- colData(skcm.metastatic) %>% as.data.frame %>% # Find max time between all days (ignoring missings) dplyr::rowwise() %>% dplyr::mutate( time = max(days_to_last_followup, days_to_death, na.rm = TRUE) ) %>% # Keep only survival variables and codes dplyr::select(patientID, status = vital_status, time) %>% # Discard individuals with survival time less or equal to 0 dplyr::filter(!is.na(time) & time > 0) %>% as.data.frame() # Get survival information ydata.raw <- colData(skcm) %>% as.data.frame %>% # Find max time between all days (ignoring missings) dplyr::rowwise() %>% dplyr::mutate( time = max(days_to_last_followup, days_to_death, na.rm = TRUE) ) %>% # Keep only survival variables and codes dplyr::select(patientID, status = vital_status, time) %>% # Discard individuals with survival time less or equal to 0 dplyr::filter(!is.na(time) & time > 0) %>% as.data.frame # Set index as the patientID rownames(ydata.raw) <- ydata.raw$patientID # keep only features that have standard deviation > 0 xdata.raw <- xdata.raw[TCGAbarcode(rownames(xdata.raw)) %in% rownames(ydata.raw),] xdata.raw <- xdata.raw %>% { (apply(., 2, sd) != 0) } %>% { xdata.raw[, .] } %>% scale # Order ydata the same as assay ydata.raw <- ydata.raw[TCGAbarcode(rownames(xdata.raw)), ] set.seed(params$seed) small.subset <- c('FOXL2', 'KLHL5', 'PCYT2', 'SLC6A10P', 'STRAP', 'TMEM33', 'WT1-AS', sample(colnames(xdata.raw), 100)) xdata <- xdata.raw[, small.subset[small.subset %in% colnames(xdata.raw)]] ydata <- ydata.raw %>% dplyr::select(time, status)

3Fit models

Fit model model penalizing by the hubs using the cross-validation function bycv.glmHub.

fitted <- cv.glmHub( xdata, Surv(ydata$time, ydata$status), family = 'cox', foldid = glmSparseNet:::balanced.cv.folds(!!ydata$status)$output, network = 'correlation', network.options = networkOptions(min.degree = .2, cutoff = .6) )

4Results of Cross Validation

Shows the results of100different parameters used to find the optimal value in 10-fold cross-validation. The two vertical dotted lines represent the best model and a model with less variables selected(genes), but within a standard error distance from the best.

plot(fitted)

coefs.v <- coef(fitted, s = 'lambda.min')[,1] %>% { .[. != 0]} coefs.v %>% { data.frame(ensembl.id = names(.), gene.name = geneNames(names(.))$external_gene_name, coefficient = ., stringsAsFactors = FALSE) } %>% arrange(gene.name) %>% knitr::kable()

geneNames(names(coefs.v)) %>% { hallmarks(.$external_gene_name)$heatmap }

## Error in curl::curl_fetch_memory(url, handle = handle): OpenSSL SSL_connect: SSL_ERROR_SYSCALL in connection to chat.lionproject.net:443 ## Request failed [ERROR]. Retrying in 1.2 seconds...

## Error in curl::curl_fetch_memory(url, handle = handle): OpenSSL SSL_connect: SSL_ERROR_SYSCALL in connection to chat.lionproject.net:443 ## Request failed [ERROR]. Retrying in 1.6 seconds...

## Cannot call Hallmark API, please try again later.

## NULL

separate2GroupsCox(as.vector(coefs.v), xdata[, names(coefs.v)], ydata, plot.title = 'Full dataset', legend.outside = FALSE)

## $pvalue ## [1] 0.0001269853 ## ## $plot

## ## $km ## Call: survfit(formula = survival::Surv(time, status) ~ group, data = prognostic.index.df) ## ## n events median 0.95LCL 0.95UCL ## Low risk 180 79 4000 2927 6164 ## High risk 179 114 2005 1524 2829

5Session Info

sessionInfo()

# # R版本4.2.0 RC (2022-04-21 r82226) # # Platform: x86_64-pc-linux-gnu (64-bit) ## Running under: Ubuntu 20.04.4 LTS ## ## Matrix products: default ## BLAS: /home/biocbuild/bbs-3.16-bioc/R/lib/libRblas.so ## LAPACK: /home/biocbuild/bbs-3.16-bioc/R/lib/libRlapack.so ## ## locale: ## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C ## [3] LC_TIME=en_GB LC_COLLATE=C ## [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 ## ## attached base packages: ## [1] grid parallel stats4 stats graphics grDevices utils ## [8] datasets methods base ## ## other attached packages: ## [1] VennDiagram_1.7.3 reshape2_1.4.4 ## [3] forcats_0.5.1 glmSparseNet_1.15.0 ## [5] glmnet_4.1-4 Matrix_1.4-1 ## [7] TCGAutils_1.17.0 curatedTCGAData_1.17.1 ## [9] MultiAssayExperiment_1.23.0 SummarizedExperiment_1.27.0 ## [11] Biobase_2.57.0 GenomicRanges_1.49.0 ## [13] GenomeInfoDb_1.33.0 IRanges_2.31.0 ## [15] S4Vectors_0.35.0 BiocGenerics_0.43.0 ## [17] MatrixGenerics_1.9.0 matrixStats_0.62.0 ## [19] futile.logger_1.4.3 survival_3.3-1 ## [21] ggplot2_3.3.5 dplyr_1.0.8 ## [23] BiocStyle_2.25.0 ## ## loaded via a namespace (and not attached): ## [1] backports_1.4.1 AnnotationHub_3.5.0 ## [3] BiocFileCache_2.5.0 plyr_1.8.7 ## [5] splines_4.2.0 BiocParallel_1.31.0 ## [7] digest_0.6.29 foreach_1.5.2 ## [9] htmltools_0.5.2 magick_2.7.3 ## [11] fansi_1.0.3 magrittr_2.0.3 ## [13] memoise_2.0.1 tzdb_0.3.0 ## [15] Biostrings_2.65.0 readr_2.1.2 ## [17] prettyunits_1.1.1 colorspace_2.0-3 ## [19] blob_1.2.3 rvest_1.0.2 ## [21] rappdirs_0.3.3 xfun_0.30 ## [23] crayon_1.5.1 RCurl_1.98-1.6 ## [25] jsonlite_1.8.0 zoo_1.8-10 ## [27] iterators_1.0.14 glue_1.6.2 ## [29] survminer_0.4.9 GenomicDataCommons_1.21.0 ## [31] gtable_0.3.0 zlibbioc_1.43.0 ## [33] XVector_0.37.0 DelayedArray_0.23.0 ## [35] car_3.0-12 shape_1.4.6 ## [37] abind_1.4-5 scales_1.2.0 ## [39] futile.options_1.0.1 DBI_1.1.2 ## [41] rstatix_0.7.0 Rcpp_1.0.8.3 ## [43] xtable_1.8-4 progress_1.2.2 ## [45] bit_4.0.4 km.ci_0.5-6 ## [47] httr_1.4.2 ellipsis_0.3.2 ## [49] pkgconfig_2.0.3 XML_3.99-0.9 ## [51] farver_2.1.0 sass_0.4.1 ## [53] dbplyr_2.1.1 utf8_1.2.2 ## [55] tidyselect_1.1.2 labeling_0.4.2 ## [57] rlang_1.0.2 later_1.3.0 ## [59] AnnotationDbi_1.59.0 munsell_0.5.0 ## [61] BiocVersion_3.16.0 tools_4.2.0 ## [63] cachem_1.0.6 cli_3.3.0 ## [65] generics_0.1.2 RSQLite_2.2.12 ## [67] ExperimentHub_2.5.0 broom_0.8.0 ## [69] evaluate_0.15 stringr_1.4.0 ## [71] fastmap_1.1.0 yaml_2.3.5 ## [73] knitr_1.38 bit64_4.0.5 ## [75] survMisc_0.5.6 purrr_0.3.4 ## [77] KEGGREST_1.37.0 mime_0.12 ## [79] formatR_1.12 xml2_1.3.3 ## [81] biomaRt_2.53.0 compiler_4.2.0 ## [83] filelock_1.0.2 curl_4.3.2 ## [85] png_0.1-7 interactiveDisplayBase_1.35.0 ## [87] ggsignif_0.6.3 tibble_3.1.6 ## [89] bslib_0.3.1 stringi_1.7.6 ## [91] highr_0.9 GenomicFeatures_1.49.0 ## [93] lattice_0.20-45 KMsurv_0.1-5 ## [95] vctrs_0.4.1 pillar_1.7.0 ## [97] lifecycle_1.0.1 BiocManager_1.30.17 ## [99] jquerylib_0.1.4 data.table_1.14.2 ## [101] bitops_1.0-7 httpuv_1.6.5 ## [103] rtracklayer_1.57.0 R6_2.5.1 ## [105] BiocIO_1.7.0 bookdown_0.26 ## [107] promises_1.2.0.1 gridExtra_2.3 ## [109] codetools_0.2-18 lambda.r_1.2.4 ## [111] assertthat_0.2.1 rjson_0.2.21 ## [113] withr_2.5.0 GenomicAlignments_1.33.0 ## [115] Rsamtools_2.13.0 GenomeInfoDbData_1.2.8 ## [117] hms_1.1.1 tidyr_1.2.0 ## [119] rmarkdown_2.14 carData_3.0-5 ## [121] ggpubr_0.4.0 pROC_1.18.0 ## [123] shiny_1.7.1 restfulr_0.0.13