An introduction to hRUV package
Taiyun Kim
School of Mathematics and Statistics, The University of SydneyhRUV.RmdIntroduction
hRUV is a package for normalisation of multiple batches of metabolomics data in a hierarchical strategy with use of samples replicates in a large-scale studies. The tool utilises 2 types of replicates: intra-batch and inter-batch replicates to estimate the unwanted variation within and between batches with RUV-III. We have designed the replicate embedding arrangements within and between batches from http://shiny.maths.usyd.edu.au/hRUV. Our novel tool is a novel hierarchical approach to removing unwanted variation by harnessing information from sample replicates embedded in the sequence of experimental runs/batches and applying signal drift correction with robust linear or non-linear smoothers.
Below is a schematic overview of hRUV framework.
Overview of hRUV methods
The purpose of this vignette is to illustrate use of hRUV.
Installation
Install the R package from GitHub using the devtools package:
if (!("devtools" %in% rownames(installed.packages())))
install.packages("devtools")
library(devtools)
devtools::install_github("SydneyBioX/hRUV", build_vignettes = TRUE)Loading packages and data
First, we will load the hRUV package and other packages required for the demonstration. If you haven’t installed the packages yet, please follow the steps above. You can install dplyr package with the command install.packages("dplyr") and SummarizedExperiment object
suppressPackageStartupMessages({
library(hRUV)
library(dplyr)
library(SummarizedExperiment)
})
# You can install dplyr package with the command
# install.packages("dplyr")
# To install SummarizedExperiment,
# if (!requireNamespace("BiocManager", quietly = TRUE))
# install.packages("BiocManager")
#
# BiocManager::install("SummarizedExperiment")For demonstration purposes, we provide a first 5 batches of BioHEART-CT metabolomics data. The data contains metabolite signals of each patients, and contains sample replicate within and between batches. The data is deposited at https://github.com/SydneyBioX/BioHEART_metabolomics.
## [1] "SummarizedExperiment"
## attr(,"package")
## [1] "SummarizedExperiment"
b1## class: SummarizedExperiment
## dim: 79 89
## metadata(1): name
## assays(1): raw
## rownames(79): 1-methylhistamine 2-Arachidonyl glycerol ... Valine
## Valine-d8
## rowData names(1): metabolite
## colnames(89): 1_Pool 0 2_26 ... 88_Pool 9 89_25
## colData names(7): Sample sample_name ... biological_sample batch_info
b2## class: SummarizedExperiment
## dim: 79 91
## metadata(1): name
## assays(1): raw
## rownames(79): 1-methylhistamine 2-Arachidonyl glycerol ... Valine
## Valine-d8
## rowData names(1): metabolite
## colnames(91): 1_Pool 1 2_115 ... 90_193 91_Pool 11
## colData names(7): Sample sample_name ... biological_sample batch_infoThe data is already formatted in to a SummarizedExperiment object.
Data cleaning
We will now perform preprocessing of all batches.
dat_list = list(
"Batch1" = b1,
"Batch2" = b2,
"Batch3" = b3,
"Batch4" = b4,
"Batch5" = b5
)
# checking the dimensions of the data
lapply(dat_list, dim)## $Batch1
## [1] 79 89
##
## $Batch2
## [1] 79 91
##
## $Batch3
## [1] 79 87
##
## $Batch4
## [1] 79 95
##
## $Batch5
## [1] 79 90
# Log transform raw assay
dat_list = lapply(dat_list, function(dat) {
assay(dat, "logRaw", withDimnames = FALSE) = log2(assay(dat, "raw") + 1)
dat
})
# Setting the order of batches
h_order = paste("Batch", 1:5, sep = "")
h_order = h_order[order(h_order)]
# checking the dimensions of the data before cleaning
lapply(dat_list, dim)## $Batch1
## [1] 79 89
##
## $Batch2
## [1] 79 91
##
## $Batch3
## [1] 79 87
##
## $Batch4
## [1] 79 95
##
## $Batch5
## [1] 79 90
dat_list = hRUV::clean(dat_list, threshold = 0.5,
method = "intersect", assay = "logRaw", newAssay = "rawImpute")
# checking the dimensions of the data after cleaning
lapply(dat_list, dim)## $Batch1
## [1] 61 89
##
## $Batch2
## [1] 61 91
##
## $Batch3
## [1] 61 87
##
## $Batch4
## [1] 61 95
##
## $Batch5
## [1] 61 90We have filtered metabolites with more than 50% of missing values per batch and selected metabolites that are quantified across all batches (intersect) with clean function. After filtering process, clean function performs k-nearest neighbour imputation. The new imputed assay is named as rawImpute. After cleaning we now have 61 intersecting metabolites quantified for each batch.
Perfrom hRUV
We can proceed to perform intra- and inter-batch normalisation with hRUV.
dat = hruv(
dat_list = dat_list,
assay = "rawImpute",
intra = "loessShort",
inter = "concatenate",
intra_k = 5, inter_k = 5,
pCtlName = "biological_sample",
negCtl = NULL,
intra_rep = "short_replicate",
inter_rep = "batch_replicate"
)Here, we perform hRUV on rawImpute assay. For intra batch normalisation, we perform loess smoothing on samples and RUV-III using short replicates with parameter k set to 5. For inter batch normalisation, we perform concatenating hierarchical structure using batch replicate samples. The RUV-III paramter k is also set to 5 for inter-batch normalisation.
We can now compare before and after hruv normalisation.
PCA
library(gridExtra)##
## Attaching package: 'gridExtra'## The following object is masked from 'package:Biobase':
##
## combine## The following object is masked from 'package:BiocGenerics':
##
## combine## The following object is masked from 'package:dplyr':
##
## combine
p1 = hRUV::plotPCA(dat, assay = "rawImpute", colour = dat$batch_info)
p2 = hRUV::plotPCA(dat, assay = "loessShort_concatenate", colour = dat$batch_info)
grid.arrange(p1, p2, nrow = 1)
From the above PCA plot, you can see that the assay before normalisation has strong batch effect. The assay after normalisation on the right panel does no longer show batch effect.
Run plot
p1 = hRUV::plotRun(dat, assay = "rawImpute", colour = dat$batch_info)
p2 = hRUV::plotRun(dat, assay = "loessShort_concatenate", colour = dat$batch_info)
grid.arrange(p1$`1-methylhistamine`, p2$`1-methylhistamine`, nrow = 2)
The run plot above shows signal drift present in the raw data for metabolite 1-methylhistamine. This is corrected after hRUV normalisation.
grid.arrange(p1$GlucosePos2, p2$GlucosePos2, nrow = 2)
Further, in metabolite GlucosePos2, the assay before normalisation shows strong batch effects, but they are corrected after hRUV normalisation.
SessionInfo
## R version 4.2.1 (2022-06-23)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Big Sur ... 10.16
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] stats4 stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] gridExtra_2.3 SummarizedExperiment_1.26.1
## [3] Biobase_2.56.0 GenomicRanges_1.48.0
## [5] GenomeInfoDb_1.32.4 IRanges_2.30.1
## [7] S4Vectors_0.34.0 BiocGenerics_0.42.0
## [9] MatrixGenerics_1.8.1 matrixStats_0.62.0
## [11] dplyr_1.0.10 hRUV_0.1.3
## [13] BiocStyle_2.24.0
##
## loaded via a namespace (and not attached):
## [1] tidyr_1.2.1 sass_0.4.2 jsonlite_1.8.2
## [4] bslib_0.4.0 TTR_0.24.3 highr_0.9
## [7] BiocManager_1.30.18 GenomeInfoDbData_1.2.8 impute_1.70.0
## [10] yaml_2.3.5 pillar_1.8.1 lattice_0.20-45
## [13] glue_1.6.2 digest_0.6.29 XVector_0.36.0
## [16] colorspace_2.0-3 htmltools_0.5.3 Matrix_1.4-1
## [19] pkgconfig_2.0.3 bookdown_0.29 zlibbioc_1.42.0
## [22] purrr_0.3.5 scales_1.2.1 tzdb_0.3.0
## [25] tibble_3.1.8 farver_2.1.1 generics_0.1.3
## [28] ggplot2_3.3.6 ellipsis_0.3.2 withr_2.5.0
## [31] cachem_1.0.6 cli_3.4.1 quantmod_0.4.20
## [34] DMwR2_0.0.2 magrittr_2.0.3 memoise_2.0.1
## [37] evaluate_0.17 fs_1.5.2 fansi_1.0.3
## [40] MASS_7.3-57 xts_0.12.2 class_7.3-20
## [43] textshaping_0.3.6 tools_4.2.1 hms_1.1.2
## [46] lifecycle_1.0.3 stringr_1.4.1 munsell_0.5.0
## [49] DelayedArray_0.22.0 compiler_4.2.1 pkgdown_2.0.6
## [52] jquerylib_0.1.4 systemfonts_1.0.4 rlang_1.0.6
## [55] grid_4.2.1 RCurl_1.98-1.9 labeling_0.4.2
## [58] bitops_1.0-7 rmarkdown_2.17 gtable_0.3.1
## [61] DBI_1.1.3 curl_4.3.3 R6_2.5.1
## [64] zoo_1.8-11 knitr_1.40 fastmap_1.1.0
## [67] utf8_1.2.2 rprojroot_2.0.3 ragg_1.2.3
## [70] readr_2.1.3 desc_1.4.2 stringi_1.7.8
## [73] vctrs_0.4.2 rpart_4.1.16 tidyselect_1.2.0
## [76] xfun_0.33