fit_Kd: Fitting sigmoidal dose-response curves
In benjbuch/summerrband: Analyse Gel-shift Assays

Description Usage Arguments Details Examples

View source: R/fit_affinity.R

Fits a Kd curve for the equilibrium [RL] <=> [R] + [L], where [RL] is measured by a proxy variable (response) in dependence of the concentration [L] (dose).

fit_Kd(
  x,
  formula,
  R0 = NaN,
  include_hill = FALSE,
  limits_lower = c(-Inf, +Inf),
  limits_upper = c(-Inf, +Inf),
  limits_hill = c(-Inf, +Inf),
  limits_K_d = c(0, 1000),
  start = "auto",
  ...,
  FUN = "minpack.lm::nlsLM"
)

`x`	A data frame to be evaluated.
`formula`	A formula using variables of `x`; given as `RL ~ L0`, where `RL` is the variable that describes the (normalized) concentration of the complex [RL]_0, and `L0` the concentration of titrated species [L]_0.
`R0`	The total concentration of the species that is held constant during the titration. If unknown (`NaN`), the constant ligand approximation is used for fitting by default. See Details.
`include_hill`	If `TRUE`, a Hill coefficient is applied onto `R`; not always appropiate.
`limits_lower, limits_upper, limits_hill, limits_K_d`	Upper and lower bounds for the paramters to be fitted.
`start`	A named list or named numeric vector of starting estimates. Guessed from the data if `"auto"`; not passed if `NULL`.
`...`	Other arguments passed to the fitting algorithm `FUN`.
`FUN`	The function used as fitting algorithm; must return a non-linear model. Can be an unquoted function, but must be a character if prefixed as "package::function".

If R0 is a numeric value, the total concentration of the species that is not titrated in the assay, the following quadratic equation is fitted:

[RL] = 1/(2 [R]_0) * ([R]_0 + [L]_0 + K_d - (([R]_0 + [L]_0 + K_d)^2 - 4 [R]_0 [L]_0)^(1/2))

else

[RL] = ([R]_0 [L]_0) / (K_d + [L]_0)

The latter makes use of the approximation [L]_0 \approx [L], which is valid only if [R]_0 << K_d over the inflection point of the curve. Else, the data is biased by ligand depletion.

In both cases, the fit can be normalized on the [RL] axis using the fitted parameters lower and upper as lower + (upper - lower) * [RL].

The choice of the designations R (as mnemonic for "receptor") and L (for "ligand") are arbitrary; the fitting works also in reversed scenarios.

Choice of limits

By default, the search space for the parameters, except for K_d, is unconstraint (c(-Inf, +Inf)). This is because the fitting algorithm may need to pass via "unreasonable" values during its search for the local/global minimum.

Choice of the fitting algorithm

By default, minpack.lm::nlsLM is used as an alternative to base R's nls. You may consider using nls.multstart::nls_multstart with additional arguments via ... to probe the dependence on different starting parameters. (See Examples.)

Choice of starting values

Starting values are guessed from the data if "auto"; by default: lower is zero, upper is the maximum of [RL] (LHS of formula), hill, if present, is 1, K_d is value of [L]_0 (RHS of formula) which is closest to the half-maximal value of [RL].

You must set start = NULL with certain FUN, e.g. when using FUN = "nls.multstart::nls_multstart".

library(dplyr)
assay_file <- system.file("extdata", "gel_01.txt", package = "summerrband")
assay_data <- iqtl_meta(assay_file, list(conc = c(2^seq(10, 0), 0)))

assay_data %>%
  filter(band_id == "band_1") %>%
  fit_Kd(vol_frac ~ conc, include_hill = FALSE, R0 = 2) %>%
  broom::tidy()

assay_data %>%
  group_by(band_id) %>%
  summerr::model_cleanly_groupwise(fit_Kd, formula = vol_frac ~ conc)

library(ggplot2)

assay_data %>%
  group_by(band_id) %>%
  summerr::model_cleanly_groupwise(fit_Kd, formula = vol_frac ~ conc) %>%
  summerr::model_display(color = band_id) +
  scale_x_log10()

# Other fitting algorithms

assay_data %>%
  filter(band_id == "band_1") %>%
  fit_Kd(vol_frac ~ conc, include_hill = FALSE, R0 = 2, FUN = "nls", algorithm = "port") %>%
  broom::tidy()

assay_data %>%
  filter(band_id == "band_1") %>%
  fit_Kd(vol_frac ~ conc, include_hill = FALSE, R0 = 2,
    FUN = "nls.multstart::nls_multstart", start = NULL,  ## important!
    iter = 500,
    start_lower = c(lower = 0,   upper = 0.5, K_d = 1e-3),
    start_upper = c(lower = 0.5, upper = 1.0, K_d = 1e3)) %>%
  broom::tidy()