params <- list(seed = 41) ## ---- eval=FALSE-------------------------------------------------------------- # if (!require("BiocManager")) # install.packages("BiocManager") # BiocManager::install("glmSparseNet") ## ----packages, message=FALSE, warning=FALSE----------------------------------- 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()) ## ----curated_data, include=FALSE---------------------------------------------- # chunk not included as it produces to many unnecessary messages brca <- curatedTCGAData(diseaseCode = "BRCA", assays = "RNASeq2GeneNorm", version = "1.1.38", dry.run = FALSE) ## ----curated_data_non_eval, eval=FALSE---------------------------------------- # brca <- curatedTCGAData(diseaseCode = "BRCA", assays = "RNASeq2GeneNorm", # version = "1.1.38", dry.run = FALSE) ## ----data, warning=FALSE, message=FALSE--------------------------------------- # keep only solid tumour (code: 01) brca.primary.solid.tumor <- TCGAutils::TCGAsplitAssays(brca, '01') xdata.raw <- t(assay(brca.primary.solid.tumor[[1]])) # Get survival information ydata.raw <- colData(brca.primary.solid.tumor) %>% as.data.frame %>% # Keep only data relative to survival or samples dplyr::select(patientID, vital_status, Days.to.date.of.Death, Days.to.Date.of.Last.Contact, days_to_death, days_to_last_followup, Vital.Status) %>% # Convert days to integer dplyr::mutate(Days.to.date.of.Death = as.integer(Days.to.date.of.Death)) %>% dplyr::mutate( Days.to.Last.Contact = as.integer(Days.to.Date.of.Last.Contact) ) %>% # Find max time between all days (ignoring missings) dplyr::rowwise() %>% dplyr::mutate( time = max(days_to_last_followup, Days.to.date.of.Death, Days.to.Last.Contact, 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 # Get matches between survival and assay data 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)), ] # Using only a subset of genes previously selected to keep this short example. set.seed(params$seed) small.subset <- c('CD5', 'CSF2RB', 'IRGC', 'NEUROG2', 'NLRC4', 'PDE11A', 'PTEN', 'TP53', 'BRAF', 'PIK3CB', 'QARS', 'RFC3', 'RPGRIP1L', 'SDC1', 'TMEM31', 'YME1L1', 'ZBTB11', sample(colnames(xdata.raw), 100)) %>% unique xdata <- xdata.raw[, small.subset[small.subset %in% colnames(xdata.raw)]] ydata <- ydata.raw %>% dplyr::select(time, status) ## ----fit---------------------------------------------------------------------- set.seed(params$seed) fitted <- cv.glmHub(xdata, Surv(ydata$time, ydata$status), family = 'cox', lambda = buildLambda(1), network = 'correlation', network.options = networkOptions(cutoff = .6, min.degree = .2)) ## ----results------------------------------------------------------------------ plot(fitted) ## ----show_coefs--------------------------------------------------------------- coefs.v <- coef(fitted, s = 'lambda.min')[,1] %>% { .[. != 0]} coefs.v %>% { data.frame(gene.name = names(.), coefficient = ., stringsAsFactors = FALSE) } %>% arrange(gene.name) %>% knitr::kable() ## ----hallmarks---------------------------------------------------------------- names(coefs.v) %>% { hallmarks(.)$heatmap } ## ----------------------------------------------------------------------------- separate2GroupsCox(as.vector(coefs.v), xdata[, names(coefs.v)], ydata, plot.title = 'Full dataset', legend.outside = FALSE) ## ----sessionInfo-------------------------------------------------------------- sessionInfo()