biotmle R package can be used to isolate biomarkers in two ways: based on the associations of genomic objects with an exposure variable of interest. In this vignette, we illustrate how to use
biotmle to isolate and visualize genes associated with an exposure, using a data set containing microarray expression measures from an Illumina platform. In the analysis described below, targeted maximum likelihood estimation (TMLE) is used to transform the microarray expression values based on the influence curve representation of the Average Treatment Effect (ATE). Following this transformation, the moderated t-statistic of Smyth (Smyth 2004) is used to test for a binary group-wise difference (based on the exposure variable), using the tools provided by the R package
For a general discussion of the framework of targeted maximum likelihood estimation and the role this approach plays in statistical causal inference, the interested reader is invited to consult van der Laan and Rose (2011) and van der Laan and Rose (2018). For a more general introduction to the principles of statistical causal inference, Pearl (2000) serves well.
First, we load the
biotmle package and the (included)
illuminaData data set:
library(dplyr) library(biotmle) library(biotmleData) library(BiocParallel) library(SuperLearner) library(SummarizedExperiment, quietly=TRUE) data(illuminaData) set.seed(13847)
In order to perform Targeted Minimum Loss-Based Estimation, we need three separate data structures: (1) W, baseline covariates that could potentially confound the association of biomarkers with the exposure of interest; (2) A, the exposure of interest; and (3) Y, the biomarkers of interest. All values in W and A ought to be discretized, in order to avoid practical violations of the assumption of positivity. With the
illuminaData data set below, we discretize the age variable in the phenotype-level data below (this can be accessed via the
colData of the
SummarizedExperiment object). To invoke the biomarker assessment function (
biomarkertmle), we also need to specify a variable of interest (or the position of said variable in the design matrix). We do both in just a few lines below:
The TMLE-based biomarker discovery process can be invoked using the
biomarkertmle function. The procedure is quite resource-intensive because it evaluates the association of each individual potential biomarker (of which there are over 20,000 in the included data set) with an exposure of interest, while accounting for potential confounding based on all other covariates included in the design matrix. We demonstrate the necessary syntax for calling the
biomarkertmle function below, on a small number of the probes:
# compute TML estimates to evaluate differentially expressed biomarkers biotmle_out <- biomarkertmle(se = illuminaData[1:20, ], varInt = benz_idx, g_lib = c("SL.mean", "SL.glm"), Q_lib = c("SL.bayesglm", "SL.ranger"), cv_folds = 2, bppar_type = SerialParam() )
## | | | 0% | |==== | 5% | |======= | 10% | |========== | 15% | |============== | 20% | |================== | 25% | |===================== | 30% | |======================== | 35% | |============================ | 40% | |================================ | 45% | |=================================== | 50% | |====================================== | 55% | |========================================== | 60% | |============================================== | 65% | |================================================= | 70% | |==================================================== | 75% | |======================================================== | 80% | |============================================================ | 85% | |=============================================================== | 90% | |================================================================== | 95% | |======================================================================| 100%
The output of
biomarkertmle is an object of class
bioTMLE, containing four new slots: (1)
call, the call to
topTable, an empty slot meant to hold the output of
limma::topTable, after a later call to
modtest_ic; and (3)
data.frame containing the point estimates of the associations of each biomarker with the exposure of interest based on the influence curve representation of the Average Treatment Effect.
The output of
biomarkertmle can be directly fed to
modtest_ic, a wrapper around
limma::topTable that outputs a
biotmle object with the slots described above completely filled in. The
modtest_ic function requires as input a
biotmle object containing a data frame in the
tmleOut field as well as a design matrix indicating the groupwise difference to be tested. The design matrix should contain an intercept term and a term for the exposure of interest (with discretized exposure levels). Based on the relevant statistical theory, it is not appropriate to include any further terms in the design matrix (n.b., this differs from standard calls to
modtmle_out <- modtest_ic(biotmle = biotmle_out)
modtest_ic, the resultant
bioTMLE object will contain all information relevant to the analytic procedure for identifying biomarkers: that is, it will contain the original call to
biomarkertmle, the result of running
limma::topTable, and the result of running
biomarkertmle. The statistical results of this procedure can be extracted from the
topTable object in the
bioTMLE object produced by
This package provides several plotting methods that can be used to visualize the results of the TMLE-based biomarker discovery process. We demonstrate the syntax for calling the generic plotting methods below but refrain from showing the plots themselves since they are not particularly informative.
plot method for a
bioTMLE object will produce a histogram of the adjusted p-values of each biomarker (based on the Benjamini-Hochberg procedure for controlling the False Discovery Rate) as generated by
plot(x = modtmle_out, type = "pvals_adj")
Setting the argument
type = "pvals_raw" will instead produce a histogram of the raw p-values (these are less informative and should, in general, not be used for inferential purposes, as the computation producing these p-values ignores the multiple testing nature of the biomarker discovery problem):
plot(x = modtmle_out, type = "pvals_raw")
Heatmaps are useful graphics for visualizing the relationship between measures on genomic objects and covariates of interest. The
heatmap_ic function provides this graphic for
bioTMLE objects, allowing for the relationship between the exposure variable and some number of “top” biomarkers (as determined by the call to
modtest_ic) to be visualized. In general, the heatmap for
bioTMLE objects expresses how the contributions of each biomarker to the Average Treatment Effect (ATE) vary across differences in the exposure variable (that is, there is a causal interpretation to the findings). The plot produced is a
ggplot2 object and can be modified in place if stored properly. For our analysis:
benz_idx <- which(names(colData(illuminaData)) %in% "benzene") designVar <- as.data.frame(colData(illuminaData))[, benz_idx] designVar <- as.numeric(designVar == max(designVar)) # build heatmap heatmap_ic(x = modtmle_out, left.label = "none", scale = TRUE, clustering.method = "hierarchical", row.dendrogram = TRUE, design = designVar, FDRcutoff = 1, top = 10)
The heatmap produced in this way is actually a form of supervised clustering, as described more generally (as supervised distance matrices) by Pollard and van der Laan (2008), wherein the notion of deriving clustering procedures from the results of supervised learning methods is formulated. Since the heatmap is based on the contributions of observations to the efficient influence function (EIF) of the target parameter, it directly visualizes the degree to which each biomarker informs the difference (due to the treatment effect) represented by the average treatment effect.
The volcano plot is standard graphical tools for examining how changes in expression relate to the raw p-value. The utility of such plots lies in their providing a convenient way to identify and systematically ignore those genomic objects that have extremely low p-values due to extremely low variance between observations. The
volcano_ic function provides much of the same interpretation, except that the fold change values displayed in the x-axis refer to changes in the contributions of each biomarker to the Average Treatment Effect (in standard practice, for microarray technology, these would be fold changes in gene expression). The plot produced is a
ggplot2 object and, as such, can be modified in place. For our analysis:
volcano_ic(biotmle = modtmle_out)
## R version 4.1.1 (2021-08-10) ## Platform: x86_64-pc-linux-gnu (64-bit) ## Running under: Ubuntu 20.04.3 LTS ## ## Matrix products: default ## BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 ## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/liblapack.so.3 ## ## locale: ##  LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C ##  LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8 ##  LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8 ##  LC_PAPER=en_US.UTF-8 LC_NAME=C ##  LC_ADDRESS=C LC_TELEPHONE=C ##  LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C ## ## attached base packages: ##  parallel stats4 splines stats graphics grDevices utils ##  datasets methods base ## ## other attached packages: ##  nloptr_126.96.36.199 quadprog_1.5-8 ##  SummarizedExperiment_1.22.0 Biobase_2.52.0 ##  GenomicRanges_1.44.0 GenomeInfoDb_1.28.4 ##  IRanges_2.26.0 S4Vectors_0.30.2 ##  BiocGenerics_0.38.0 MatrixGenerics_1.4.3 ##  matrixStats_0.61.0 SuperLearner_2.0-28 ##  gam_1.20 foreach_1.5.1 ##  nnls_1.4 BiocParallel_1.26.2 ##  biotmleData_1.16.0 biotmle_1.17.0 ##  dplyr_1.0.7 BiocStyle_2.20.2 ## ## loaded via a namespace (and not attached): ##  minqa_1.2.4 cubature_188.8.131.52 colorspace_2.0-2 ##  ellipsis_0.3.2 rprojroot_2.0.2 htmlTable_2.3.0 ##  XVector_0.32.0 ggdendro_0.1.22 base64enc_0.1-3 ##  fs_1.5.0 rstudioapi_0.13 farver_2.1.0 ##  listenv_0.8.0 MatrixModels_0.5-0 fansi_0.5.0 ##  ranger_0.13.1 codetools_0.2-18 np_0.60-11 ##  cachem_1.0.6 knitr_1.36 Formula_1.2-4 ##  jsonlite_1.7.2 drtmle_184.108.40.20600 cluster_2.1.2 ##  png_0.1-7 BiocManager_1.30.16 compiler_4.1.1 ##  backports_1.2.1 assertthat_0.2.1 Matrix_1.3-4 ##  fastmap_1.1.0 limma_3.48.3 htmltools_0.5.2 ##  quantreg_5.86 tools_4.1.1 coda_0.19-4 ##  gtable_0.3.0 glue_1.4.2 GenomeInfoDbData_1.2.6 ##  Rcpp_1.0.7 jquerylib_0.1.4 pkgdown_1.6.1 ##  vctrs_0.3.8 nlme_3.1-153 conquer_1.0.2 ##  iterators_1.0.13 xfun_0.26 stringr_1.4.0 ##  globals_0.14.0 lme4_1.1-27.1 lifecycle_1.0.1 ##  future_1.22.1 MASS_7.3-54 zlibbioc_1.38.0 ##  scales_1.1.1 ragg_1.1.3 SparseM_1.81 ##  RColorBrewer_1.1-2 yaml_2.2.1 gridExtra_2.3 ##  memoise_2.0.0 ggplot2_3.3.5 sass_0.4.0 ##  rpart_4.1-15 latticeExtra_0.6-29 stringi_1.7.5 ##  highr_0.9 desc_1.4.0 checkmate_2.0.0 ##  boot_1.3-28 rlang_0.4.11 pkgconfig_2.0.3 ##  systemfonts_1.0.2 bitops_1.0-7 superheat_0.1.0 ##  arm_1.11-2 evaluate_0.14 lattice_0.20-45 ##  purrr_0.3.4 labeling_0.4.2 htmlwidgets_1.5.4 ##  tidyselect_1.1.1 parallelly_1.28.1 ggsci_2.9 ##  magrittr_2.0.1 bookdown_0.24 R6_2.5.1 ##  generics_0.1.0 Hmisc_4.6-0 DelayedArray_0.18.0 ##  DBI_1.1.1 pillar_1.6.3 foreign_0.8-81 ##  abind_1.4-5 survival_3.2-13 RCurl_1.98-1.5 ##  nnet_7.3-16 tibble_3.1.5 future.apply_1.8.1 ##  crayon_1.4.1 utf8_1.2.2 rmarkdown_2.11 ##  jpeg_0.1-9 grid_4.1.1 data.table_1.14.2 ##  digest_0.6.28 textshaping_0.3.5 munsell_0.5.0 ##  bslib_0.3.1
Pearl, Judea. 2000. Causality: Models, Reasoning, and Inference. Cambridge University Press.
Pollard, Katherine S, and Mark J van der Laan. 2008. “Supervised Distance Matrices.” Statistical Applications in Genetics and Molecular Biology 7 (1).
Smyth, Gordon K. 2004. “Linear Models and Empirical Bayes Methods for Assessing Differential Expression in Microarray Experiments.” Statistical Applications in Genetics and Molecular Biology 3 (1): 1–25. https://doi.org/10.2202/1544-6115.1027.
van der Laan, Mark J., and Sherri Rose. 2011. Targeted Learning: Causal Inference for Observational and Experimental Data. Springer Science & Business Media.
van der Laan, Mark J, and Sherri Rose. 2018. Targeted Learning in Data Science: Causal Inference for Complex Longitudinal Studies. Springer Science & Business Media.