This function takes the simulation output from simulateDE and computes receiver operator characteristics (ROC) and area under the curve values (AUC) with the help of functions implemented in calculate_performance.

evaluateROC(simRes,
alpha.type=c("adjusted","raw"),
MTC=c('BY', 'BH', 'Storey', 'IHW',
'holm', 'hochberg', 'hommel', 'bonferroni'),
alpha.nominal = 0.1,
target.by=c("lfc", "effectsize"), delta=0,
raw = FALSE, verbose = TRUE)

Arguments

simRes

The result from simulateDE.

alpha.type

A string to represent the way to call DE genes. Available options are "adjusted" i.e. applying multiple testing correction and "raw" i.e. using p-values. Default is "adjusted".

MTC

Multiple testing correction method to use. Available options are 1) see p.adjust.methods, 2) Storey's qvalue see qvalue and 3) Independent Hypothesis Weighting considering mean expression as covariate (see ihw). Default is BY, i.e. Benjamini-Yekutieli FDR correction method.

alpha.nominal

The nomial value of significance. Default is 0.1.

target.by

A string to specify the method to define "biologically important" DE genes. Available options are (1) "lfc": interesting genes are defined by absolute log2 fold changes. (2) "effectsize": interesting genes are defined by absolute log2 fold changes divided by the square root of 1/(mean+dispersion).

delta

A threshold used for defining "biologically important" genes. Genes with absolute log2 fold changes (when target.by is "lfc") or effect sizes (when target.by is "effectsize") greater than this value are deemed DE in error rates calculations. If the default delta=0, then no threshold is applied.

raw

Logical vector. The default FALSE removes intermediate calculations from the output since they can be quite large.

verbose

Logical vector to indicate whether to show progress report of evaluation. Default is TRUE.

Value

A list with the following entries:

Performances

The output of calculate_performance of aspect "fdrtprcurve" calculating the proportions of TN, FN, TP and FP and related rates which uses prediction and performance internally.

TPRvsFDR

The output of calculate_performance of aspect "fdrtpr" calculating the proportions of TN, FN, TP and FP and associated TPR and observed FDR for each nominal level from 0.01 to 1 in steps of 0.01.

Scores

The mean +/- standard error of performance measures and area under the curve. These are calculated once for conservative (extension "conv") and once for liberal nominal (extension "lib") FDR levels. Measures include: accuracy (ACC), F-measure of precision and recall (F1score), Matthew's correlation coefficient (MCC), partial Area under the Curve of TPR vs FDR (TPRvsFDR_pAUC), area under the ROC-curve (TPRvsFPR_AUC) and area under the PR-curve (TPRvsPPV_AUC).

Raw

If raw is TRUE, then the intermediate calculations of the above three list entries is also included in the output.

Settings

Reiterating chosen evaluation parameters.

See also

estimateParam for parameter estimation needed for simulations, Setup for setting up simulation parameters and simulateDE for simulating differential expression.

Examples

if (FALSE) { # estimate gene parameters data("SmartSeq2_Gene_Read_Counts") estparam_gene <- estimateParam(countData = SmartSeq2_Gene_Read_Counts, readData = NULL, batchData = NULL, spikeData = NULL, spikeInfo = NULL, Lengths = NULL, MeanFragLengths = NULL, RNAseq = 'singlecell', Protocol = 'Read', Distribution = 'ZINB', Normalisation = "scran", GeneFilter = 0.1, SampleFilter = 3, sigma = 1.96, NCores = NULL, verbose = TRUE) # define log2 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 = TRUE, sim.seed = 83596, 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 evalrocres <- evaluateROC(simRes = simres, alpha.type = "adjusted", MTC = 'BH', alpha.nominal = 0.05, raw = FALSE) # plot evaluation plotEvalROC(evalRes = evalrocres) }