% Generated by roxygen2: do not edit by hand
% Please edit documentation in R/Simulate.R
\name{simulateDE}
\alias{simulateDE}
\title{Simulate Differential Expression Pipeline}
\usage{
simulateDE(SetupRes,
Prefilter = NULL,
Imputation = NULL,
Normalisation = c("TMM", "MR", "PosCounts", "UQ",
"scran", "Linnorm", "sctransform",
"SCnorm", "Census", "depth"),
Label = "none",
DEmethod = c("T-Test", "edgeR-LRT", "edgeR-QL",
"edgeR-zingeR", "edgeR-ZINB-WaVE",
"limma-voom", "limma-trend",
"DESeq2", "DESeq2-zingeR", "DESeq2-ZINB-WaVE",
"ROTS", "baySeq", "NOISeq", "EBSeq",
"MAST", "BPSC", "scDD", "DECENT"),
DEFilter = FALSE,
Counts = FALSE,
NCores = NULL,
verbose = TRUE)
}
\arguments{
\item{SetupRes}{This object specifies the simulation setup.
This must be the return object from \code{\link{Setup}}.}

\item{Prefilter}{A character vector specifying the gene expression filtering method
to be used prior to normalisation (and possibly imputation).
Default is \code{NULL}, i.e. no filtering.
Please consult the Details section for available options.}

\item{Imputation}{A character vector specifying the gene expression imputation method
to be used prior to normalisation.
Default is \code{NULL}, i.e. no imputation.
Please consult the Details section for available options.}

\item{Normalisation}{Normalisation method to use.
Please consult the Details section for available options.}

\item{Label}{A character vector to define whether information about
group labels should be used for normalisation.
This is only implemented for scran and SCnorm.
Possible options include the default \code{"none"}
which means that no sample group information is considered for normalisation;
\code{"known"} means that the simulated group labels are used and \code{"clustering"}
which applies an unsupervised hierarchical clustering to determine the group labels
(see \code{\link[scran]{quickCluster}} for details).}

\item{DEmethod}{A character vector specifying the DE detection method to be used.
Please consult the Details section for available options.}

\item{DEFilter}{A logical vector indicating whether to run DE testing on
filtered and/or imputed count data.
Default is \code{FALSE}.}

\item{Counts}{A logical vector indicating whether the simulated count matrix is also provided as output.
Default is \code{FALSE} since the output can be quite large. Note that if DEFilter is \code{TRUE},
then the returned count matrix will countain the filtered and/or imputed count data.}

\item{NCores}{integer positive number of cores for parallel processing.
Default is \code{NULL}, i.e. 1 core.}

\item{verbose}{Logical vector to indicate whether to show progress report of simulations.
Default is \code{TRUE}.}
}
\value{
\strong{SimulateRes: Results of DE simulations}
\describe{
\item{pvalue, fdr}{3D array (ngenes * N * nsims) for p-values and FDR from each simulation.
Note that FDR values will be empty and the calculation will be done by \code{\link{evaluateDE}} whenever applicable.}
\item{mu,disp,dropout}{3D (ngenes * N * nsims) array for mean, dispersion and dropout of library size factor normalized read counts.}
\item{true.mu,true.disp,true.dropout}{3D (ngenes * N * nsims) array for true mean, dispersion and dropout of simulated read counts.}
\item{true.depth,est.depth}{True simulated and processed (after prefiltering and/or imputation if defined) sequencing depth per sample.}
\item{est.sf}{Global library size factor estimates per sample.}
\item{est.gsf}{3D array (ngenes * N * nsims) for size factor estimates. These are gene- and sample-wise estimates and only for SCnorm and Linnorm normalisation.}
\item{elfc,rlfc}{3D array (ngenes * N * nsims) for log fold changes (LFC):
elfc is for the DE tool estimated LFC; rlfc is for the LFC estimated from the normalised read counts.}
\item{time.taken}{The time taken given by \code{\link[base]{proc.time}} for each simulation, given for preprocessing, normalisation, differential expression testing and moment estimation.}
}
\strong{SetupRes: Simulation specifications}
\describe{
\item{DESetup - ... - estSpikeRes}{Reiterating the simulated setup defined by \code{\link{Setup}}.}
\item{Pipeline}{A list of chosen pipeline tools defined by above arguments.}
}
\strong{Counts: Simulated Count Matrices}
\describe{
\item{Counts}{List of lists object where \code{Counts[[Sample Size Setup]][[Simulation Run]]} containing simulated counts. Note that this will only be returned when \code{Counts} is \code{TRUE}. In addition, if \code{DEFilter} is \code{TRUE} then the filtered/imputed counts are returned.}
}
}
\description{
simulateDE is the main function to simulate differential expression for RNA-seq experiments.
The simulation parameters are specified with \code{\link{Setup}}.
The user now needs to specify the RNA-seq Analysis Pipeline including preprocessing, normalisation and differential testing method.
The return object contains DE test results from all simulations as well as descriptive statistics.
The error matrix calculations will be conducted with \code{\link{evaluateDE}}.\cr
}
\details{
Here you can find detailed information about preprocessing, imputation, normalisation and differential testing choices.
For recommendations concerning single cell RNA-sequencing pipelines,
we kindly refer the user to \href{https://www.nature.com/articles/s41467-019-12266-7}{Vieth, et al (2019). A systematic evaluation of single cell RNA-seq analysis pipelines. Nature Communications, 10(1), 4667}.\cr
\strong{Prefiltering}\cr
\describe{
\item{CountFilter}{removes genes that have a mean expression below 0.2.}
\item{FreqFilter}{removes genes that have more than 80 \% dropouts.}
}
\strong{Imputation}\cr
\describe{
\item{scImpute}{apply the imputation as implemented in \code{\link[scImpute]{scimpute}}}
\item{DrImpute}{apply the imputation as implemented in \code{\link[DrImpute]{DrImpute}}.}
\item{SAVER}{apply the imputation as implemented in \code{\link[SAVER]{saver}}.}
\item{scone}{apply the imputation as implemented in \code{\link[scone]{scone}}, defining 'house keeping genes' for the FNR model estimation as those that have less than 20 \% dropouts and small variance (i.e. in the lower 20th quartile). If less than 25 genes could be identified, the genes with less than 5 \% dropouts are used. If spike-in data is provided then these are used for the FNR model estimation.}
\item{MAGIC}{apply the imputation as implemented in \code{\link[Rmagic]{magic}}. Please note that the python package MAGIC needs to be installed to use this implementation.}
}
\strong{Normalisation}\cr
\describe{
\item{TMM, UQ}{employ the edgeR style normalization of weighted trimmed mean of M-values and upperquartile
as implemented in \code{\link[edgeR]{calcNormFactors}}, respectively.}
\item{MR, PosCounts}{employ the DESeq2 style normalization of median ratio method and a modified geometric mean method
as implemented in \code{\link[DESeq2]{estimateSizeFactors}}, respectively.}
\item{scran, SCnorm}{apply the deconvolution and quantile regression normalization methods developed for sparse RNA-seq data
as implemented in \code{\link[scran]{computeSumFactors}} and \code{\link[SCnorm]{SCnorm}}, respectively. Spike-ins can also be supplied for both methods via \code{spikeData}. Note, however that this means for scran that the normalisation as implemented in \code{\link[scran]{calculateSumFactors}} is also applied to genes (\code{general.use=TRUE}).}
\item{Linnorm}{apply the normalization method for sparse RNA-seq data
as implemented in \code{\link[Linnorm]{Linnorm.Norm}}.
For \code{Linnorm}, the user can also supply \code{spikeData}.}
\item{sctransform}{apply the normalization method developed for single-cell
UMI RNA-seq data as implemented in \code{\link[sctransform]{vst}}. }
\item{Census}{converts relative measures of TPM/FPKM values into mRNAs per cell (RPC) without the need of spike-in standards.
Census at least needs \code{Lengths} for single-end data and preferably \code{MeanFragLengths} for paired-end data.
The authors state that Census should not be used for UMI data.}
\item{depth}{Sequencing depth normalisation.}
}
\strong{Differential testing}\cr
\describe{
\item{T-Test}{A T-Test per gene is applied using log transformed and normalized expression values (i.e. CPM or TPM).}
\item{limma-trend, limma-voom}{apply differential testing as implemented in \code{\link[limma]{lmFit}}
followed by \code{\link[limma]{eBayes}} on counts transformed by \code{\link[limma]{voom}} or by applying mean-variance trend on log CPM values in \code{\link[limma]{eBayes}}.}
\item{edgeR-LRT, edgeR-QL}{apply differential testing as implemented in \code{\link[edgeR]{glmFit}}, \code{\link[edgeR]{glmLRT}} and\code{\link[edgeR]{glmQLFit}}, \code{\link[edgeR]{glmQLFTest}}, respectively.}
\item{DESeq2}{applies differential testing as implemented in \code{\link[DESeq2]{DESeq}}.}
\item{ROTS}{applies differential testing as implemented in \code{\link[ROTS]{ROTS}} with 100 permutations on transformed counts (\code{\link[limma]{voom}}).}
\item{baySeq}{applies differential testing as implemented in \code{\link[baySeq]{getLikelihoods}} based on negative binomial prior estimates (\code{\link[baySeq]{getPriors.NB}}).}
\item{NOISeq}{applies differential testing as implemented in \code{\link[NOISeq]{noiseqbio}} based on CPM values.}
\item{EBSeq}{applies differential testing as implemented in \code{\link[EBSeq]{EBTest}}.}
\item{MAST}{applies differential testing as implemented in \code{\link[MAST]{zlm}} for zero-inflated model fitting followed by \code{\link[MAST]{lrTest}} on log CPM values.}
\item{BPSC}{applies differential testing as implemented in \code{\link[BPSC]{BPglm}} on CPM values.}
\item{scDD}{applies differential testing as implemented in \code{\link[scDD]{scDD}} on CPM values.}
\item{DECENT}{applies differential testing as implemented in \code{\link[DECENT]{decent}}.}
\item{edgeR-zingeR, DESeq2-zingeR}{In a first step, the posterior probabilities of the zero-inflated negative binomial component are estimated (see \code{\link[zingeR]{zeroWeightsLS}}) and used to define a weight matrix for dispersion estimation in \code{\link[edgeR]{estimateDisp}}. For the edgeR approach, the generalized model as implemented in \code{\link[edgeR]{glmFit}} is fitted. This is followed by an adapted LRT for differential testing to account for the weighting (see \code{\link[zingeR]{glmWeightedF}}). For DESeq2, the generalized linear model coefficients are estimated using \code{\link[DESeq2]{nbinomWaldTest}} and the weighting is done by setting the degrees of freedom for the T distribution.}
\item{edgeR-ZINB-WaVE, DESeq2-ZINB-WaVE}{In a first step, a zero-inflated negative binomial regression model  is fitted (see \code{\link[zinbwave]{zinbFit}}) to estimate observational weights (see \code{\link[zinbwave]{computeObservationalWeights}}) used for dispersion estimation in \code{\link[edgeR]{estimateDisp}}. For the edgeR approach, the generalized model as implemented in \code{\link[edgeR]{glmFit}} is fitted. This is followed by an adapted LRT for differential testing to account for the weighting (see \code{\link[zinbwave]{glmWeightedF}}). For DESeq2, the generalized linear model coefficients are estimated using \code{\link[DESeq2]{nbinomWaldTest}} and the weighting is done by setting the degrees of freedom for the T distribution.}
}
}
\examples{
\dontrun{
# estimate gene parameters
data("SmartSeq2_Gene_Read_Counts")
Batches = data.frame(Batch = sapply(strsplit(colnames(SmartSeq2_Gene_Read_Counts),
                                    "_"), "[[", 1),
                     stringsAsFactors = F,
                     row.names = colnames(SmartSeq2_Gene_Read_Counts))
data("GeneLengths_mm10")
estparam_gene <- estimateParam(countData = SmartSeq2_Gene_Read_Counts,
                               readData = NULL,
                              batchData = Batches,
                               spikeData = NULL, spikeInfo = NULL,
                               Lengths = GeneLengths_mm10, MeanFragLengths = NULL,
                               RNAseq = 'singlecell', Protocol = 'Read',
                               Distribution = 'ZINB', Normalisation = "scran",
                               GeneFilter = 0.1, SampleFilter = 3,
                               sigma = 1.96, NCores = NULL, verbose = TRUE)
# estimate spike parameters
data("SmartSeq2_SpikeIns_Read_Counts")
data("SmartSeq2_SpikeInfo")
Batches = data.frame(Batch = sapply(strsplit(colnames(SmartSeq2_SpikeIns_Read_Counts),
                                    "_"), "[[", 1),
                     stringsAsFactors = F,
                     row.names = colnames(SmartSeq2_SpikeIns_Read_Counts))
estparam_spike <- estimateSpike(spikeData = SmartSeq2_SpikeIns_Read_Counts,
                                spikeInfo = SmartSeq2_SpikeInfo,
                                MeanFragLength = NULL,
                                batchData = Batches,
                                Normalisation = 'depth')
# define log fold change
p.lfc <- function(x) sample(c(-1,1), size=x,replace=T)*rgamma(x, shape = 1, rate = 2)
# set up simulations
setupres <- Setup(ngenes = 10000, nsims = 10,
                  p.DE = 0.1, pLFC = p.lfc,
                  n1 = c(20,50,100), n2 = c(30,60,120),
                  Thinning = c(1,0.9,0.8), LibSize = 'given',
                  estParamRes = estparam_gene,
                  estSpikeRes = estparam_spike,
                  DropGenes = FALSE,
                  sim.seed = 52679, verbose = TRUE)
# run simulation
simres <- simulateDE(SetupRes = setupres,
                     Prefilter = "FreqFilter", Imputation = NULL,
                     Normalisation = 'scran', Label = 'none',
                     DEmethod = "limma-trend", DEFilter = FALSE,
                     NCores = NULL, verbose = TRUE)
# quick evaluation
evalderes <- evaluateDE(simRes = simres)
# plot evaluation
plotEvalDE(evalRes = evalderes, rate = "marginal",
           quick = TRUE, Annot = FALSE)
}
}
\seealso{
\code{\link{estimateParam}} and \code{\link{estimateSpike}},  for parameter specifications;\cr
 \code{\link{Setup}} for simulation setup;\cr
 \code{\link{evaluateDE}}, \code{\link{evaluateROC}} and \code{\link{evaluateSim}} for evaluation of simulation results.
}
\author{
Beate Vieth
}
