Function for fitting four-parameter dose response curves for each group (precursor, peptide or protein). In addition it can annotate data based on completeness, the completeness distribution and statistical testing using ANOVA. Filtering by the function is only performed based on completeness if selected.
fit_drc_4p(
data,
sample,
grouping,
response,
dose,
filter = "post",
replicate_completeness = 0.7,
condition_completeness = 0.5,
n_replicate_completeness = NULL,
n_condition_completeness = NULL,
complete_doses = NULL,
anova_cutoff = 0.05,
correlation_cutoff = 0.8,
log_logarithmic = TRUE,
include_models = FALSE,
retain_columns = NULL
)Arguments
- data
a data frame that contains at least the input variables.
- sample
a character column in the
datadata frame that contains the sample names.- grouping
a character column in the
datadata frame that contains the precursor, peptide or protein identifiers.- response
a numeric column in the
datadata frame that contains the response values, e.g. log2 transformed intensities.- dose
a numeric column in the
datadata frame that contains the dose values, e.g. the treatment concentrations.- filter
a character value that can either be
"pre","post"or"none". The data is annotated for completeness, ANOVA significance and the completeness distribution along the doses ("pre"and"post"). The combined output of this filtering step can be found in thepassed_filtercolumn and depends on the cutoffs provided to the function. Note that this is only an annotation and nothing is removed from the output. If"pre"is selected then, in addition to the annotation, the data is filtered for completeness based on the condition completeness prior to the curve fitting and ANOVA calculation and p-value adjustment. This has the benefit that less curves need to be fitted and that the ANOVA p-value adjustment is done only on the relevant set of tests. If"none"is selected the data will be neither annotated nor filtered.- replicate_completeness
![[Deprecated]](figures/lifecycle-deprecated.svg)
please use n_replicate_completenessinstead. A numeric value which similar tocompletenss_MARof theassign_missingnessfunction sets a threshold for the completeness of data. In contrast toassign_missingnessit only determines the completeness for one condition and not the comparison of two conditions. The threshold is used to calculate a minimal degree of data completeness. The value provided to this argument has to be between 0 and 1, default is 0.7. It is multiplied with the number of replicates and then adjusted downward. The resulting number is the minimal number of observations that a condition needs to have to be considered "complete enough" for thecondition_completenessargument.- condition_completeness
- please use
n_condition_completenessinstead. A numeric value which determines how many conditions need to at least fulfill the "complete enough" criteria set withreplicate_completeness. The value provided to this argument has to be between 0 and 1, default is 0.5. It is multiplied with the number of conditions and then adjusted downward. The resulting number is the minimal number of conditions that need to fulfill thereplicate_completenessargument for a peptide to pass the filtering. - n_replicate_completeness
a numeric value that defines the minimal number of observations that a condition (concentration) needs to have to be considered "complete enough" for the
n_condition_completenessargument. E.g. if each concentration has 4 replicates this argument could be set to 3 to allow for one replicate to be missing for the completeness criteria.- n_condition_completeness
a numeric value that defines the minimal number of conditions that need to fulfill the
n_replicate_completenessargument for a feature to pass the filtering. E.g. if an experiment has 12 concentrations, this argument could be set to 6 to define that at least 6 of 12 concentrations need to make the replicate completeness cutoff.- complete_doses
an optional numeric vector that supplies all the actually used doses (concentrations) to the function. Usually the function extracts this information from the supplied data. However, for incomplete datasets the total number of assumed doses might be wrong. Therefore, it becomes important to provide this argument when the dataset is small and potentially incomplete. This information is only used for the missing not at random (MNAR) estimations.
- anova_cutoff
a numeric value that specifies the ANOVA adjusted p-value cutoff used for data filtering. Any fits with an adjusted ANOVA p-value bellow the cutoff will be considered for scoring. The default is
0.05.- correlation_cutoff
a numeric value that specifies the correlation cutoff used for data filtering. Any fits with a correlation above the cutoff will be considered for scoring.
- log_logarithmic
a logical value that indicates if a logarithmic or log-logarithmic model is fitted. If response values form a symmetric curve for non-log transformed dose values, a logarithmic model instead of a log-logarithmic model should be used. Usually biological dose response data has a log-logarithmic distribution, which is the reason this is the default. Log-logarithmic models are symmetric if dose values are log transformed.
- include_models
a logical value that indicates if model fit objects should be exported. These are usually very large and not necessary for further analysis.
- retain_columns
a vector that specifies columns that should be retained from the input data frame. Default is not retaining additional columns
retain_columns = NULL. Specific columns can be retained by providing their names (not in quotations marks, just like other column names, but in a vector).
Value
If include_models = FALSE a data frame is returned that contains correlations
of predicted to measured values as a measure of the goodness of the curve fit, an associated
p-value and the four parameters of the model for each group. Furthermore, input data for plots
is returned in the columns plot_curve (curve and confidence interval) and plot_points
(measured points). If include_models = TURE, a list is returned that contains:
fit_objects: The fit objects of typedrcfor each group.correlations: The correlation data frame described above
Details
If data filtering options are selected, data is annotated based on multiple criteria.
If "post" is selected the data is annotated based on completeness, the completeness distribution, the
adjusted ANOVA p-value cutoff and a correlation cutoff. Completeness of features is determined based on
the n_replicate_completeness and n_condition_completeness arguments. The completeness distribution determines
if there is a distribution of not random missingness of data along the dose. For this it is checked if half of a
features values (+/-1 value) pass the replicate completeness criteria and half do not pass it. In order to fall into
this category, the values that fulfill the completeness cutoff and the ones that do not fulfill it
need to be consecutive, meaning located next to each other based on their concentration values. Furthermore,
the values that do not pass the completeness cutoff need to be lower in intensity. Lastly, the difference
between the two groups is tested for statistical significance using a Welch's t-test and a
cutoff of p <= 0.1 (we want to mainly discard curves that falsely fit the other criteria but that
have clearly non-significant differences in mean). This allows curves to be considered that have
missing values in half of their observations due to a decrease in intensity. It can be thought
of as conditions that are missing not at random (MNAR). It is often the case that those entities
do not have a significant p-value since half of their conditions are not considered due to data
missingness. The ANOVA test is performed on the features by concentration. If it is significant it is
likely that there is some response. However, this test would also be significant even if there is one
outlier concentration so it should only be used only in combination with other cutoffs to determine
if a feature is significant. The passed_filter column is TRUE for all the
features that pass the above mentioned criteria and that have a correlation greater than the cutoff
(default is 0.8) and the adjusted ANOVA p-value below the cutoff (default is 0.05).
The final list is ranked based on a score calculated on entities that pass the filter.
The score is the negative log10 of the adjusted ANOVA p-value scaled between 0 and 1 and the
correlation scaled between 0 and 1 summed up and divided by 2. Thus, the highest score an
entity can have is 1 with both the highest correlation and adjusted p-value. The rank is
corresponding to this score. Please note, that entities with MNAR conditions might have a
lower score due to the missing or non-significant ANOVA p-value. If no score could be calculated
the usual way these cases receive a score of 0. You should have a look at curves that are TRUE
for dose_MNAR in more detail.
If the "pre" option is selected for the filter argument then the data is filtered for completeness
prior to curve fitting and the ANOVA test. Otherwise annotation is performed exactly as mentioned above.
We recommend the "pre" option because it leaves you with not only the likely hits of your treatment, but
also with rather high confidence true negative results. This is because the filtered data has a high
degree of completeness making it unlikely that a real dose-response curve is missed due to data missingness.
Please note that in general, curves are only fitted if there are at least 5 conditions with data points present to ensure that there is potential for a good curve fit. This is done independent of the selected filtering option.
Examples
# \donttest{
# Load libraries
library(dplyr)
set.seed(123) # Makes example reproducible
# Create example data
data <- create_synthetic_data(
n_proteins = 2,
frac_change = 1,
n_replicates = 3,
n_conditions = 8,
method = "dose_response",
concentrations = c(0, 1, 10, 50, 100, 500, 1000, 5000),
additional_metadata = FALSE
)
# Perform dose response curve fit
drc_fit <- fit_drc_4p(
data = data,
sample = sample,
grouping = peptide,
response = peptide_intensity_missing,
dose = concentration,
n_replicate_completeness = 2,
n_condition_completeness = 5,
retain_columns = c(protein, change_peptide)
)
glimpse(drc_fit)
#> Rows: 45
#> Columns: 19
#> $ rank <dbl> 1, 2, 3, 4, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
#> $ protein <chr> "protein_2", "protein_2", "protein_1", "protein_2", …
#> $ change_peptide <lgl> TRUE, TRUE, TRUE, TRUE, TRUE, FALSE, FALSE, FALSE, F…
#> $ peptide <chr> "peptide_2_1", "peptide_2_3", "peptide_1_1", "peptid…
#> $ hill_coefficient <dbl> 2.226694e+00, 3.280387e+00, 1.909346e+00, -3.310549e…
#> $ min_model <dbl> 14.87823, 12.46267, 15.98568, 19.14721, 17.38775, 14…
#> $ max_model <dbl> 19.12933, 16.18910, 17.09790, 19.80934, 17.91039, 15…
#> $ ec_50 <dbl> 2.453723e+02, 1.417578e+02, 9.518399e+01, 6.293053e+…
#> $ correlation <dbl> 0.9986256, 0.9938863, 0.9668006, 0.9280152, 0.769789…
#> $ pval <dbl> 1.131609e-29, 8.423311e-15, 1.595858e-14, 6.625673e-…
#> $ plot_curve <list> [<tbl_df[100 x 4]>], [<tbl_df[100 x 4]>], [<tbl_df[…
#> $ plot_points <list> [<tbl_df[24 x 2]>], [<tbl_df[16 x 2]>], [<tbl_df[24…
#> $ enough_conditions <lgl> TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE…
#> $ dose_MNAR <lgl> FALSE, TRUE, FALSE, FALSE, FALSE, FALSE, FALSE, FALS…
#> $ anova_pval <dbl> 3.666126e-20, 1.471594e-05, 1.966510e-08, 2.858086e-…
#> $ anova_adj_pval <dbl> 1.613095e-18, 1.618753e-04, 4.326322e-07, 4.191860e-…
#> $ anova_significant <lgl> TRUE, TRUE, TRUE, TRUE, TRUE, FALSE, FALSE, FALSE, F…
#> $ passed_filter <lgl> TRUE, TRUE, TRUE, TRUE, FALSE, FALSE, FALSE, FALSE, …
#> $ score <dbl> 1.00000000, 0.46644068, 0.37304758, 0.05860592, NA, …
head(drc_fit, n = 10)
#> # A tibble: 10 × 19
#> rank protein change_peptide peptide hill_coefficient min_model max_model
#> <dbl> <chr> <lgl> <chr> <dbl> <dbl> <dbl>
#> 1 1 protein_2 TRUE peptide_… 2.23 14.9 19.1
#> 2 2 protein_2 TRUE peptide_… 3.28 12.5 16.2
#> 3 3 protein_1 TRUE peptide_… 1.91 16.0 17.1
#> 4 4 protein_2 TRUE peptide_… -3.31 19.1 19.8
#> 5 NA protein_1 TRUE peptide_… -4.79 17.4 17.9
#> 6 NA protein_1 FALSE peptide_… 2.71 14.8 15.2
#> 7 NA protein_1 FALSE peptide_… -0.0110 15.3 14.7
#> 8 NA protein_1 FALSE peptide_… 3.59 16.7 17.2
#> 9 NA protein_1 FALSE peptide_… -0.0000174 16.1 15.4
#> 10 NA protein_1 FALSE peptide_… -1.33 16.2 16.5
#> # ℹ 12 more variables: ec_50 <dbl>, correlation <dbl>, pval <dbl>,
#> # plot_curve <list>, plot_points <list>, enough_conditions <lgl>,
#> # dose_MNAR <lgl>, anova_pval <dbl>, anova_adj_pval <dbl>,
#> # anova_significant <lgl>, passed_filter <lgl>, score <dbl>
# }