Contents

0.1 Instalation

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

1 Required Packages

library(dplyr)
library(ggplot2)
library(survival)
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())

2 Load data

The data is loaded from an online curated dataset downloaded from TCGA using curatedTCGAData bioconductor 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.

brca <- curatedTCGAData(diseaseCode = "BRCA", assays = "RNASeq2GeneNorm", FALSE)
brca <- TCGAutils::splitAssays(brca, c('01','11'))
## Selecting 'Primary Solid Tumor' samples
## Selecting 'Solid Tissue Normal' samples
xdata.raw <- t(cbind(assay(brca[[1]]), assay(brca[[2]])))

# Get matches between survival and assay data
class.v        <- TCGAbiospec(rownames(xdata.raw))$sample_definition %>% factor
names(class.v) <- rownames(xdata.raw)

# keep features with standard deviation > 0
xdata.raw <- xdata.raw %>% 
  { (apply(., 2, sd) != 0) } %>% 
  { xdata.raw[, .] } %>%
  scale

set.seed(params$seed)
small.subset <- c('CD5', 'CSF2RB', 'HSF1', 'IRGC', 'LRRC37A6P', 'NEUROG2', 
                  'NLRC4', 'PDE11A', 'PIK3CB', 'QARS', 'RPGRIP1L', 'SDC1', 
                  'TMEM31', 'YME1L1', 'ZBTB11', 
                  sample(colnames(xdata.raw), 100))

xdata <- xdata.raw[, small.subset[small.subset %in% colnames(xdata.raw)]]
ydata <- class.v

3 Fit models

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

fitted <- cv.glmHub(xdata, ydata, 
                    family  = 'binomial',
                    network = 'correlation', 
                    nlambda = 1000,
                    network.options = networkOptions(cutoff = .6, 
                                                     min.degree = .2))

4 Results of Cross Validation

Shows the results of 1000 different 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)

4.1 Coefficients of selected model from Cross-Validation

Taking the best model described by lambda.min

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()
## Warning: Prefixing `UQ()` with the rlang namespace is deprecated as of rlang 0.3.0.
## Please use the non-prefixed form or `!!` instead.
## 
##   # Bad:
##   rlang::expr(mean(rlang::UQ(var) * 100))
## 
##   # Ok:
##   rlang::expr(mean(UQ(var) * 100))
## 
##   # Good:
##   rlang::expr(mean(!!var * 100))
## 
## This warning is displayed once per session.
ensembl.id gene.name coefficient
(Intercept) (Intercept) -6.3376277
CD5 ABCD2 -0.8421766
HSF1 ACSL1 -0.5736074
QARS ADSS 0.2015771
YME1L1 AGR2 -0.0075664
LOC143666 AMOT -0.1533910
CLDN8 ASB8 0.1750873
RBM12B C19orf23 -0.0919641
CDH26 C8orf33 -0.2110669
C19orf23 CD5 -0.4630254
GOLT1B CDH26 -0.3952287
VIM CLDN8 0.3416265
RPL17 ENTPD3 0.8270199
ASB8 FAM13C 0.6173100
ITGB1BP3 GOLT1B 1.3414562
UTP23 HLA-DPA1 -1.1789746
HLA-DPA1 HSF1 -0.1204987
USE1 ITGB1BP3 -0.5523879
WDR90 KCNS2 -0.7995671
ZNF69 LOC143666 -0.1268551
ABCD2 MUC15 6.5050567
ADSS NISCH -0.3468637
ENTPD3 OR6B1 -0.8446263
FAM13C PGAM2 1.0084372
SLC13A2 PLCB2 0.2721161
OR6B1 QARS -0.3034341
C8orf33 RBM12B -0.5743207
ACSL1 RPL17 1.5034740
TSPYL2 SLC13A2 0.6332315
MUC15 TSPYL2 0.0501196
KCNS2 TTC35 -0.9306531
TTC35 USE1 -0.1870829
AGR2 UTP23 -0.2088197
NISCH VIM 0.7051475
AMOT WDR90 0.1709998
PGAM2 YME1L1 -1.7735337
PLCB2 ZNF69 -0.1947588

4.2 Hallmarks of Cancer

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

4.3 Accuracy

## [INFO] Misclassified (6)
## [INFO]   * False primary solid tumour: 4
## [INFO]   * False normal              : 2

Histogram of predicted response