Contents

1 Introduction

omicRexposome is an R package designed to work join with rexposome. The aim of omicRexposome is to perform analysis joining exposome and omic data with the goal to find the relationship between a single or set of exposures (external exposome) and the behavior of a gene, a group of CpGs, the level of a protein, etc. Also to provide a series of tools to analyse exposome and omic data using standard methods from Biocondcutor.

1.1 Installation

omicRexposome is currently in development and not available from CRAN nor Bioconductor. Anyway, the package can be installed by using devtools R package and taking the source from Bioinformatic Research Group in Epidemiology’s GitHub repository.

This can be done by opening an R session and typing the following code:

devtools::install_github("isglobal-brge/omicRexposome")

User must take into account that this sentence do not install the packages’ dependencies.

1.2 Pipeline

Two different types of analyses can be done with omicRexposome:

Analysis omicRexposome function
Association Study association
Integration Study crossomics

Both association and integration studies are based in objects of class MultiDataSet. A MultiDataSet object is a contained for multiple layers of sample information. Once the exposome data and the omics data are encapsulated in a MultiDataSet the object can be used for both association and integration studies.

The method association requires a MultiDataSet object having to types of information: the exposome data from an ExposomeSet object and omic information from objects of class ExpressionSet, MethylationSet, SummarizedExperiment or others. ExposomeSet objects are created with functions read_exposome and load_exposome from rexposome R package (see next section Loading Exposome Data) and encapsulates exposome data. The method crossomics expects a MultiDataSet with any number of different data-sets (at last two). Compared with association method, crossomics do not requires an ExposomeSet.

1.3 Exposome and Omic Data

In order to illustrate the capabilities of omicRexposome and the exposome-omic analysis pipeline, we will use the data from BRGdata package. This package includes different omic-sets including methylation, transcriptome and proteome data-sets and an exposome-data set.

2 Analysis

omicRexposome and MultiDataSet R packages are loaded using the standard library command:

library(omicRexposome)
library(MultiDataSet)

2.1 Association Studies

The association studies are performed using the method association. This method requires, at last four, augments:

  1. Argument object should be filled with a MultiDataSet object.
  2. Argument formula should be filled with an expression containing the covariates used to adjust the model.
  3. Argument expset should be filled with the name that the exposome-set receives in the MultiDataSet object.
  4. Argument omicset should be filled with the name that the omic-set receives in the MultiDataSet object.

The argument formula should follow the pattern: ~sex+age. The method association will fill the formula placing the exposures in the ExposomeSetm between ~ and the covariates sex+age.

association implements the limma pipeline using lmFit and eBayes in the extraction methods from MultiDataSet. The method takes care of the missing data in exposures, outcomes and omics data and locating and is subsets both data-sets, exposome data and omic data, by common samples. The argument method allows to select the fitting method used in lmFit. By default it takes the value "ls" for least squares but it can also takes "robust" for robust regression.

The following subsections illustrates the usage of association with different types of omics data: methylome, transcriptome and proteome.

2.1.1 Exposome - Transcriptome Data Association

First we get the exposome data from BRGdata package that we will use in the whole section.

data("brge_expo", package = "brgedata")
class(brge_expo)
## [1] "ExposomeSet"
## attr(,"package")
## [1] "rexposome"

The aim of this analysis is to perform an association test between the gene expression levels and the exposures. So the first point is to obtain the transcriptome data from the brgedata package.

data("brge_gexp", package = "brgedata")

The association studies between exposures and transcriptome are done in the same way that the ones with methylome. The method used is association, that takes as input an object of MultiDataSet class with both exposome and expression data.

mds <- createMultiDataSet()
mds <- add_genexp(mds, brge_gexp)
mds <- add_exp(mds, brge_expo)

gexp <- association(mds, formula=~Sex+Age, 
    expset = "exposures", omicset = "expression")

We can have a look to the number of hits and the lambda score of each analysis with the methods tableHits and tableLambda, seen in the previous section.

hit <- tableHits(gexp, th=0.001)
lab <- tableLambda(gexp)
merge(hit, lab, by="exposure")
##    exposure hits    lambda
## 1     BPA_p   19 0.9072377
## 2    BPA_t1   27 0.8807316
## 3    BPA_t3   56 0.9391129
## 4     Ben_p   19 0.8013466
## 5    Ben_t1   12 0.8234104
## 6    Ben_t2    9 0.8393350
## 7    Ben_t3   21 0.8301203
## 8     NO2_p   32 1.0281960
## 9    NO2_t1   16 0.7942881
## 10   NO2_t2   35 1.1482314
## 11   NO2_t3   31 0.8770931
## 12   PCB118   59 0.9308472
## 13   PCB138   38 1.0726221
## 14   PCB153   51 1.1743989
## 15   PCB180   17 0.9790750

Since most of all models have a lambda under one, we should consider use Surrogate Variable Analysis. This can be done using the same association method but by setting the argument sva to "fast" so the pipeline of isva and SmartSVA R packages is applied. If sva is set to "slow" the applied. pipeline is the one from sva R package.

gexp <- association(mds, formula=~Sex+Age, 
    expset = "exposures", omicset = "expression", sva = "fast")

We can re-check the results creating the same table than before:

hit <- tableHits(gexp, th=0.001)
lab <- tableLambda(gexp)
merge(hit, lab, by="exposure")
##    exposure hits    lambda
## 1     BPA_p   50 0.9874152
## 2    BPA_t1   51 0.9453795
## 3    BPA_t3   61 0.9842216
## 4     Ben_p   76 1.0117733
## 5    Ben_t1   64 1.0115515
## 6    Ben_t2   71 1.0089834
## 7    Ben_t3   59 0.9969123
## 8     NO2_p   78 1.0117151
## 9    NO2_t1   68 1.0056950
## 10   NO2_t2   69 1.0210000
## 11   NO2_t3   49 0.9802407
## 12   PCB118  129 1.0518170
## 13   PCB138   67 1.0094139
## 14   PCB153   58 0.9924482
## 15   PCB180   67 0.9973381

The objects of class ResultSet have a method called plotAssociation that allows to create QQ Plots (that are another useful way to see if there are some inflation/deflation in the P-Values).

gridExtra::grid.arrange(
    plotAssociation(gexp, rid="Ben_p", type="qq") + 
        ggplot2::ggtitle("Transcriptome - Pb Association"),
    plotAssociation(gexp, rid="BPA_p", type="qq") + 
        ggplot2::ggtitle("Transcriptome - THM Association"),
    ncol=2
)

Following this line, the same method plotAssociation can be used to create volcano plots.

gridExtra::grid.arrange(
    plotAssociation(gexp, rid="Ben_p", type="volcano", tPV=-log10(1e-04)) + 
        ggplot2::ggtitle("Transcriptome - Pb Association"),
    plotAssociation(gexp, rid="BPA_p", type="volcano", tPV=-log10(1e-04)) + 
        ggplot2::ggtitle("Transcriptome - THM Association"),
    ncol=2
)

2.1.2 Exposome - Proteome Data Association

The proteome data-set included in brgedata has 47 proteins for 90 samples.

data("brge_prot", package="brgedata")
brge_prot
## ExpressionSet (storageMode: lockedEnvironment)
## assayData: 47 features, 90 samples 
##   element names: exprs 
## protocolData: none
## phenoData
##   sampleNames: x0001 x0002 ... x0090 (90 total)
##   varLabels: age sex
##   varMetadata: labelDescription
## featureData
##   featureNames: Adiponectin_ok Alpha1AntitrypsinAAT_ok ...
##     VitaminDBindingProte_ok (47 total)
##   fvarLabels: chr start end
##   fvarMetadata: labelDescription
## experimentData: use 'experimentData(object)'
## Annotation:

The association analysis between exposures and proteome is also done using association.

mds <- createMultiDataSet()
mds <- add_eset(mds, brge_prot, dataset.type  ="proteome")
mds <- add_exp(mds, brge_expo)

prot <- association(mds, formula=~Sex+Age,
    expset = "exposures", omicset = "proteome")

The tableHits indicates that no association was found between the 47 proteins and the exposures.

tableHits(prot, th=0.001)
##        exposure hits
## Ben_p     Ben_p    0
## Ben_t1   Ben_t1    0
## Ben_t2   Ben_t2    0
## Ben_t3   Ben_t3    0
## BPA_p     BPA_p    0
## BPA_t1   BPA_t1    0
## BPA_t3   BPA_t3    0
## NO2_p     NO2_p    0
## NO2_t1   NO2_t1    1
## NO2_t2   NO2_t2    0
## NO2_t3   NO2_t3    0
## PCB118   PCB118    0
## PCB138   PCB138    0
## PCB153   PCB153    0
## PCB180   PCB180    0

This is also seen in the Manhattan plot for proteins that can be obtained from plotAssociation.

gridExtra::grid.arrange(
    plotAssociation(prot, rid="Ben_p", type="protein") + 
        ggplot2::ggtitle("Proteome - Cd Association") +
        ggplot2::geom_hline(yintercept = 1, color = "LightPink"),
    plotAssociation(prot, rid="NO2_p", type="protein") + 
        ggplot2::ggtitle("Proteome - Cotinine Association") +
        ggplot2::geom_hline(yintercept = 1, color = "LightPink"),
    ncol=2
)

NOTE: A real Manhattan plot can be draw with plot method for ResultSet objects by setting the argument type to "manhattan".

##            used  (Mb) gc trigger   (Mb)  max used   (Mb)
## Ncells  6981140 372.9   12527447  669.1  12527447  669.1
## Vcells 19537041 149.1  360795230 2752.7 450988177 3440.8

2.2 Integration Analysis

omicRexposome allows to study the relation between exposures and omic-features from another perspective, different from the association analyses. The integration analysis can be done, in omicRexposome using multi canonical correlation analysis or using multiple co-inertia analysis. The first methods is implemented in R package PMA (CRAN) and the second in omicade4 R package (Bioconductor). The two methods are encapsulated in the crossomics method.

The differences between association and crossomics are that the first method test association between two complete data-sets, by removing the samples having missing values in any of the involved data-sets, and the second try to find latent relationships between two or more sets.

Hence, we need to explore the missing data in the exposome data-set. This can be done using the methods plotMissings and tableMissings from rexposome R package.

library(rexposome)
plotMissings(brge_expo, set = "exposures")

From the plot we can see that more of the exposures have up to 25% of missing values. Hence the first step in the integration analysis is to avoid missing values. so, we perform a fast imputation on the exposures side:

brge_expo <- imputation(brge_expo)

crossomics function expects to obtain the different data-sets in a single labelled-list, in the argument called list. The argument method from crossomics function can be set to mcia (for multiple co-inertia analysis) or to mcca (for multi canonical correlation analysis).

The following code shows how to perform the integration of the exposome and the proteome. The method crossomics request a MultiDataSet object as input, containing the data-set to be integrated.

mds <- createMultiDataSet()
mds <- add_genexp(mds, brge_gexp)
mds <- add_eset(mds, brge_prot, dataset.type = "proteome")
mds <- add_exp(mds, brge_expo)

cr_mcia <- crossomics(mds, method = "mcia", verbose = TRUE)
cr_mcia
## Object of class 'ResultSet'
##  . created with: crossomics 
##  . sva:   
##     . method: mcia  ( omicade4 )
##  . #results: 1 ( error: 0 )
##  . featureData:  3 
##     . expression: 67528x11
##     . proteome: 47x3
##     . exposures: 15x12

As can be seen, crossomics returns an object of class ResultSet. In the integration process, the different data-sets are subset by common samples. This is done taking advantage of MultiDataSet capabilities.

The same is done when method is set to mcca.

cr_mcca <- crossomics(mds, method = "mcca", permute=c(4, 2))
cr_mcca

We used an extra argument (permute) into the previous call to crossomics using multi canonical correlation analysis. This argument allows to set the internal argument corresponding to permutations and iterations, that are used to tune-up internal parameters.

When a ResultSet is generated using crossomics the methods plotHits, plotLambda and plotAssociation can NOT be used. But the plotIntegration will help us to understand what was done. This method allows to provide the colors to be used on the plots:

colors <- c("green", "blue", "red")
names(colors) <- names(mds)

The graphical representation of the results from a multiple co-inertia analysis is a composition of four different plots.

plotIntegration(cr_mcia, colors=colors)

The first plot (first row, first column) is the samples space. It illustrates how the different data-sets are related in terms of intra-sample variability (each data-set has a different color). The second plot (first row, second column) shows the feature space. The features of each set are drawn on the same components so the relation between each data-set can be seen (the features are colored depending of the set were they belong).

The third plot (second row, first column) shows the inertia of each component. The two first plots are drawn on the first and second component. Finally, the fourth plot shows the behavior of the data-sets.

A radar plots is obtained when plotIntegration is used on a ResultSet created though multi canonical correlation analysis.

plotIntegration(cr_mcca, colors=colors)

This plot shows the features of the three data-sets in the same 2D space.The relation between the features can be understood by proximity. This means that the features that clusters, or that are in the same quadrant are related and goes in a different direction than the features in the other quadrants.

rm(cr_mcia, cr_mcca)

Session info

## R version 4.0.3 (2020-10-10)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows Server 2012 R2 x64 (build 9600)
## 
## Matrix products: default
## 
## locale:
## [1] LC_COLLATE=C                          
## [2] LC_CTYPE=English_United States.1252   
## [3] LC_MONETARY=English_United States.1252
## [4] LC_NUMERIC=C                          
## [5] LC_TIME=English_United States.1252    
## 
## attached base packages:
## [1] parallel  stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
## [1] rexposome_1.12.1     MultiDataSet_1.18.0  omicRexposome_1.12.1
## [4] Biobase_2.50.0       BiocGenerics_0.36.0  BiocStyle_2.18.1    
## 
## loaded via a namespace (and not attached):
##   [1] backports_1.2.1             circlize_0.4.12            
##   [3] Hmisc_4.4-2                 corrplot_0.84              
##   [5] qqman_0.1.4                 plyr_1.8.6                 
##   [7] gmm_1.6-5                   splines_4.0.3              
##   [9] BiocParallel_1.24.1         omicade4_1.30.0            
##  [11] GenomeInfoDb_1.26.2         fastICA_1.2-2              
##  [13] ggplot2_3.3.3               made4_1.64.0               
##  [15] pryr_0.1.4                  sva_3.38.0                 
##  [17] digest_0.6.27               foreach_1.5.1              
##  [19] htmltools_0.5.1.1           magick_2.6.0               
##  [21] magrittr_2.0.1              checkmate_2.0.0            
##  [23] memoise_1.1.0               cluster_2.1.0              
##  [25] limma_3.46.0                annotate_1.68.0            
##  [27] matrixStats_0.57.0          isva_1.9                   
##  [29] imputeLCMD_2.0              sandwich_3.0-0             
##  [31] prettyunits_1.1.1           jpeg_0.1-8.1               
##  [33] colorspace_2.0-0            blob_1.2.1                 
##  [35] ggrepel_0.9.1               xfun_0.20                  
##  [37] dplyr_1.0.3                 crayon_1.3.4               
##  [39] RCurl_1.98-1.2              lme4_1.1-26                
##  [41] genefilter_1.72.1           impute_1.64.0              
##  [43] survival_3.2-7              zoo_1.8-8                  
##  [45] iterators_1.0.13            glue_1.4.2                 
##  [47] gtable_0.3.0                zlibbioc_1.36.0            
##  [49] XVector_0.30.0              DelayedArray_0.16.1        
##  [51] shape_1.4.5                 scales_1.1.1               
##  [53] mvtnorm_1.1-1               DBI_1.1.1                  
##  [55] edgeR_3.32.1                Rcpp_1.0.6                 
##  [57] SmartSVA_0.1.3              xtable_1.8-4               
##  [59] progress_1.2.2              htmlTable_2.1.0            
##  [61] clue_0.3-58                 flashClust_1.01-2          
##  [63] foreign_0.8-81              bit_4.0.4                  
##  [65] Formula_1.2-4               stats4_4.0.3               
##  [67] DT_0.17                     glmnet_4.1                 
##  [69] htmlwidgets_1.5.3           httr_1.4.2                 
##  [71] gplots_3.1.1                RColorBrewer_1.1-2         
##  [73] calibrate_1.7.7             ellipsis_0.3.1             
##  [75] farver_2.0.3                pkgconfig_2.0.3            
##  [77] XML_3.99-0.5                nnet_7.3-14                
##  [79] locfit_1.5-9.4              labeling_0.4.2             
##  [81] tidyselect_1.1.0            rlang_0.4.10               
##  [83] reshape2_1.4.4              AnnotationDbi_1.52.0       
##  [85] munsell_0.5.0               tools_4.0.3                
##  [87] generics_0.1.0              RSQLite_2.2.2              
##  [89] ade4_1.7-16                 evaluate_0.14              
##  [91] stringr_1.4.0               yaml_2.2.1                 
##  [93] knitr_1.30                  bit64_4.0.5                
##  [95] caTools_1.18.1              purrr_0.3.4                
##  [97] nlme_3.1-151                JADE_2.0-3                 
##  [99] leaps_3.1                   PMA_1.2.1                  
## [101] compiler_4.0.3              rstudioapi_0.13            
## [103] png_0.1-7                   statmod_1.4.35             
## [105] tibble_3.0.5                stringi_1.5.3              
## [107] RSpectra_0.16-0             lattice_0.20-41            
## [109] Matrix_1.3-2                nloptr_1.2.2.2             
## [111] tmvtnorm_1.4-10             vctrs_0.3.6                
## [113] norm_1.0-9.5                pillar_1.4.7               
## [115] lifecycle_0.2.0             BiocManager_1.30.10        
## [117] GlobalOptions_0.1.2         data.table_1.13.6          
## [119] bitops_1.0-6                GenomicRanges_1.42.0       
## [121] qvalue_2.22.0               pcaMethods_1.82.0          
## [123] R6_2.5.0                    latticeExtra_0.6-29        
## [125] bookdown_0.21               KernSmooth_2.23-18         
## [127] gridExtra_2.3               IRanges_2.24.1             
## [129] codetools_0.2-18            boot_1.3-25                
## [131] MASS_7.3-53                 gtools_3.8.2               
## [133] assertthat_0.2.1            SummarizedExperiment_1.20.0
## [135] S4Vectors_0.28.1            GenomeInfoDbData_1.2.4     
## [137] mgcv_1.8-33                 hms_1.0.0                  
## [139] grid_4.0.3                  rpart_4.1-15               
## [141] minqa_1.2.4                 rmarkdown_2.6              
## [143] MatrixGenerics_1.2.0        lsr_0.5                    
## [145] scatterplot3d_0.3-41        base64enc_0.1-3            
## [147] FactoMineR_2.4