TOAST.R
In TOAST: Tools for the analysis of heterogeneous tissues

## ----install, eval = FALSE----------------------------------------------------
#  if (!requireNamespace("BiocManager", quietly = TRUE))
#      install.packages("BiocManager")
#  
#  BiocManager::install("TOAST")

## ----quick_start, eval = FALSE------------------------------------------------
#  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)

## ----loadData-----------------------------------------------------------------
library(TOAST)
data("RA_100samples")
Y_raw <- RA_100samples$Y_raw
Pheno <- RA_100samples$Pheno
Blood_ref <- RA_100samples$Blood_ref

## ----checkData----------------------------------------------------------------
dim(Y_raw) 
Y_raw[1:4,1:4]

## ----checkPheno---------------------------------------------------------------
dim(Pheno)
head(Pheno, 3)

## ----checkRef-----------------------------------------------------------------
dim(Blood_ref)
head(Blood_ref, 3)

## ----loadDataCBS--------------------------------------------------------------
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

## ----checkDataCBS-------------------------------------------------------------
dim(CBS_mix) 
CBS_mix[1:4,1:4]
head(CBS_trueProp, 3)

## ----checkLM5-----------------------------------------------------------------
head(LM_5, 3)

## ----checkPrior---------------------------------------------------------------
prior_alpha
prior_sigma

## ----loadDatabeta-------------------------------------------------------------
data("beta_emp")
Ybeta = beta_emp$Y.raw
ref_m = beta_emp$ref.m

## ----checkDatabeta------------------------------------------------------------
dim(Ybeta) 
Ybeta[1:4,1:4]

## ----checkDataRef-------------------------------------------------------------
head(ref_m, 3)

## ----SelFeature---------------------------------------------------------------
refinx <- findRefinx(Y_raw, nmarker = 1000)

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

## ----DB2----------------------------------------------------------------------
library(EpiDISH)
outT <- epidish(beta.m = Y, ref.m = Ref, method = "RPC")
estProp_RB <- outT$estF

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

## ----DF, message=FALSE, warning=FALSE-----------------------------------------
library(RefFreeEWAS)

## ----DF2----------------------------------------------------------------------
refinx <- findRefinx(Y_raw, nmarker = 1000)
Y <- Y_raw[refinx,]

## ---- DF3, results='hide', message=FALSE, warning=FALSE-----------------------
K <- 6
outT <- RefFreeCellMix(Y, mu0=RefFreeCellMixInitialize(Y, K = K))
estProp_RF <- outT$Omega

## ----compareRFRB--------------------------------------------------------------
# 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)))

## ----IRB-RFCM1, message=FALSE, warning=FALSE----------------------------------
library(TOAST)

## ----IRB-RFCM2, results='hide', message=FALSE, warning=FALSE------------------
K=6
set.seed(1234)
outRF1 <- csDeconv(Y_raw, K, TotalIter = 30, bound_negative = TRUE) 

## ----IRB-RFCM3, message=FALSE, warning=FALSE----------------------------------
## check the accuracy of deconvolution
estProp_RF_improved <- assignCellType(input=outRF1$estProp,
                                      reference=estProp_RB) 
mean(diag(cor(estProp_RF_improved, estProp_RB)))

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

## ---- eval = FALSE------------------------------------------------------------
#  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)

## ----chooseMarker-------------------------------------------------------------
## 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)

## ----PRFwithoutPrior----------------------------------------------------------
resCBS0 <- MDeconv(CBS_mix, myMarker,
                epsilon = 1e-3, verbose = FALSE)
diag(cor(CBS_trueProp, t(resCBS0$H)))
mean(abs(as.matrix(CBS_trueProp) - t(resCBS0$H)))

## ----manualalpha--------------------------------------------------------------
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) 

## ----getprior-----------------------------------------------------------------
thisprior <- GetPrior("human pbmc")
thisprior

## ----PRFwithPrior-------------------------------------------------------------
resCBS1 <- MDeconv(CBS_mix, myMarker,
                alpha = prior_alpha,
                sigma = prior_sigma,
                epsilon = 1e-3, verbose = FALSE)
diag(cor(CBS_trueProp, t(resCBS1$H)))
mean(abs(as.matrix(CBS_trueProp) - t(resCBS1$H)))

## ----PRFwithPrior2------------------------------------------------------------
resCBS2 <- MDeconv(CBS_mix, myMarker,
                   alpha = "human pbmc",
                   epsilon = 1e-3, verbose = FALSE)
diag(cor(CBS_trueProp, t(resCBS2$H)))
mean(abs(as.matrix(CBS_trueProp) - t(resCBS2$H)))

## ----expKRef,results='hide',message=FALSE, warning=FALSE----------------------
out = Tsisal(Ybeta,K = 4, knowRef = ref_m)
out$estProp[1:3,1:4]
head(out$selMarker)

## ----expAllNULL,results='hide',message=FALSE, warning=FALSE-------------------
out = Tsisal(Ybeta,K = NULL, knowRef = NULL, possibleCellNumber = 2:5)
out$estProp[1:3,1:out$K]
head(out$selMarker)
out$K

## ----expKNULL,results='hide',message=FALSE, warning=FALSE---------------------
out = Tsisal(Ybeta, K = NULL, knowRef = ref_m, possibleCellNumber = 2:5)
out$estProp[1:3,1:out$K]
head(out$selMarker)
out$K

## ----csDE2--------------------------------------------------------------------
head(Pheno, 3)
design <- data.frame(disease = as.factor(Pheno$disease))

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


## ----csDE3--------------------------------------------------------------------
Design_out <- makeDesign(design, Prop)

## ----csDE4--------------------------------------------------------------------
fitted_model <- fitModel(Design_out, Y_raw)
# print all the cell type names
fitted_model$all_cell_types
# print all phenotypes
fitted_model$all_coefs

## -----------------------------------------------------------------------------
res_table <- csTest(fitted_model, 
                    coef = "disease", 
                    cell_type = "Gran")
head(res_table, 3)
Disease_Gran_res <- res_table

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

## ----joint, eval = FALSE------------------------------------------------------
#  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")

## ----general2-----------------------------------------------------------------
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

## ----general3-----------------------------------------------------------------
Design_out <- makeDesign(design, Prop)

## ----general4-----------------------------------------------------------------
fitted_model <- fitModel(Design_out, Y_raw)
# print all the cell type names
fitted_model$all_cell_types
# print all phenotypes
fitted_model$all_coefs

## ----general5-----------------------------------------------------------------
res_table <- csTest(fitted_model, 
                    coef = "age", 
                    cell_type = "Gran")
head(res_table, 3)

## ----general6-----------------------------------------------------------------
res_table <- csTest(fitted_model, 
                    coef = "disease", 
                    cell_type = "Bcell")
head(res_table, 3)

## -----------------------------------------------------------------------------
res_table <- csTest(fitted_model, 
                    coef = c("gender", 2, 1), 
                    cell_type = "CD4T")
head(res_table, 3)

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

## ---- eval = FALSE------------------------------------------------------------
#  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")

## ----crossCellType2-----------------------------------------------------------
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

## ----crossCellType3, eval = FALSE---------------------------------------------
#  design <- data.frame(disease = as.factor(rep(0,100)))

## ----crossCellType4-----------------------------------------------------------
Design_out <- makeDesign(design, Prop)

## ----crossCellType5-----------------------------------------------------------
fitted_model <- fitModel(Design_out, Y_raw)
# print all the cell type names
fitted_model$all_cell_types
# print all phenotypes
fitted_model$all_coefs

## -----------------------------------------------------------------------------
test <- csTest(fitted_model, 
               coef = c("disease", 1), 
               cell_type = c("CD8T", "Bcell"), 
               contrast_matrix = NULL)
head(test, 3)

## -----------------------------------------------------------------------------
test <- csTest(fitted_model, 
               coef = c("disease", 0), 
               cell_type = c("CD8T", "Bcell"), 
               contrast_matrix = NULL)
head(test, 3)

## -----------------------------------------------------------------------------
test <- csTest(fitted_model, 
               coef = "joint", 
               cell_type = c("Gran", "CD4T"), 
               contrast_matrix = NULL)
head(test, 3)

## -----------------------------------------------------------------------------
test <- csTest(fitted_model, 
               coef = NULL, 
               cell_type = c("Gran", "CD4T"), 
               contrast_matrix = NULL)
lapply(test, head, 3)

## -----------------------------------------------------------------------------
test <- csTest(fitted_model, 
               coef = "disease", 
               cell_type = c("Gran", "CD4T"), 
               contrast_matrix = NULL)
head(test, 3)

## -----------------------------------------------------------------------------
test <- csTest(fitted_model, 
               coef = "gender", 
               cell_type = c("Gran", "CD4T"), 
               contrast_matrix = NULL)
head(test, 3)

## ----sessionInfo, echo=FALSE--------------------------------------------------
sessionInfo()