Contents

library(dplyr)
library(biotmle)
library(biotmleData)
suppressMessages(library(SummarizedExperiment))
"%ni%" = Negate("%in%")

0.1 Introduction

Here, we briefly work through how to use the biotmle package with data generated by next-generation sequencing technologies, which, in contrast to microarray technologies, produce measurements in the form of discrete counts.

0.2 Simulation: Data Structure

set.seed(6423709)
n <- 50
g <- 2500
cases_pois <- 50
controls_pois <- 10

ngs_cases <- as.data.frame(matrix(replicate(n, rpois(g, cases_pois)), g))
ngs_controls <- as.data.frame(matrix(replicate(n, rpois(g, controls_pois)), g))

ngs_data <- as.data.frame(cbind(ngs_cases, ngs_controls))
exp_var <- c(rep(1, n), rep(0, n))
batch <- rep(1:2, n)
covar <- rep(1, n * 2)
design <- as.data.frame(cbind(exp_var, batch, covar))

head(ngs_data[, 1:7])
##   V1 V2 V3 V4 V5 V6 V7
## 1 69 61 39 46 50 57 43
## 2 42 59 50 54 44 53 54
## 3 41 49 52 54 50 58 34
## 4 30 44 49 44 46 36 61
## 5 50 56 45 44 46 61 58
## 6 61 50 47 53 49 54 67

0.3 Processing: Data Transformation

se <- SummarizedExperiment(assays = list(counts = DataFrame(ngs_data)),
                           colData = DataFrame(design))
se
## class: SummarizedExperiment 
## dim: 2500 100 
## metadata(0):
## assays(1): counts
## rownames: NULL
## rowData names(0):
## colnames(100): V1 V2 ... V49.1 V50.1
## colData names(3): exp_var batch covar

0.4 Analysis: Assessing the Effect of Exposure

rnaseqTMLEout <- biomarkertmle(se = se,
                               varInt = 1,
                               type = "exposure",
                               ngscounts = TRUE,
                               parallel = TRUE,
                               family = "gaussian",
                               g_lib = c("SL.mean", "SL.glm", "SL.randomForest"),
                               Q_lib = c("SL.mean", "SL.glm", "SL.randomForest",
                                         "SL.nnet")
                              )
head(rnaseqTMLEout@tmleOut$E[, seq_len(6)])
##                [,1]       [,2]      [,3]       [,4]      [,5]       [,6]
## result.1 -329.83655 -212.28070 160.39144   32.96832 -36.05786 -146.77124
## result.2   91.43805 -129.00196 -39.20680  -47.18052  44.29473  -30.83862
## result.3  102.01072 -294.65959 -77.65629 -376.36991 -61.66602 -441.74551
## result.4  323.94443   93.76761  13.80324   93.71017  47.56326  224.33603
## result.5   11.71456 -148.33216  93.43461   47.79556  61.82086 -229.99698
## result.6 -183.73280   -4.00971  45.02578  -53.02987  -4.00971  -69.37162
limmaTMLEout <- modtest_ic(biotmle = rnaseqTMLEout)
head(limmaTMLEout@topTable)
##               logFC    AveExpr          t    P.Value adj.P.Val         B
## result.1 -32.072167 -32.072167 -1.7430776 0.08297431 0.9303074 -4.595093
## result.2 -25.903359 -25.903359 -1.4289259 0.15470453 0.9342061 -4.595106
## result.3 -15.865179 -15.865179 -0.6675294 0.50526227 0.9731485 -4.595127
## result.4  20.394523  20.394523  0.9768991 0.32988929 0.9401033 -4.595120
## result.5 -23.932669 -23.932669 -1.2276717 0.22112330 0.9401033 -4.595113
## result.6  -8.540187  -8.540187 -0.4839883 0.62896415 0.9746949 -4.595130
##               IDs
## result.1 result.1
## result.2 result.2
## result.3 result.3
## result.4 result.4
## result.5 result.5
## result.6 result.6

0.5 Results: Data Visualization

plot(x = limmaTMLEout, type = "pvals_adj")

plot(x = limmaTMLEout, type = "pvals_raw")

varInt_index <- which(names(colData(se)) %in% "exp_var")
designVar <- as.data.frame(colData(se))[, varInt_index]
design <- as.numeric(designVar == max(designVar))

heatmap_ic(x = limmaTMLEout, design = design, FDRcutoff = 1.0, top = 10)

volcano_ic(biotmle = limmaTMLEout)

0.6 Session Information

## R version 3.4.2 (2017-09-28)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 16.04.3 LTS
## 
## Matrix products: default
## BLAS: /home/biocbuild/bbs-3.6-bioc/R/lib/libRblas.so
## LAPACK: /home/biocbuild/bbs-3.6-bioc/R/lib/libRlapack.so
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=C              
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] parallel  stats4    stats     graphics  grDevices utils     datasets 
## [8] methods   base     
## 
## other attached packages:
##  [1] bindrcpp_0.2               SummarizedExperiment_1.8.0
##  [3] DelayedArray_0.4.0         matrixStats_0.52.2        
##  [5] Biobase_2.38.0             GenomicRanges_1.30.0      
##  [7] GenomeInfoDb_1.14.0        IRanges_2.12.0            
##  [9] S4Vectors_0.16.0           BiocGenerics_0.24.0       
## [11] biotmleData_1.2.0          biotmle_1.3.0             
## [13] dplyr_0.7.4                BiocStyle_2.6.0           
## 
## loaded via a namespace (and not attached):
##  [1] Rcpp_0.12.13            lattice_0.20-35         listenv_0.6.0          
##  [4] assertthat_0.2.0        rprojroot_1.2           digest_0.6.12          
##  [7] foreach_1.4.3           R6_2.2.2                wesanderson_0.3.2      
## [10] plyr_1.8.4              nnls_1.4                backports_1.1.1        
## [13] evaluate_0.10.1         ggplot2_2.2.1           zlibbioc_1.24.0        
## [16] rlang_0.1.2             lazyeval_0.2.1          Matrix_1.2-11          
## [19] rmarkdown_1.6           labeling_0.3            BiocParallel_1.12.0    
## [22] stringr_1.2.0           RCurl_1.95-4.8          munsell_0.4.3          
## [25] compiler_3.4.2          pkgconfig_2.0.1         superheat_0.1.0        
## [28] globals_0.10.3          htmltools_0.3.6         tibble_1.3.4           
## [31] GenomeInfoDbData_0.99.1 bookdown_0.5            codetools_0.2-15       
## [34] doFuture_0.6.0          future_1.6.2            tmle_1.2.0-5           
## [37] MASS_7.3-47             bitops_1.0-6            grid_3.4.2             
## [40] gtable_0.2.0            magrittr_1.5            scales_0.5.0           
## [43] stringi_1.1.5           XVector_0.18.0          limma_3.34.0           
## [46] ggdendro_0.1-20         iterators_1.0.8         tools_3.4.2            
## [49] glue_1.2.0              yaml_2.1.14             colorspace_1.3-2       
## [52] SuperLearner_2.0-22     knitr_1.17              bindr_0.1