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.5.0 (2018-04-23)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 16.04.4 LTS
## 
## Matrix products: default
## BLAS: /home/biocbuild/bbs-3.7-bioc/R/lib/libRblas.so
## LAPACK: /home/biocbuild/bbs-3.7-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.2              SummarizedExperiment_1.10.0
##  [3] DelayedArray_0.6.0          BiocParallel_1.14.0        
##  [5] matrixStats_0.53.1          Biobase_2.40.0             
##  [7] GenomicRanges_1.32.0        GenomeInfoDb_1.16.0        
##  [9] IRanges_2.14.0              S4Vectors_0.18.0           
## [11] BiocGenerics_0.26.0         biotmleData_1.3.0          
## [13] biotmle_1.4.0               dplyr_0.7.4                
## [15] BiocStyle_2.8.0            
## 
## loaded via a namespace (and not attached):
##  [1] Rcpp_0.12.16           lattice_0.20-35        listenv_0.7.0         
##  [4] assertthat_0.2.0       rprojroot_1.3-2        digest_0.6.15         
##  [7] foreach_1.4.4          R6_2.2.2               plyr_1.8.4            
## [10] nnls_1.4               backports_1.1.2        evaluate_0.10.1       
## [13] ggplot2_2.2.1          pillar_1.2.2           zlibbioc_1.26.0       
## [16] rlang_0.2.0            lazyeval_0.2.1         Matrix_1.2-14         
## [19] rmarkdown_1.9          labeling_0.3           stringr_1.3.0         
## [22] RCurl_1.95-4.10        munsell_0.4.3          compiler_3.5.0        
## [25] xfun_0.1               pkgconfig_2.0.1        superheat_0.1.0       
## [28] globals_0.11.0         htmltools_0.3.6        tibble_1.4.2          
## [31] GenomeInfoDbData_1.1.0 bookdown_0.7           codetools_0.2-15      
## [34] doFuture_0.6.0         future_1.8.0           tmle_1.3.0-1          
## [37] MASS_7.3-50            bitops_1.0-6           grid_3.5.0            
## [40] gtable_0.2.0           magrittr_1.5           scales_0.5.0          
## [43] stringi_1.1.7          XVector_0.20.0         limma_3.36.0          
## [46] ggdendro_0.1-20        ggsci_2.8              iterators_1.0.9       
## [49] tools_3.5.0            glue_1.2.0             yaml_2.1.18           
## [52] colorspace_1.3-2       SuperLearner_2.0-23    knitr_1.20            
## [55] bindr_0.1.1