params <- list(seed = 20188102) ## ---- eval=FALSE----------------------------------------------------------- # if (!require("BiocManager")) # install.packages("BiocManager") # BiocManager::install("glmSparseNet") ## ----packages, message=FALSE, warning=FALSE-------------------------------- library(dplyr) library(Matrix) library(ggplot2) library(forcats) library(parallel) library(STRINGdb) library(reshape2) library(survival) library(VennDiagram) library(loose.rock) 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 <- loose.rock::show.message(FALSE) # Setting ggplot2 default theme as minimal theme_set(ggplot2::theme_minimal()) ## ---- eval=FALSE----------------------------------------------------------- # # Not evaluated in vignette as it takes too long to download and process # all.interactions.700 <- stringDBhomoSapiens(score_threshold = 700) # string.network <- buildStringNetwork(all.interactions.700, # use.names = 'external') ## ---- include=FALSE-------------------------------------------------------- data('string.network.700.cache', package = 'glmSparseNet') string.network <- string.network.700.cache ## -------------------------------------------------------------------------- string.network.undirected <- string.network + Matrix::t(string.network) string.network.undirected <- (string.network.undirected != 0) * 1 ## ---- echo=FALSE, collapse=TRUE-------------------------------------------- flog.info('Directed graph (score_threshold = %d)', 700) flog.info(' * total edges: %d', sum(string.network != 0)) flog.info(' * unique protein: %d', nrow(string.network)) flog.info(' * edges per protein: %f', sum(string.network != 0) / nrow(string.network) ) flog.info('') flog.info('Undirected graph (score_threshold = %d)', 700) flog.info(' * total edges: %d', sum(string.network.undirected != 0) / 2) flog.info(' * unique protein: %d', nrow(string.network.undirected)) flog.info(' * edges per protein: %f', sum(string.network.undirected != 0)/2/nrow(string.network.undirected)) ## ---- echo=FALSE----------------------------------------------------------- string.network.binary <- (string.network.undirected != 0) * 1 degree.network <- (Matrix::rowSums(string.network.binary) + Matrix::colSums(string.network.binary)) / 2 flog.info('Summary of degree:', summary(degree.network), capture = TRUE) ## ---- warning=FALSE-------------------------------------------------------- qplot(degree.network[degree.network <= quantile(degree.network, probs = .99999)], geom = 'histogram', fill = my.colors(2), bins = 100) + theme(legend.position = 'none') + xlab('Degree (up until 99.999% quantile)') ## ---- include=FALSE-------------------------------------------------------- # chunk not included as it produces to many unnecessary messages brca <- curatedTCGAData(diseaseCode = "BRCA", assays = "RNASeq2GeneNorm", FALSE) ## ---- eval=FALSE----------------------------------------------------------- # brca <- curatedTCGAData(diseaseCode = "BRCA", assays = "RNASeq2GeneNorm", FALSE) ## ----data.show, warning=FALSE, error=FALSE--------------------------------- # keep only solid tumour (code: 01) brca.primary.solid.tumor <- TCGAutils::splitAssays(brca, '01') xdata.raw <- t(assay(brca.primary.solid.tumor[[1]])) # Get survival information ydata.raw <- colData(brca.primary.solid.tumor) %>% as.data.frame %>% # Convert days to integer mutate(Days.to.date.of.Death = as.integer(Days.to.date.of.Death), Days.to.Last.Contact = as.integer(Days.to.Date.of.Last.Contact)) %>% # Find max time between all days (ignoring missings) rowwise %>% 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 select(patientID, status = vital_status, time) %>% # Discard individuals with survival time less or equal to 0 filter(!is.na(time) & time > 0) %>% as.data.frame # Set index as the patientID rownames(ydata.raw) <- ydata.raw$patientID # keep only features that are in degree.network and have standard deviation > 0 valid.features <- colnames(xdata.raw)[colnames(xdata.raw) %in% names(degree.network[degree.network > 0])] xdata.raw <- xdata.raw[TCGAbarcode(rownames(xdata.raw)) %in% rownames(ydata.raw), valid.features] xdata.raw <- scale(xdata.raw) # 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('AAK1', 'ADRB1', 'AK7', 'ALK', 'APOBEC3F', 'ARID1B', 'BAMBI', 'BRAF', 'BTG1', 'CACNG8', 'CASP12', 'CD5', 'CDA', 'CEP72', 'CPD', 'CSF2RB', 'CSN3', 'DCT', 'DLG3', 'DLL3', 'DPP4', 'DSG1', 'EDA2R', 'ERP27', 'EXD1', 'GABBR2', 'GADD45A', 'GBP1', 'HTR1F', 'IFNK', 'IRF2', 'IYD', 'KCNJ11', 'KRTAP5-6', 'MAFA', 'MAGEB4', 'MAP2K6', 'MCTS1', 'MMP15', 'MMP9', 'NFKBIA', 'NLRC4', 'NT5C1A', 'OPN4', 'OR13C5', 'OR13C8', 'OR2T6', 'OR4K2', 'OR52E6', 'OR5D14', 'OR5H1', 'OR6C4', 'OR7A17', 'OR8J3', 'OSBPL1A', 'PAK6', 'PDE11A', 'PELO', 'PGK1', 'PIK3CB', 'PMAIP1', 'POLR2B', 'POP1', 'PPFIA3', 'PSME1', 'PSME2', 'PTEN', 'PTGES3', 'QARS', 'RABGAP1', 'RBM3', 'RFC3', 'RGPD8', 'RPGRIP1L', 'SAV1', 'SDC1', 'SDC3', 'SEC16B', 'SFPQ', 'SFRP5', 'SIPA1L1', 'SLC2A14', 'SLC6A9', 'SPATA5L1', 'SPINT1', 'STAR', 'STXBP5', 'SUN3', 'TACC2', 'TACR1', 'TAGLN2', 'THPO', 'TNIP1', 'TP53', 'TRMT2B', 'TUBB1', 'VDAC1', 'VSIG8', 'WNT3A', 'WWOX', 'XRCC4', 'YME1L1', 'ZBTB11', 'ZSCAN21') %>% sample(size = 50) %>% sort # make sure we have 100 genes small.subset <- c(small.subset, sample(colnames(xdata.raw), 51)) %>% unique %>% sort xdata <- xdata.raw[, small.subset[small.subset %in% colnames(xdata.raw)]] ydata <- ydata.raw %>% select(time, status) %>% filter(!is.na(time) | time < 0) ## -------------------------------------------------------------------------- # # Add degree 0 to genes not in STRING network my.degree <- degree.network[small.subset] my.string <- string.network.binary[small.subset, small.subset] ## ---- echo=FALSE, warning=FALSE-------------------------------------------- qplot(my.degree, bins = 100, fill = my.colors(3)) + theme(legend.position = 'none') ## -------------------------------------------------------------------------- set.seed(params$seed) foldid <- balanced.cv.folds(!!ydata$status)$output ## ---- include=FALSE-------------------------------------------------------- # List that will store all selected genes selected.genes <- list() ## ---- warning=FALSE, error=FALSE------------------------------------------- cv.hub <- cv.glmHub(xdata, Surv(ydata$time, ydata$status), family = 'cox', foldid = foldid, network = my.string, network.options = networkOptions(min.degree = 0.2)) ## ---- echo=FALSE----------------------------------------------------------- glmSparseNet::separate2GroupsCox(as.vector(coef(cv.hub, s = 'lambda.min')[,1]), xdata, ydata, plot.title = 'Full dataset', legend.outside = FALSE) selected.genes[['Hub']] <- coef(cv.hub, s = 'lambda.min')[,1] %>% { .[. != 0] } %>% names %>% geneNames %>% { .[['external_gene_name']]} ## ---- warning=FALSE, error=FALSE------------------------------------------- cv.orphan <- cv.glmOrphan(xdata, Surv(ydata$time, ydata$status), family = 'cox', foldid = foldid, network = my.string, network.options = networkOptions(min.degree = 0.2)) ## ---- echo=FALSE----------------------------------------------------------- glmSparseNet::separate2GroupsCox(as.vector(coef(cv.orphan, s = 'lambda.min')[,1]), xdata, ydata, plot.title = 'Full dataset', legend.outside = FALSE) selected.genes[['Orphan']] <- coef(cv.orphan, s = 'lambda.min')[,1] %>% { .[. != 0] } %>% names %>% geneNames %>% { .[['external_gene_name']]} ## ---- warning=FALSE, error=FALSE------------------------------------------- cv.glm <- cv.glmnet(xdata, Surv(ydata$time, ydata$status), family = 'cox', foldid = foldid) ## ---- echo=FALSE----------------------------------------------------------- glmSparseNet::separate2GroupsCox(as.vector(coef(cv.glm, s = 'lambda.min')[,1]), xdata, ydata, plot.title = 'Full dataset', legend.outside = FALSE) selected.genes[['GLMnet']] <- coef(cv.glm, s = 'lambda.min')[,1] %>% { .[. != 0] } %>% names %>% geneNames %>% { .[['external_gene_name']]} ## ---- echo=FALSE, warning=FALSE-------------------------------------------- venn.plot <- venn.diagram(selected.genes , NULL, fill = c(my.colors(5), my.colors(3), my.colors(4)), alpha = c(0.3,0.3, .3), cex = 2, cat.fontface = 4, category.names = names(selected.genes)) grid.draw(venn.plot) ## ---- echo=FALSE, warning=FALSE-------------------------------------------- melt(selected.genes) %>% mutate(Degree = my.degree[value], value = factor(value), L1 = factor(L1)) %>% mutate(value = fct_reorder(value, Degree)) %>% as.data.frame %>% ggplot() + geom_point(aes(value, L1, size = Degree), shape = my.symbols(3), color = my.colors(3)) + theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0), legend.position = 'top') + ylab('Model') + xlab('Gene') + scale_size_continuous(trans = 'log10')