1 Introduction

High-throughput technologies have revolutionized the genomics research. The early applications of the technologies were largely on cell lines. However, there is an increasing number of larger-scale, population level clinical studies in recent years, hoping to identify diagnostic biomarkers and therapeutic targets. The samples collected in these studies, such as blood, tumor, or brain tissue, are mixtures of a number of different cell types. The sample mixing complicates data analysis because the experimental data from the high-throughput experiments are weighted average of signals from multiple cell types. For these data, traditional analysis methods that ignores the cell mixture will lead to results with low resolution, biased, or even errorneous results. For example, it has been discovered that in epigenome-wide association studies (EWAS), the mixing proportions can be confounded with the experimental factor of interest (such as age). Ignoring the cell mixing will lead to false positives. On the other hand, cell type specific changes under different conditions could be associated with disease pathogenesis and progressions, which are of great interests to researchers.

For heterogeneous samples, it is possible to profile the pure cell types through experimental techniques. They are, however, laborious and expensive that cannot be applied to large scale studies. Computational tools for analzying the mixed data have been developed for proportion estimation and cell type specific signal detection. There are two fundamental questions in this type of analyses:

  1. How to estimate mixing proportions?

There are a number of existing methods devoted to solve this question. These methods mainly can be categorized to two groups: reference-based (require pure cell type profiles) and reference-free (does not require pure cell type profiles). It has been found that reference-based deconvolution is more accurate and reliable than reference-free deconvolution. However, the reference panels required for reference-based deconvolution can be difficult to obtain, thus reference-free method has wider application.

  1. with available mixing proportions, how to detect cell-type specific DE/DM?

TOAST is a package designed to answer these questions and serve the research communities with tools for the analysis of heterogenuous tissues. Currently TOAST provides functions to detect cell-type specific DE/DM, as well as differences across different cell types. TOAST also has functions to improve the accuracy of reference-free deconvolutions through better feature selection.

2 Installation and quick start

2.1 Install TOAST

To install this package, start R (version “3.6”) and enter:

if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

BiocManager::install("TOAST")

2.2 How to get help for TOAST

Any TOAST questions should be posted to the GitHub Issue section of TOAST homepage at https://github.com/ziyili20/TOAST/issues.

2.3 Quick start on detecting cell type-specific differential signals

Here we show the key steps for a cell type-specific different analysis. This code chunk assumes you have an expression or DNA methylation matrix called Y_raw, a data frame of sample information called design, and a table of cellular composition (i.e. mixing proportions) called prop. Instead of a data matrix, Y_raw could also be a SummarizedExperiment object. If the cellular composition is not available, the following sections will discuss about how to obtain mixing proportions using reference-free deconvolution or reference-based deconvolution.

Design_out <- makeDesign(design, Prop)
fitted_model <- fitModel(Design_out, Y_raw)
fitted_model$all_coefs # list all phenotype names
fitted_model$all_cell_types # list all cell type names
# coef should be one of above listed phenotypes
# cell_type should be one of above listed cell types
res_table <- csTest(fitted_model, coef = "age", 
                    cell_type = "Neuron", contrast_matrix = NULL)
head(res_table)

3 Example dataset

TOAST provides an example dataset of 450K DNA methylation. The following sections will use this example dataset to demonstrate function usages.

We obtain and process this dataset based on the raw data provided by GSE42861. This is a DNA methylation 450K data for Rheumatoid Arthiritis patients and controls. The original dataset has 485577 features and 689 samples. We have reduced the dataset to 3000 CpGs for randomly selected 50 RA patients and 50 controls.

library(TOAST)
## Loading required package: RefFreeEWAS
## Loading required package: quadprog
## Loading required package: EpiDISH
data("RA_100samples")
Y_raw <- RA_100samples$Y_raw
Pheno <- RA_100samples$Pheno
Blood_ref <- RA_100samples$Blood_ref

Check matrix including beta values for 3000 CpG by 100 samples.

dim(Y_raw) 
## [1] 3000  100
Y_raw[1:4,1:4]
##            GSM1051525 GSM1051526 GSM1051527  GSM1051528
## cg14521995  0.8848926  0.8654487  0.8172092 0.004429362
## cg11738485  0.9306579  0.9189274  0.5486962 0.039545301
## cg06193597  0.1388632  0.7127654  0.6925506 0.677185017
## cg14323910  0.8282483  0.8528023  0.8449638 0.828873689

Check phenotype of these 100 samples.

dim(Pheno)
## [1] 100   3
head(Pheno, 3)
##            age gender disease
## GSM1051525  67      2       1
## GSM1051526  49      2       1
## GSM1051527  53      2       1

Our example dataset also contain blood reference matrix for the matched 3000 CpGs (obtained from bioconductor package FlowSorted.Blood.450k.

dim(Blood_ref)
## [1] 3000    6
head(Blood_ref, 3)
##                 CD8T      CD4T        NK     Bcell      Mono      Gran
## cg14521995 0.9321049 0.9245206 0.9184654 0.9178081 0.8902820 0.9314544
## cg11738485 0.3745548 0.2916655 0.3144788 0.2985633 0.3027369 0.2911957
## cg06193597 0.4157621 0.4292540 0.4104737 0.4335429 0.4789953 0.4747334

4 Estimate mixing proportions

If you have mixing proportions available, you can directly go to Section 5.

In many situations, mixing proportions are not readily available. There are a number of deconvolution methods available to solve this problem. To name a few:

  • For DNA methylation: RefFreeEWAS (Houseman et al. 2016) is reference-free,
    and EpiDISH (Teschendorff et al. 2017) is reference-based.

  • For gene expression: qprog (Gong et al. 2011), deconf (Repsilber et al. 2010),
    lsfit (Abbas et al. 2009) and DSA (Zhong et al. 2013).

In addition, CellMix package has summarized a number of deconvolution methods and is a good resource to look up.

Here we demonstrate two ways to estimate mixing proportions, one using RefFreeEWAS (Houseman et al. 2016), representing the class of reference-free methods, and the other using EpiDISH (Teschendorff et al. 2017) as a representation of reference-based methods.

We also provide function to improve reference-free deconvolution performance in Section 4.3. Note that in the following example, we have only 3000 features in the Y_raw dataset, thus the proportion estimation is not very accurate. Real 450K dataset should have around 485,000 features. More features generally lead to better estimation, because there are more information in the data.

4.1 Reference-based deconvolution using least square method

  1. Select the top 1000 most variant features by findRefinx(). To select the top features with largest coefficients of variations, one can use findRefinx(..., sortBy = "cv"). Default sortBy argument is "var". Here, instead of a data matrix, Y_raw could also be a SummarizedExperiment object.
refinx <- findRefinx(Y_raw, nmarker = 1000)
  1. Subset data and reference panel.
Y <- Y_raw[refinx,]
Ref <- as.matrix(Blood_ref[refinx,])
  1. Use EpiDISH to solve cellular proportions and use post-hoc constraint.
library(EpiDISH)
outT <- epidish(beta.m = Y, ref.m = Ref, method = "RPC")
estProp_RB <- outT$estF

A word about Step 1

For step 1, one can also use findRefinx(..., sortBy = "cv") to select features based on coefficient of variantion. The choice of sortby = "cv" and sortBy = "var" depends on whether the feature variances of your data correlates with the means. For RNA-seq counts, the variance-mean correlation is strong, thus sortBy = "cv" is recommended. For log-counts, the variance-mean correlation largely disappears, so both sortBy = "cv" and sortBy = "var" would work similarly. In DNA methylation data, this correlation is not strong, either sortBy = "cv" or sortBy = "var" can be used. In this case, we recommend sortBy = "var" because we find it has better feature selection for DNA methylation data than sortBy = "cv" (unpublished results).

refinx = findRefinx(Y_raw, nmarker=1000, sortBy = "var")

4.2 Reference-free deconvolution using RefFreeEWAS

  1. Load RefFreeEWAS.
library(RefFreeEWAS)
  1. Similar to Reference-based deconvolution we also select the top 1000 most variant features by findRefinx(). And then subset data.
refinx <- findRefinx(Y_raw, nmarker = 1000)
Y <- Y_raw[refinx,]
  1. Do reference-free deconvolution on the RA dataset.
K <- 6
outT <- RefFreeCellMix(Y, mu0=RefFreeCellMixInitialize(Y, K = K))
estProp_RF <- outT$Omega
  1. Comparing the reference-free method versus reference-base method
# first we align the cell types from RF 
# and RB estimations using pearson's correlation
estProp_RF <- assignCellType(input=estProp_RF,
                             reference=estProp_RB) 
mean(diag(cor(estProp_RF, estProp_RB)))
## [1] 0.1967946

4.3 Improve reference-free deconvolution with cross-cell type differential analysis

Feature selection is an important step before RF deconvolution and is directly related with the estimation quality of cell composition. findRefinx() and findRefinx(..., sortBy = "var") simply select the markers with largest CV or largest variance, which may not always result in a good selection of markers. Here, we propose to improve RF deconvolution marker selection through cross-cell type differential analysis. We implement two versions of such improvement, one is for DNA methylation microarray data using RefFreeCellMix from RefFreeEWAS package, the other one is for gene expression microarray data using deconf from CellMix package. To implement both, RefFreeEWAS and CellMix need to be installed first.

4.3.1 Improved-RF with RefFreeCellMix

  1. Load TOAST package.
library(TOAST)
  1. Do reference-free deconvolution using improved-RF implemented with RefFreeCellMix. The default deconvolution function implemented in csDeconv() is RefFreeCellMix_wrapper(). Here, instead of a data matrix, Y_raw could also be a SummarizedExperiment object.
K=6
set.seed(1234)
outRF1 <- csDeconv(Y_raw, K, TotalIter = 30, bound_negative = TRUE)
  1. Comparing udpated RF estimations versus RB results.
## check the accuracy of deconvolution
estProp_RF_improved <- assignCellType(input=outRF1$estProp,
                                      reference=estProp_RB) 
mean(diag(cor(estProp_RF_improved, estProp_RB)))
## [1] 0.2254084

A word about Step 2

For step 2, initial features (instead of automatic selection by largest variation) can be provided to function RefFreeCellMixT(). For example

refinx <- findRefinx(Y_raw, nmarker = 1000, sortBy = "cv")
InitNames <- rownames(Y_raw)[refinx]
csDeconv(Y_raw, K = 6, nMarker = 1000, 
         InitMarker = InitNames, TotalIter = 30)

A word about bounding the negative estimators

Since all the parameters represent the mean observation levels for each cell type, it may not be reasonable to have negative estimators. As such, we provide options to bound negative estimated parameters to zero through the bound_negative argument in csDeconv() function. Although we find bounding negative estimators has minimum impact on the performance, the users could choose to bound or not bound the negative values in the function. The default value for bound_negative is FALSE.

4.3.2 Improved-RF with use-defined RF function

In order to use other RF functions, users can wrap the RF function a bit first to make it accept Y (raw data) and K (number of cell types) as input, and return a N (number of cell types) by K proportion matrix. We take RefFreeCellMix() as an example. Other deconvolution methods can be used similarly.

mydeconv <- function(Y, K){
     if (is(Y, "SummarizedExperiment")) {
          se <- Y
          Y <- assays(se)$counts
     } else if (!is(Y, "matrix")) {
          stop("Y should be a matrix
               or a SummarizedExperiment object!")
     }
     
     if (K<0 | K>ncol(Y)) {
         stop("K should be between 0 and N (samples)!")
     }
     outY = RefFreeEWAS::RefFreeCellMix(Y, 
               mu0=RefFreeEWAS::RefFreeCellMixInitialize(Y, 
               K = K))
     Prop0 = outY$Omega
     return(Prop0)
}
set.seed(1234)
outT <- csDeconv(Y_raw, K, FUN = mydeconv, bound_negative = TRUE)

5 Detect cell type-specific and cross-cell type differential signals

The csDE/csDM detection function requires a table of microarray or RNA-seq measurements from all samples, a table of mixing proportions, and a design vector representing the status of subjects.

We demonstrate the usage of TOAST in three common settings.

5.1 Detect cell type-specific differential signals under two-group comparison

  1. Assuming you have TOAST library and dataset loaded, the first step is to generate the study design based on the phenotype matrix. Note that all the binary (e.g. disease = 0, 1) or categorical variable (e.g. gender = 1, 2) should be transformed to factor class. Here we use the proportions estimated from step 4.3.1 as input proportion.
head(Pheno, 3)
##            age gender disease
## GSM1051525  67      2       1
## GSM1051526  49      2       1
## GSM1051527  53      2       1
design <- data.frame(disease = as.factor(Pheno$disease))

Prop <- estProp_RF_improved
colnames(Prop) <- colnames(Ref) 
## columns of proportion matrix should have names
  1. Make model design using the design (phenotype) data frame and proportion matrix.
Design_out <- makeDesign(design, Prop)
  1. Fit linear models for raw data and the design generated from Design_out(). Y_raw here is a data matrix with dimension P (features) by N (samples). Instead of a data matrix, Y_raw could also be a SummarizedExperiment object.
fitted_model <- fitModel(Design_out, Y_raw)
# print all the cell type names
fitted_model$all_cell_types
## [1] "CD8T"  "CD4T"  "NK"    "Bcell" "Mono"  "Gran"
# print all phenotypes
fitted_model$all_coefs
## [1] "disease"

TOAST allows a number of hypotheses to be tested using csTest() in two group setting.

5.1.1 Testing one parameter (e.g. disease) in one cell type.

For example, testing disease (patient versus controls) effect in Gran.

res_table <- csTest(fitted_model, 
                    coef = "disease", 
                    cell_type = "Gran")
## Test the effect of disease in Gran.
head(res_table, 3)
##                  beta   beta_var         mu effect_size f_statistics
## cg03999583 -0.9336162 0.03923521  1.2339749   -1.216966     22.21573
## cg04021544 -0.6798894 0.02899288  0.8520754   -1.327570     15.94356
## cg07755735  0.7132178 0.03332731 -0.1296480    3.142470     15.26315
##                 p_value        fdr
## cg03999583 9.065204e-06 0.02719561
## cg04021544 1.349669e-04 0.15106724
## cg07755735 1.831638e-04 0.15106724
Disease_Gran_res <- res_table

5.1.2 Testing one parameter in all cell types.

For example, testing the joint effect of age in all cell types:

res_table <- csTest(fitted_model, 
                    coef = "disease", 
                    cell_type = "joint")
head(res_table, 3)

Specifying cell_type as NULL or not specifying cell_type will test the effect in each cell type and the joint effect in all cell types.

res_table <- csTest(fitted_model, 
                    coef = "disease", 
                    cell_type = NULL)
lapply(res_table, head, 3)

## this is exactly the same as
res_table <- csTest(fitted_model, coef = "disease")

5.2 Detect cell type-specific differential signals from a general experimental design

  1. Assuming you have TOAST library and dataset loaded, generate the study design based on the phenotype matrix. Note that all the binary variable (e.g. disease = 0, 1) or categorical variable (e.g. gender = 1, 2) should be transformed to factor class.
design <- data.frame(age = Pheno$age,
                     gender = as.factor(Pheno$gender),
                     disease = as.factor(Pheno$disease))

Prop <- estProp_RF_improved
colnames(Prop) <- colnames(Ref)  
## columns of proportion matrix should have names
  1. Make model design using the design (phenotype) data frame and proportion matrix.
Design_out <- makeDesign(design, Prop)
  1. Fit linear models for raw data and the design generated from Design_out().
fitted_model <- fitModel(Design_out, Y_raw)
# print all the cell type names
fitted_model$all_cell_types
## [1] "CD8T"  "CD4T"  "NK"    "Bcell" "Mono"  "Gran"
# print all phenotypes
fitted_model$all_coefs
## [1] "age"     "gender"  "disease"

TOAST allows a number of hypotheses to be tested using csTest() in two group setting.

5.2.1 Testing one parameter in one cell type

For example, testing age effect in Gran.

res_table <- csTest(fitted_model, 
                    coef = "age", 
                    cell_type = "Gran")
## Test the effect of age in Gran.
head(res_table, 3)
##                   beta     beta_var       mu effect_size f_statistics
## cg10785373 -0.03998455 8.795653e-05 2.502370 -0.01610736     18.17675
## cg05364038 -0.03051950 5.682604e-05 2.210576 -0.01390210     16.39108
## cg16611967 -0.02758025 5.284133e-05 1.658134 -0.01677280     14.39536
##                 p_value       fdr
## cg10785373 0.0000572205 0.1716615
## cg05364038 0.0001229534 0.1844301
## cg16611967 0.0002956521 0.2956521

We can test disease effect in Bcell.

res_table <- csTest(fitted_model, 
                    coef = "disease", 
                    cell_type = "Bcell")
## Test the effect of disease in Bcell.
head(res_table, 3)
##                  beta   beta_var          mu effect_size f_statistics
## cg07075387 -0.7321689 0.03946395 -0.06071048    1.715505     13.58382
## cg13293535 -0.4776350 0.01825102  0.03636267    2.359218     12.49986
## cg15300101 -0.4525076 0.01625961 -0.13260057    1.260978     11.37536
##                 p_value       fdr
## cg07075387 0.0004255055 0.9330659
## cg13293535 0.0006969280 0.9330659
## cg15300101 0.0011733899 0.9330659

Instead of using the names of single coefficient, you can specify contrast levels, i.e. the comparing levels in this coefficient. For example, using male (gender = 1) as reference, testing female (gender = 2) effect in CD4T:

res_table <- csTest(fitted_model, 
                    coef = c("gender", 2, 1), 
                    cell_type = "CD4T")
## Test the effect of gender level 2 vs. level 1 in CD4T.
head(res_table, 3)
##            f_statistics      p_value       fdr
## cg17959722     15.41598 0.0001881780 0.4919237
## cg15473904     14.16328 0.0003279492 0.4919237
## cg09651654     12.39320 0.0007319324 0.7319324

5.2.2 Testing the joint effect of single parameter in all cell types.

For example, testing the joint effect of age in all cell types:

res_table <- csTest(fitted_model, 
                    coef = "age", 
                    cell_type = "joint")
head(res_table, 3)

Specifying cell_type as NULL or not specifying cell_type will test the effect in each cell type and the joint effect in all cell types.

res_table <- csTest(fitted_model, 
                    coef = "age", 
                    cell_type = NULL)
lapply(res_table, head, 3)

## this is exactly the same as
res_table <- csTest(fitted_model, 
                    coef = "age")

5.3 Detect cross-cell type differential signals

  1. Assuming you have TOAST library and dataset loaded, first step is to generate the study design based on the phenotype matrix. We allow general design matrix such as the following:
design <- data.frame(age = Pheno$age,
                     gender = as.factor(Pheno$gender),
                     disease = as.factor(Pheno$disease))

Prop <- estProp_RF_improved
colnames(Prop) <- colnames(Ref)  ## columns of proportion matrix should have names

Note that if all subjects belong to one group, we also allow detecting cross-cell type differences. In this case, the design matrix can be specified as:

design <- data.frame(disease = as.factor(rep(0,100)))
  1. Make model design using the design (phenotype) data frame and proportion matrix.
Design_out <- makeDesign(design, Prop)
  1. Fit linear models for raw data and the design generated from Design_out().
fitted_model <- fitModel(Design_out, Y_raw)
# print all the cell type names
fitted_model$all_cell_types
## [1] "CD8T"  "CD4T"  "NK"    "Bcell" "Mono"  "Gran"
# print all phenotypes
fitted_model$all_coefs
## [1] "age"     "gender"  "disease"

For cross-cell type differential signal detection, TOAST also allows multiple ways for testing. For example

5.3.1 Testing cross-cell type differential signals in cases (or in controls).

For example, testing the differences between CD8T and B cells in case group

test <- csTest(fitted_model, 
               coef = c("disease", 1), 
               cell_type = c("CD8T", "Bcell"), 
               contrast_matrix = NULL)
## Test the differences of CD8T vs. Bcell in disease:1.
head(test, 3)
##            f_statistics      p_value       fdr
## cg22692549     12.81337 0.0006037238 0.5101987
## cg10543947     12.64785 0.0006512032 0.5101987
## cg00079898     12.62378 0.0006584236 0.5101987

Or testing the differences between CD8T and B cells in control group

test <- csTest(fitted_model, 
               coef = c("disease", 0), 
               cell_type = c("CD8T", "Bcell"), 
               contrast_matrix = NULL)
## Test the differences of CD8T vs. Bcell in disease:0.
head(test, 3)
##            f_statistics     p_value       fdr
## cg00079898     14.85597 0.000240920 0.6052438
## cg06888460     11.64838 0.001033055 0.6052438
## cg16661522     11.52106 0.001096201 0.6052438

5.3.2 Testing the overall cross-cell type differences in all samples.

For example, testing the overall differences between Gran and CD4T in all samples, regardless of phenotypes.

test <- csTest(fitted_model, 
               coef = "joint", 
               cell_type = c("Gran", "CD4T"), 
               contrast_matrix = NULL)
## Test the joint effect of Gran vs. CD4T.
head(test, 3)
##            f_statistics      p_value       fdr
## cg05364038     17.55640 7.448377e-05 0.1770472
## cg14036627     16.48538 1.180315e-04 0.1770472
## cg10785373     15.37910 1.912542e-04 0.1912542

If you do not specify coef but only the two cell types to be compared, TOAST will test the differences of these two cell types in each coef parameter and the overall effect.

test <- csTest(fitted_model, 
               coef = NULL, 
               cell_type = c("Gran", "CD4T"), 
               contrast_matrix = NULL)
## Test the difference of Gran vs. CD4T in different values of age.
## Test the difference of Gran vs. CD4T in different values of gender.
## Test the difference of Gran vs. CD4T in different values of disease.
## Test the joint effect of Gran vs. CD4T.
lapply(test, head, 3)
## $age
##            f_statistics      p_value       fdr
## cg14036627     16.73995 0.0001057335 0.1820143
## cg05364038     15.77666 0.0001606611 0.1820143
## cg10785373     15.49181 0.0001820143 0.1820143
## 
## $gender
##            f_statistics     p_value       fdr
## cg23622369     10.71461 0.001601029 0.9997172
## cg19801747     10.66346 0.001640256 0.9997172
## cg19962075     10.32512 0.001926138 0.9997172
## 
## $disease
##            f_statistics      p_value        fdr
## cg04021544     22.20766 1.083840e-05 0.03251519
## cg12036633     16.95815 9.624687e-05 0.14437031
## cg09982942     12.65167 6.500670e-04 0.41857616
## 
## $joint
##            f_statistics      p_value       fdr
## cg05364038     17.55640 7.448377e-05 0.1770472
## cg14036627     16.48538 1.180315e-04 0.1770472
## cg10785373     15.37910 1.912542e-04 0.1912542

5.3.3 Testing the differences of two cell types over different values of one phenotype (higher-order test).

For example, testing the differences between Gran and CD4T in disease patients versus in controls.

test <- csTest(fitted_model, 
               coef = "disease", 
               cell_type = c("Gran", "CD4T"), 
               contrast_matrix = NULL)
## Test the difference of Gran vs. CD4T in different values of disease.
head(test, 3)
##            f_statistics      p_value        fdr
## cg04021544     22.20766 1.083840e-05 0.03251519
## cg12036633     16.95815 9.624687e-05 0.14437031
## cg09982942     12.65167 6.500670e-04 0.41857616

For another example, testing the differences between Gran and CD4T in males versus females.

test <- csTest(fitted_model, 
               coef = "gender", 
               cell_type = c("Gran", "CD4T"), 
               contrast_matrix = NULL)
## Test the difference of Gran vs. CD4T in different values of gender.
head(test, 3)
##            f_statistics     p_value       fdr
## cg23622369     10.71461 0.001601029 0.9997172
## cg19801747     10.66346 0.001640256 0.9997172
## cg19962075     10.32512 0.001926138 0.9997172

5.4 A few words about variance bound and Type I error.

5.4.1 Variance bound

There is an argument in csTest() called var_shrinkage. var_shrinkage is whether to apply shrinkage on estimated mean squared errors (MSEs) from the regression. Based on our experience, extremely small variance estimates sometimes cause unstable test statistics. In our implementation, use the 10% quantile value to bound the smallest MSEs. We recommend to use the default opinion var_shrinkage = TRUE.

5.4.2 Type I error

For all the above tests, we implement them using F-test. In our own experiments, we observe inflated type I errors from using F-test. As a result, we recommend to perform a permutation test to validate the significant signals identified are “real”.

Session info

## R version 3.6.1 (2019-07-05)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 18.04.3 LTS
## 
## Matrix products: default
## BLAS:   /home/biocbuild/bbs-3.10-bioc/R/lib/libRblas.so
## LAPACK: /home/biocbuild/bbs-3.10-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] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] TOAST_1.0.0      EpiDISH_2.2.0    RefFreeEWAS_2.2  quadprog_1.5-7  
## [5] BiocStyle_2.14.0
## 
## loaded via a namespace (and not attached):
##  [1] Rcpp_1.0.2                  compiler_3.6.1             
##  [3] BiocManager_1.30.9          GenomeInfoDb_1.22.0        
##  [5] XVector_0.26.0              bitops_1.0-6               
##  [7] class_7.3-15                tools_3.6.1                
##  [9] zlibbioc_1.32.0             digest_0.6.22              
## [11] lattice_0.20-38             evaluate_0.14              
## [13] rlang_0.4.1                 Matrix_1.2-17              
## [15] DelayedArray_0.12.0         yaml_2.2.0                 
## [17] parallel_3.6.1              csSAM_1.2.4                
## [19] xfun_0.10                   GenomeInfoDbData_1.2.2     
## [21] e1071_1.7-2                 stringr_1.4.0              
## [23] knitr_1.25                  S4Vectors_0.24.0           
## [25] IRanges_2.20.0              stats4_3.6.1               
## [27] grid_3.6.1                  Biobase_2.46.0             
## [29] BiocParallel_1.20.0         rmarkdown_1.16             
## [31] bookdown_0.14               magrittr_1.5               
## [33] matrixStats_0.55.0          htmltools_0.4.0            
## [35] MASS_7.3-51.4               BiocGenerics_0.32.0        
## [37] GenomicRanges_1.38.0        SummarizedExperiment_1.16.0
## [39] stringi_1.4.3               RCurl_1.95-4.12