Contents

library(data.table)
library(DESeq2)
library(ggplot2)
library(dplyr)
library(magrittr)
library(plyr)
library(cowplot)
library(reshape2)
options(stringsAsFactors = FALSE)
library(graper)
library(tidyverse)
library("wesanderson")
library(ggbeeswarm)
source("../util_defs.R") # contains color schemes and methods to include

1 Perparations

Set input/output paths.

datadir <- "data"
outdir <- "2018-11-05"
knitr::opts_chunk$set(fig.path = "figs/", dev = c('png',"pdf"))

Load data used for fitting.

load(file.path(datadir, "dataGTEx_t5_PCs.RData"))
dim(dataGTEx$X)
## [1] 251 250

Load summarised results.

load(file.path(outdir,"result_GTEx.RData"))

Colors for tissues

cols4groups <- c(wes_palette("GrandBudapest1"),
                 wes_palette("GrandBudapest2")[1])
names(cols4groups) <- unique(dataGTEx$annot)
names(cols4groups) <- sub(" ", "\n ", names(cols4groups))

2 Make plots

2.1 RMSE

Compare prediciton performance between the methods in terms of root mean squared error.

df_RMSE <- melt(sapply(resultList, function(l) l$RMSE),
                varnames = c("method", "run"), value.name="RMSE")
df_RMSE %<>% mutate(method = make_nicenames(method))
df_RMSE %<>% mutate(method_type = ifelse(method%in% methods2compare_sparse,
                                         "sparse", "dense"))  
df_RMSE %<>% filter(method %in% methods2compare_sparse |
                      method %in% methods2compare_dense) 

# make plot 
ggRMSE <- ggplot(df_RMSE, aes(x=method, y=RMSE, fill=method)) +
  geom_boxplot(alpha=0.5, outlier.shape = NA) + 
  # ggtitle ("Prediction of donors' age")  +
  geom_beeswarm(cex = 2.5) + facet_wrap(~method_type, scale="free_x") +
  scale_fill_manual(values = cols4methods) +
  theme(axis.text.x = element_text(angle=60, vjust=1, hjust=1),
        axis.title.x = element_blank()) + guides(fill=FALSE) 
ggRMSE

2.2 Penalty Factors

## Penalty Factors and sparsity level (only in sparse graper)
df_pf <- melt(lapply(resultList, function(l) l$pf_mat[, "graper_SS", drop=F]),
              varnames = c("tissue", "method"))
df_pf$value <- as.numeric(df_pf$value)
df_pf %<>% mutate(tissue=ifelse(tissue == "Adipose Tissue", "Adipose\n Tissue",
                                ifelse(tissue == "Blood Vessel", "Blood\n Vessel", as.character(tissue))))
gg1 <- ggplot(df_pf,aes(x=tissue, y=value)) + 
    geom_boxplot(alpha=0.5, outlier.shape = NA) +
  ylab(expression(hat(gamma))) + geom_beeswarm(aes(col=tissue), cex=2.5)+
  guides(col=FALSE) + xlab("tissue type") + 
  theme_bw(base_size = 15) +
  scale_color_manual(values = cols4groups)

df_pf <- melt(lapply(resultList, function(l) l$sparsity_mat[, "graper_SS", drop=F]),
              varnames = c("tissue", "method"))
df_pf %<>% mutate(tissue=ifelse(tissue == "Adipose Tissue", "Adipose\n Tissue",
                                ifelse(tissue == "Blood Vessel", "Blood\n Vessel", as.character(tissue))))
gg2 <- ggplot(df_pf,aes(x=tissue, y=value)) +
    geom_boxplot(alpha=0.5, outlier.shape = NA) +
  ylab(expression(hat(pi))) + geom_beeswarm(aes(col=tissue), cex=2.5)+
  xlab("tissue type") + guides(col=FALSE) +
  theme_bw(base_size = 15) + scale_color_manual(values = cols4groups)

ggPF <- plot_grid(gg1,gg2, rel_widths = c(1,1.25), nrow=1, align = "hv", axis = "lb")
ggPF

2.3 Joint plot

plot_grid(ggRMSE, ggPF, ncol=1, labels = letters[1:2],
          label_size = 18, rel_heights = c(2,1))

3 Runtime

time_mat <- sapply(resultList, function(l) l$runtime)
df_time <- melt(time_mat,
                varname=c("method", "fold"), value.name="time")

# only compare to a single iteration for each method (fitted using 3)
df_time %<>% mutate(time=ifelse(grepl("graper", method), time/3/60, time/60))

# nice method names
df_time %<>% mutate(method = make_nicenames(method))
df_time %<>% filter(method %in% c(methods2compare_dense, methods2compare_sparse))

ggplot(df_time, aes(x=method, y=time, col=method)) +
  stat_summary() + ylab("time [min]") + 
  theme_bw(base_size = 15) +
  theme(axis.title.x = element_blank(),
        axis.text.x = element_text(angle=60, vjust=1, hjust=1))+
  scale_color_manual(values = cols4methods)
## No summary function supplied, defaulting to `mean_se()

df_time %>% group_by(method) %>% dplyr::summarize(mean(time))
## # A tibble: 12 x 2
##    method                  `mean(time)`
##    <chr>                          <dbl>
##  1 adaptive Lasso               0.0103 
##  2 elastic net                  0.00618
##  3 graper (dense, multiv.)      0.0844 
##  4 graper (dense)               0.00230
##  5 graper (sparse)              0.00239
##  6 group Lasso                  0.00991
##  7 GRridge                      0.607  
##  8 IPF-Lasso                   36.9    
##  9 Lasso                        0.00684
## 10 ridge regression             0.00676
## 11 sparse group Lasso           4.77   
## 12 varbvs                       0.0299

4 SessionInfo

sessionInfo()
## R version 3.5.3 (2019-03-11)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS Mojave 10.14.4
## 
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_IE.UTF-8/en_IE.UTF-8/en_IE.UTF-8/C/en_IE.UTF-8/en_IE.UTF-8
## 
## attached base packages:
## [1] parallel  stats4    stats     graphics  grDevices utils     datasets 
## [8] methods   base     
## 
## other attached packages:
##  [1] RColorBrewer_1.1-2          ggbeeswarm_0.6.0           
##  [3] wesanderson_0.3.6           forcats_0.3.0              
##  [5] stringr_1.4.0               purrr_0.3.2                
##  [7] readr_1.3.1                 tidyr_0.8.3                
##  [9] tibble_2.1.1                tidyverse_1.2.1            
## [11] graper_0.99.0               reshape2_1.4.3             
## [13] cowplot_0.9.4               plyr_1.8.4                 
## [15] magrittr_1.5                dplyr_0.8.0.1              
## [17] ggplot2_3.1.0               DESeq2_1.22.2              
## [19] SummarizedExperiment_1.12.0 DelayedArray_0.8.0         
## [21] BiocParallel_1.16.5         matrixStats_0.54.0         
## [23] Biobase_2.42.0              GenomicRanges_1.34.0       
## [25] GenomeInfoDb_1.18.1         IRanges_2.16.0             
## [27] S4Vectors_0.20.1            BiocGenerics_0.28.0        
## [29] data.table_1.12.0           BiocStyle_2.10.0           
## 
## loaded via a namespace (and not attached):
##  [1] colorspace_1.4-1       htmlTable_1.13.1       XVector_0.22.0        
##  [4] base64enc_0.1-3        rstudioapi_0.9.0       bit64_0.9-7           
##  [7] AnnotationDbi_1.44.0   fansi_0.4.0            lubridate_1.7.4       
## [10] xml2_1.2.0             splines_3.5.3          geneplotter_1.60.0    
## [13] knitr_1.21             Formula_1.2-3          jsonlite_1.6          
## [16] broom_0.5.1            annotate_1.60.0        cluster_2.0.7-1       
## [19] BiocManager_1.30.4     compiler_3.5.3         httr_1.4.0            
## [22] backports_1.1.3        assertthat_0.2.1       Matrix_1.2-15         
## [25] lazyeval_0.2.2         cli_1.1.0              acepack_1.4.1         
## [28] htmltools_0.3.6        tools_3.5.3            gtable_0.2.0          
## [31] glue_1.3.1             GenomeInfoDbData_1.2.0 Rcpp_1.0.1            
## [34] cellranger_1.1.0       nlme_3.1-137           xfun_0.4              
## [37] rvest_0.3.2            XML_3.98-1.16          zlibbioc_1.28.0       
## [40] scales_1.0.0           hms_0.4.2              yaml_2.2.0            
## [43] memoise_1.1.0          gridExtra_2.3          rpart_4.1-13          
## [46] latticeExtra_0.6-28    stringi_1.4.3          RSQLite_2.1.1         
## [49] genefilter_1.64.0      checkmate_1.9.1        rlang_0.3.2           
## [52] pkgconfig_2.0.2        bitops_1.0-6           evaluate_0.12         
## [55] lattice_0.20-38        htmlwidgets_1.3        labeling_0.3          
## [58] bit_1.1-14             tidyselect_0.2.5       bookdown_0.9          
## [61] R6_2.4.0               generics_0.0.2         Hmisc_4.2-0           
## [64] DBI_1.0.0              pillar_1.3.1           haven_2.0.0           
## [67] foreign_0.8-71         withr_2.1.2            survival_2.43-3       
## [70] RCurl_1.95-4.12        nnet_7.3-12            modelr_0.1.2          
## [73] crayon_1.3.4           utf8_1.1.4             rmarkdown_1.11        
## [76] locfit_1.5-9.1         grid_3.5.3             readxl_1.2.0          
## [79] blob_1.1.1             digest_0.6.18          xtable_1.8-3          
## [82] munsell_0.5.0          beeswarm_0.2.3         vipor_0.4.5