Allows for the computation of differential analyses using
proDA
(Ahlmann-Eltze & Anders, 2020), a dropout-aware model that treats the missing values as
left-censored observations instead of replacing them. Includes means, Fold Changes, and p-values.
# S3 method for class 'analyses.proDA'
diff(
DEprot.object,
contrast.list,
include.rep.model = FALSE,
replicate.column = NULL,
linear.FC.th = 2,
linear.FC.unresp.range = c(1/1.1, 1.1),
padj.th = 0.05,
padj.method = "BH",
moderate.location = TRUE,
moderate.variance = TRUE,
min.detected = 1,
n.subsample = NULL,
max.iter = 20,
max.attempts = 3,
missingness.object = NULL,
check.missingness = TRUE,
up.color = "indianred",
down.color = "steelblue",
unresponsive.color = "purple",
null.color = "gray",
which.data = "normalized",
seed = NULL,
verbose = FALSE,
overwrite.analyses = FALSE
)An object of class DEprot.
List of 3-elements vectors indicating (in order): metadata_column, variable_1, variable_2.
Logical value indicating whether the model should include the replicate information (added as a blocking fixed effect). Requires a 'replicate.column'. Default: FALSE.
String indicating the name of a column from the metadata table in which are stored the replicate IDs. This column is used only if include.rep.model = TRUE. Default: NULL.
Number indicating the (absolute) fold change threshold (linear scale) to use to define differential proteins. Default: 2.
A numeric 2-elements vector indicating the range (linear scale) used to define the unresponsive fold changes. Default: c(1/1.1, 1.1).
Numeric value indicating the p.adjusted threshold to apply to the differential analyses. Default: 0.05.
String indicating the method to use to correct the p-values. One among: "holm", "hochberg", "hommel", "bonferroni", "BH", "BY", "fdr", "none", "fdrtool". Default: BH.
Logical value indicating whether the protein means should be moderated towards a common location prior. Default: TRUE.
Logical value indicating whether the protein variances should be moderated towards a common variance prior (empirical Bayes, as in limma). Default: TRUE.
Integer indicating the minimum number of observed (non-missing) values that a protein must have, over all the samples included in the model, to be fitted. Proteins below this threshold, as well as the proteins never quantified in any of the samples of a specific contrast, are kept in the results table with NA statistics. Default: 1.
Integer indicating the number of proteins to use to estimate the hyper-parameters (dropout curves, location and variance priors). When NULL all the proteins are used. Rows are shuffled before the subsampling. Default: NULL.
Integer indicating the maximum number of iterations of the hyper-parameter estimation loop. Default: 20.
Integer indicating the number of times a new subsample of proteins is drawn when the estimation of the hyper-parameters fails. Used only when n.subsample is not NULL; the last attempt is always performed on all the proteins available for the contrast. Default: 3.
An object generated by missingness.diagnostic used to check the censoring assumption. When NULL (default) the diagnostic is computed internally if check.missingness = TRUE. Default: NULL.
Logical value indicating whether the missingness pattern should be checked before fitting the model. Default: TRUE.
String indicating the color to use for up-regulated proteins in the plots. Default: "indianred".
String indicating the color to use for up-regulated proteins in the plots. Default: "steelblue".
String indicating the color to use for unresponsive proteins in the plots. Default: "purple".
String indicating the color to use for null proteins in the plots. Default: "gray".
String indicating which type of counts should be used. One among: 'raw', 'normalized', 'norm', 'randomized', 'random', 'imputed', 'imp'. Default: "normalized".
Numeric value to use as seed when n.subsample is defined. Default: NULL (random seed).
Logical value indicating whether the messages of the model fitting should be printed. Default: FALSE.
Logical value to indicate whether overwrite analyses already generated. Default: FALSE.
An object of class DEprot.analyses
Differently from the other differential functions of the package, this one must be applied to counts that
still contain the missing values (normalized, not imputed): a missing value is not replaced but read as
the statement "this protein was below the detection limit in this sample", exactly as a left-censored
observation in a survival model. The probability of observing a value at a given intensity is described by
a sigmoidal dropout curve estimated separately for each sample, and the group means are inferred on the
mixture of measured and censored observations. The uncertainty on the censored values is therefore kept
inside the model and propagated to lfcSE, so that a protein detected in one replicate out of three
receives a wide standard error instead of the artificially small variance produced by the imputation of
values drawn from the bottom of the distribution.
The model is fitted on all the samples sharing the metadata column of the contrast (design
~ 0 + <metadata.column>, plus the replicate when include.rep.model = TRUE), and the contrast
is tested afterwards with a Wald test: the behaviour is the same as diff.analyses.limma, and it
means that the samples of the other groups of the same column contribute to the estimation of the
hyper-parameters and to the residual degrees of freedom. Contrasts defined on the same metadata column
share a single fit whenever they are computed on the same set of proteins.
The proteins used for the fit are defined contrast by contrast: a protein that has not been quantified in
any of the samples of the two groups compared carries no information on that specific contrast, the
censored likelihood becomes flat and the optimizer fails. These proteins, together with the ones below
min.detected, are excluded from the fit and re-introduced in the results table with NA fold
change, p-value and differential status, so that the output always keeps one row per protein, in the order
of the counts table. Being untested, they are not counted in n.diff and are not drawn in the volcano
and MA plots. Proteins detected in only one of the two groups are instead retained, since that asymmetry is
exactly what the dropout-aware model is meant to describe.
The log2(FoldChange) is the contrast estimate of the censored model, and not the difference between
the observed means: for a protein missing in one of the two groups the two quantities differ substantially,
since the model extrapolates below the detection limit instead of using the values that happened to be
observed. The log2.mean.<group> columns report the model-based (censoring-aware) group means,
computed as the average of the fitted values over the samples of each group; their difference coincides
with the reported log2(FoldChange) for balanced designs, but not necessarily when
include.rep.model = TRUE on an unbalanced design. The sd.<group> and sem.<group>
columns are instead purely descriptive and are computed on the observed values only, so they are affected
by the survivor bias of the missing values.
Three extra columns are appended at the end of the results table: n.detected.<group>, the number of
samples in which the protein was actually quantified in each group, and n.approx, the effective
sample size estimated by proDA over all the samples of the model, which counts how much information
the censored observations carry (a value close to n.obs means that the missing values are not
informative).
When set to "fdrtool", the moderated t-statistics are passed to fdrtool::fdrtool (as normal
deviates, statistic = "normal"), which estimates the proportion of null features and an empirical
null: the resulting tail-area-based q-values are stored in the padj column, while the local false
discovery rate (lfdr) is added as an extra column of the results table. If fewer than 200 proteins
with a finite statistic are available for a contrast, or the fit fails, a Benjamini-Hochberg correction is
applied for that contrast (a warning is raised).
The censoring assumption is verified before fitting the model: when missingness.diagnostic finds no
significant dependency between the intensity of a protein and its probability of being missing, the data
are not left-censored and there is nothing to gain from this model over diff.analyses.limma applied
to imputed counts (a warning is raised, the analyses are performed anyway). The dropout curve position
estimated by proDA for each sample is stored in the results list and is the per-sample equivalent of
the LOD50 returned by the diagnostic: strongly diverging values indicate that one or more samples
have a detection limit of their own and should be inspected.
# Dropout-aware analyses: the counts must still contain the missing values
dpo <- diff.analyses.proDA(DEprot.object = DEprot::test.toolbox$dpo.norm,
contrast.list = list(c("condition", "FBS", "6h.DMSO"),
c("condition", "6h.10nM.E2", "6h.DMSO")),
linear.FC.th = 1.5)
#> Counts available: raw, normalized | counts used: normalized.
#> Error in which(unlist(lapply(X = fit.cache, FUN = function(x) { identical(x$id, cache.id) & identical(sort(x$prot), sort(prot.fit))}))): argument to 'which' is not logical
# Including the replicate as a blocking factor
dpo <- diff.analyses.proDA(DEprot.object = DEprot::test.toolbox$dpo.norm,
contrast.list = list(c("condition", "FBS", "6h.DMSO"),
c("condition", "6h.10nM.E2", "6h.DMSO")),
include.rep.model = TRUE,
replicate.column = "replicate",
linear.FC.th = 1.5)
#> Counts available: raw, normalized | counts used: normalized.
#> Error in which(unlist(lapply(X = fit.cache, FUN = function(x) { identical(x$id, cache.id) & identical(sort(x$prot), sort(prot.fit))}))): argument to 'which' is not logical