if (!require("BiocManager")) {
install.packages("BiocManager")
}
BiocManager::install("glmSparseNet")
library(dplyr)
library(Matrix)
library(ggplot2)
library(forcats)
library(parallel)
library(reshape2)
library(survival)
library(VennDiagram)
library(futile.logger)
library(curatedTCGAData)
library(MultiAssayExperiment)
library(TCGAutils)
#
library(glmSparseNet)
#
#
# Some general options for futile.logger the debugging package
flog.layout(layout.format("[~l] ~m"))
options(
"glmSparseNet.show_message" = FALSE,
"glmSparseNet.base_dir" = withr::local_tempdir()
)
# Setting ggplot2 default theme as minimal
theme_set(ggplot2::theme_minimal())
This vignette uses the STRING database (https://string-db.org/) of protein-protein interactions as the network-based penalizer in generalized linear models using Breast invasive carcinoma sample dataset.
The degree vector is calculated manually to account for genes that are not present in the STRING database, as these will not have any interactions, i.e. edges.
Retrieve all interactions from STRING databse. We have included a helper function that retrieves the Homo sapiens known interactions.
For this vignette, we use a cached version of all interaction with
score_threshold = 700
Note: Text-based interactions are excluded from the network.
# Not evaluated in vignette as it takes too long to download and process
allInteractions700 <- stringDBhomoSapiens(scoreThreshold = 700)
stringNetwork <- buildStringNetwork(allInteractions700, "external")
Build a sparse matrix object that contains the network.
stringNetworkUndirected <- stringNetwork + Matrix::t(stringNetwork)
stringNetworkUndirected <- (stringNetworkUndirected != 0) * 1
## [INFO] Directed graph (score_threshold = 700)
## [INFO] * total edges: 225330
## [INFO] * unique protein: 11033
## [INFO] * edges per protein: 20.423276
## [INFO]
## [INFO] Undirected graph (score_threshold = 700)
## [INFO] * total edges: 225209
## [INFO] * unique protein: 11033
## [INFO] * edges per protein: 20.412309
## [INFO] Summary of degree:
##
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1.00 3.00 13.00 40.83 43.00 4175.00
qplot(
degreeNetworkVector[
degreeNetworkVector <= quantile(degreeNetworkVector, probs = .99999)
],
geom = "histogram", fill = my.colors(2), bins = 100
) +
theme(legend.position = "none") + xlab("Degree (up until 99.999% quantile)")
glmSparseNet
brca <- curatedTCGAData(
diseaseCode = "BRCA", assays = "RNASeq2GeneNorm",
version = "1.1.38", dry.run = FALSE
)
Build the survival data from the clinical columns.
xdata
and ydata
# keep only solid tumour (code: 01)
brcaPrimarySolidTumor <- TCGAutils::TCGAsplitAssays(brca, "01")
xdataRaw <- t(assay(brcaPrimarySolidTumor[[1]]))
# Get survival information
ydataRaw <- colData(brcaPrimarySolidTumor) |>
as.data.frame() |>
# Convert days to integer
dplyr::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)
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(ydataRaw) <- ydataRaw$patientID
# keep only features that are in degreeNetworkVector and have
# standard deviation > 0
validFeatures <- colnames(xdataRaw)[
colnames(xdataRaw) %in% names(degreeNetworkVector[degreeNetworkVector > 0])
]
xdataRaw <- xdataRaw[
TCGAbarcode(rownames(xdataRaw)) %in% rownames(ydataRaw), validFeatures
]
xdataRaw <- scale(xdataRaw)
# Order ydata the same as assay
ydataRaw <- ydataRaw[TCGAbarcode(rownames(xdataRaw)), ]
# Using only a subset of genes previously selected to keep this short example.
set.seed(params$seed)
smallSubset <- 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
smallSubset <- c(smallSubset, sample(colnames(xdataRaw), 51)) |>
unique() |>
sort()
xdata <- xdataRaw[, smallSubset[smallSubset %in% colnames(xdataRaw)]]
ydata <- ydataRaw |>
dplyr::select(time, status) |>
dplyr::filter(!is.na(time) | time < 0)
#
# Add degree 0 to genes not in STRING network
myDegree <- degreeNetworkVector[smallSubset]
myString <- stringNetworkBinary[smallSubset, smallSubset]
Degree distribution for sample set of gene features (in xdata).
set.seed(params$seed)
foldid <- glmSparseNet:::balancedCvFolds(ydata$status)$output
Penalizes using the Hub heuristics, see hubHeuristic
function definition for
more details.
resultCVHub <- cv.glmHub(xdata,
Surv(ydata$time, ydata$status),
family = "cox",
foldid = foldid,
network = myString,
network.options = networkOptions(minDegree = 0.2)
)
Kaplan-Meier estimator separating individuals by low and high risk (based on model’s coefficients)
## Warning: The `plot.title` argument of `separate2GroupsCox()` is deprecated as of
## glmSparseNet 1.21.0.
## ℹ Please use the `plotTitle` argument instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## Warning: The `legend.outside` argument of `separate2GroupsCox()` is deprecated as of
## glmSparseNet 1.21.0.
## ℹ Please use the `legendOutside` argument instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## $pvalue
## [1] 9.628346e-05
##
## $plot
##
## $km
## Call: survfit(formula = survival::Surv(time, status) ~ group, data = prognosticIndexDf)
##
## n events median 0.95LCL 0.95UCL
## Low risk - 1 983 129 3959 3738 7455
## High risk - 1 97 23 2911 2469 NA
Penalizes using the Orphan heuristics, see orphanHeuristic
function
definition for more details.
resultCVOrphan <- cv.glmOrphan(xdata,
Surv(ydata$time, ydata$status),
family = "cox",
foldid = foldid,
network = myString,
network.options = networkOptions(minDegree = 0.2)
)
Kaplan-Meier estimator separating individuals by low and high risk (based on model’s coefficients)
## $pvalue
## [1] 0
##
## $plot
##
## $km
## Call: survfit(formula = survival::Surv(time, status) ~ group, data = prognosticIndexDf)
##
## n events median 0.95LCL 0.95UCL
## Low risk - 1 540 34 7455 6456 NA
## High risk - 1 540 118 2866 2573 3472
Uses regular glmnet model as simple baseline
library(glmnet)
## Loaded glmnet 4.1-8
resultCVGlmnet <- cv.glmnet(xdata,
Surv(ydata$time, ydata$status),
family = "cox",
foldid = foldid
)
Kaplan-Meier estimator separating individuals by low and high risk (based on model’s coefficients)
## $pvalue
## [1] 0
##
## $plot
##
## $km
## Call: survfit(formula = survival::Surv(time, status) ~ group, data = prognosticIndexDf)
##
## n events median 0.95LCL 0.95UCL
## Low risk - 1 540 38 6593 4398 NA
## High risk - 1 540 114 2854 2551 3461
Venn diagram of overlapping genes.