R/Evaluate.R
evaluateSim.RdThis function takes the simulation output from simulateDE
and computes several metrics that give an indication of the simulation setup performance.
evaluateSim(simRes)
| simRes | The result from |
|---|
A list with the following entries:
The absolute mean error (MAE), root mean square error (RMSE)
and root mean square residual error of a robust linear model (rRMSE, rlm)
of log fold change differences between estimated and simulated.
Furthermore, the fraction of missing entries (NAFraction) for MAE annd RMSE.
The median absolute deviation (MAD) between the simulated and estimated size factors,
the root mean square residual error of a robust linear model (rRMSE, rlm) and
the group-specific ratio between simulated and estimated size factors (GroupX).
estimateParam for negative binomial parameters,
Setup for setting up simulation parameters and
simulateDE for simulating differential expression
if (FALSE) { # 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) # 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 = NULL, DropGenes = FALSE, sim.seed = 66437, 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) # evaluation evalsimres <- evaluateSim(simRes = simres) plotEvalSim(evalRes = evalsimres, Annot = TRUE) plotTime(evalRes = evalsimres, Annot = TRUE) }