TOAST 1.4.0

- 1 Introduction
- 2 Installation and quick start
- 3 Example dataset
- 4 Estimate mixing proportions
- 4.1 Reference-based deconvolution using least square method
- 4.2 Reference-free deconvolution using RefFreeEWAS
- 4.3 Improve reference-free deconvolution with cross-cell type differential analysis
- 4.4 Partial reference-free deconvolution (TOAST/-P and TOAST/+P)
- 4.5 Complete deconvolution using a geometric approach

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

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 averages 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:

- 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.

- 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. If cell type-specific markers (or prior knowledge of cell compositions) are available, TOAST provides partial reference-free deconvolution function, which is more accuracte than RF methods and works well even for very small sample size (e.g.<10).

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

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

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

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)
```

TOAST provides two sample dataset.

The first example dataset is 450K DNA methylation data. 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`

`## Loading required package: limma`

`## Loading required package: nnls`

```
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
```

The second example dataset is microarray gene expression data. We obtain and process this dataset based on the raw data provided by GSE65133. This microarary data is from 20 PBMC samples. The original dataset has 47323 probes. We mapped the probes into 21626 genes and then further reduced the dataset to 511 genes by selecting the genes that have matches in reference panel.

```
data("CBS_PBMC_array")
CBS_mix <- CBS_PBMC_array$mixed_all
LM_5 <- CBS_PBMC_array$LM_5
CBS_trueProp <- CBS_PBMC_array$trueProp
prior_alpha <- CBS_PBMC_array$prior_alpha
prior_sigma <- CBS_PBMC_array$prior_sigma
```

Check the PBMC microarray gene expression data and true proportions

`dim(CBS_mix) `

`## [1] 511 20`

`CBS_mix[1:4,1:4]`

```
## X17.002 X17.006 X17.019 X17.023
## ABCB4 96.0 107.50 110.00 92.3
## ABCB9 98.3 109.75 103.85 92.1
## ACAP1 196.8 217.80 351.00 140.7
## ACHE 92.7 97.20 87.10 87.1
```

`head(CBS_trueProp, 3)`

```
## Bcells CD8T CD4T NKcells Monos
## 17-002 0.16201354 0.2364636 0.2453469 0.17533841 0.18083756
## 17-006 0.05279853 0.3279443 0.4698660 0.09106723 0.05832395
## 17-019 0.21897143 0.2041143 0.1468571 0.33874286 0.09131429
```

Check reference matrix for 5 immune cell types

`head(LM_5, 3)`

```
## BCells CD8T CD4T NK cells Monocytes
## ABCB4 283.22884 4.31128 6.685426 9.119776 6.202496
## ABCB9 18.84917 24.22372 34.725041 19.129933 20.309426
## ACAP1 268.46349 1055.61338 1017.711114 450.109326 190.879024
```

Check prior knowledge for the 5 cell types

`prior_alpha`

`## [1] 0.09475324 0.23471057 0.33231687 0.09689958 0.24131974`

`prior_sigma`

`## [1] 0.0996325 0.1441782 0.1602440 0.1006351 0.1455614`

The third example dataset is a list containing two matrices, one of which is methylation 450K array data of 3000 CpG sites on 50 samples, the other is methylation 450K array data of 3000 matched CpG sites on three immune cell types. The first dataset is generated by simulation. It originally has 459226 features and 50 samples.We reduce it to 3000 CpGs by random selection.

```
data("beta_emp")
Ybeta = beta_emp$Y.raw
ref_m = beta_emp$ref.m
```

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

`dim(Ybeta) `

`## [1] 3000 50`

`Ybeta[1:4,1:4]`

```
## [,1] [,2] [,3] [,4]
## cg08752431 0.76611838 0.76117075 0.79835014 0.79239766
## cg14555682 0.09677872 0.08096398 0.11906131 0.10553938
## cg23086843 0.88824253 0.88740138 0.92060008 0.92240033
## cg20308511 0.02863965 0.03976921 0.04269377 0.02678175
```

Check reference matrix for 3000 CpGs by three immune cell types

`head(ref_m, 3)`

```
## CD4T CD8T BCell
## cg08752431 0.77118808 0.7546620 0.7589832
## cg14555682 0.09960192 0.1201291 0.1006569
## cg23086843 0.90647715 0.8810441 0.9065662
```

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, which works for both gene expression data and DNA methylation data. The example in Section 4.3 demonstrates the usage of this. Note that we have only 3000 features in the Y_raw from RA_100samples 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.

In Secion 4.4, we demonstrate the usage of partial reference-free (PRF) deconvolution. Compared to RB methods, PRF does not require reference panel thus can be more wdiely applied. Compared to RF methods, PRF uses additional biological information, which improves the estimation accuracy and automatically assign cell type labels.

- 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)`

- Subset data and reference panel.

```
Y <- Y_raw[refinx,]
Ref <- as.matrix(Blood_ref[refinx,])
```

- 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")`

- Load
*RefFreeEWAS*.

`library(RefFreeEWAS)`

- 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,]
```

- Do reference-free deconvolution on the RA dataset.

```
K <- 6
outT <- RefFreeCellMix(Y, mu0=RefFreeCellMixInitialize(Y, K = K))
estProp_RF <- outT$Omega
```

- 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`

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.

- Load TOAST package.

`library(TOAST)`

- 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)
```

- 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.

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)
```

Similar to DSA, our PRF method requires
the knowledge of **cell type-specific markers**.
Such markers can be selected from pure
cell type gene expression
profiles from same or different platforms (through
function `ChooseMarker()`

). They can also
be manually specified
(see function manual `?MDeconv`

for more explanation).
The **prior knowledge of cell compositions**
are optional, but
highly recommended. We find prior
knowledge of cell compositions (`alpha`

and `sigma`

)
help calibrate the scales of the estimations,
and reduce estimation bias. Such information
can be estimated from previous cell sorting
experiments or single cell study.
We currently provide prior knowledge for five
tissue types: “human pbmc”,“human
liver”, “human brain”, “human pancreas”, “human skin”,
which can be directly specified in `MDeconv()`

function.

We provide functions to choose cell type-specific markers from pure cell type profiles or single cell RNA-seq data. Here we demonstrate how to select markers from PBMC pure cell type gene expression profile.

```
## create cell type list:
CellType <- list(Bcells = 1,
CD8T = 2,
CD4T = 3,
NK = 4,
Monocytes = 5)
## choose (up to 20) significant markers
## per cell type
myMarker <- ChooseMarker(LM_5,
CellType,
nMarkCT = 20,
chooseSig = TRUE,
verbose = FALSE)
lapply(myMarker, head, 3)
```

```
## $Bcells
## [1] "BANK1" "MS4A1" "IGLL3P"
##
## $CD8T
## [1] "CD8B" "CD8A" "GZMK"
##
## $CD4T
## [1] "IL9" "CTLA4" "IL3"
##
## $NK
## [1] "KIR3DL2" "IL18RAP" "KLRF1"
##
## $Monocytes
## [1] "FCN1" "P2RY13" "NCF2"
```

```
resCBS0 <- MDeconv(CBS_mix, myMarker,
epsilon = 1e-3, verbose = FALSE)
```

`## Deconvolution without prior information.`

`diag(cor(CBS_trueProp, t(resCBS0$H)))`

`## [1] 0.5333925 0.5647335 0.7027891 0.6484607 0.7116513`

`mean(abs(as.matrix(CBS_trueProp) - t(resCBS0$H)))`

`## [1] 0.1313198`

We allow manually input the prior knowledge of all cell types, or select from currently supported tissues (“human pbmc”,“human liver”, “human brain”, “human pancreas”, “human skin”). Note that order of cell types in prior knowledge here should match the order in marker list.

Here is an example of manually specifying alpha and sigma:

```
prior_alpha <- c(0.09475, 0.23471, 0.33232, 0.0969, 0.24132)
prior_sigma <- c(0.09963, 0.14418, 0.16024, 0.10064, 0.14556)
names(prior_alpha) <- c("B cells", "CD8T", "CD4T",
"NK cells", "Monocytes")
names(prior_sigma) <- names(prior_alpha)
```

Here is to see alpha and sigma for supported tisuses using
`GetPrior()`

:

```
thisprior <- GetPrior("human pbmc")
thisprior
```

```
## $alpha_prior
## B cells CD8T CD4T NK cells Monocytes
## 0.09475 0.23471 0.33232 0.09690 0.24132
##
## $sigma_prior
## B cells CD8T CD4T NK cells Monocytes
## 0.09963 0.14418 0.16024 0.10064 0.14556
```

Deconvolution using manually input alpha and sigma:

```
resCBS1 <- MDeconv(CBS_mix, myMarker,
alpha = prior_alpha,
sigma = prior_sigma,
epsilon = 1e-3, verbose = FALSE)
```

`## Deconvolution with prior infromation.`

`diag(cor(CBS_trueProp, t(resCBS1$H)))`

`## [1] 0.5545308 0.5627714 0.6302584 0.6841369 0.7101932`

`mean(abs(as.matrix(CBS_trueProp) - t(resCBS1$H)))`

`## [1] 0.079472`

For supported tissues, you can directly specify tissue type as alpha input:

```
resCBS2 <- MDeconv(CBS_mix, myMarker,
alpha = "human pbmc",
epsilon = 1e-3, verbose = FALSE)
```

`## Deconvolution with prior infromation.`

`diag(cor(CBS_trueProp, t(resCBS2$H)))`

`## [1] 0.5545308 0.5627714 0.6302584 0.6841369 0.7101932`

`mean(abs(as.matrix(CBS_trueProp) - t(resCBS2$H)))`

`## [1] 0.079472`

Tsisal is a complete deconvolution method which estimates cell compositions from DNA methylation data without prior knowledge of cell types and their proportions. Tsisal is a full pipeline to estimate number of cell types, cell compositions, find cell-type-specific CpG sites, and assign cell type labels when (full or part of) reference panel is available.

Here is an example of manually specifying K and reference panel:

```
out = Tsisal(Ybeta,K = 4, knowRef = ref_m)
out$estProp[1:3,1:4]
head(out$selMarker)
```

Here is an example where both K and reference panel are unknown:

```
out = Tsisal(Ybeta,K = NULL, knowRef = NULL, possibleCellNumber = 2:5)
out$estProp[1:3,1:out$K]
head(out$selMarker)
out$K
```

Here is an example where K is unknown and reference panel is known:

```
out = Tsisal(Ybeta, K = NULL, knowRef = ref_m, possibleCellNumber = 2:5)
out$estProp[1:3,1:out$K]
head(out$selMarker)
out$K
```

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.

- 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
```

- Make model design using the design (phenotype) data frame and proportion matrix.

`Design_out <- makeDesign(design, Prop)`

- 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.

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`

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")
```

- 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
```

- Make model design using the design (phenotype) data frame and proportion matrix.

`Design_out <- makeDesign(design, Prop)`

- 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.

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
```

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")
```

- 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)))`

- Make model design using the design (phenotype) data frame and proportion matrix.

`Design_out <- makeDesign(design, Prop)`

- 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

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
```

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
```

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
```

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`

.

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”.

```
## R version 4.0.3 (2020-10-10)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 18.04.5 LTS
##
## Matrix products: default
## BLAS: /home/biocbuild/bbs-3.12-bioc/R/lib/libRblas.so
## LAPACK: /home/biocbuild/bbs-3.12-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.4.0 nnls_1.4 limma_3.46.0 EpiDISH_2.6.0
## [5] RefFreeEWAS_2.2 quadprog_1.5-8 BiocStyle_2.18.0
##
## loaded via a namespace (and not attached):
## [1] compiler_4.0.3 BiocManager_1.30.10
## [3] GenomeInfoDb_1.26.0 XVector_0.30.0
## [5] MatrixGenerics_1.2.0 bitops_1.0-6
## [7] class_7.3-17 tools_4.0.3
## [9] zlibbioc_1.36.0 digest_0.6.27
## [11] lattice_0.20-41 evaluate_0.14
## [13] rlang_0.4.8 Matrix_1.2-18
## [15] DelayedArray_0.16.0 yaml_2.2.1
## [17] parallel_4.0.3 csSAM_1.2.4
## [19] xfun_0.18 GenomeInfoDbData_1.2.4
## [21] e1071_1.7-4 stringr_1.4.0
## [23] knitr_1.30 S4Vectors_0.28.0
## [25] IRanges_2.24.0 grid_4.0.3
## [27] stats4_4.0.3 Biobase_2.50.0
## [29] rmarkdown_2.5 bookdown_0.21
## [31] corpcor_1.6.9 locfdr_1.1-8
## [33] magrittr_1.5 matrixStats_0.57.0
## [35] htmltools_0.5.0 MASS_7.3-53
## [37] splines_4.0.3 BiocGenerics_0.36.0
## [39] GenomicRanges_1.42.0 SummarizedExperiment_1.20.0
## [41] stringi_1.5.3 RCurl_1.98-1.2
```