Nothing
#' @importFrom magrittr %>%
#' @importFrom magrittr %<>%
scramblase_assay_calculations <- function(
x,
scale_to,
ppr_scale_factor = 0.65,
generation_of_algorithm = 2,
force_through_origin = TRUE,
split_by_experiment = TRUE,
n_averaging = 10,
r_bar = 88,
sigma_r_bar = 28){
# Set parameters ----------------------------------------------------------
nlsControl <- list(
minFactor = 1/20480,
maxit = 100)
formulae <- list(
pre_fit = list(
"TRUE" = y ~ b * ( 1 - exp( -x / a ) ),
"FALSE" = y ~ b - c * exp( -x / a )),
first_generation_algorithm = list(
"TRUE" = y ~ b * ( 1 - exp( -x / a ) ),
"FALSE" = y ~ b - c * exp( -x / a )),
second_generation_algorithm = list(
"TRUE" = y ~ b * ( 1 - ( 1 / sqrt( 1 + sigma_r_bar^2 * a * x ) ) * exp( ( -a * r_bar^2 / 2 * x ) / ( 1 + sigma_r_bar^2 * a * x ) ) ),
"FALSE" = y ~ b - c * ( ( 1 / sqrt( 1 + sigma_r_bar^2 * a * x ) ) * exp( ( -a * r_bar^2 / 2 * x ) / ( 1 + sigma_r_bar^2 * a * x ) ) )))
generation_of_algorithm %<>%
switch(
"1" = "first_generation_algorithm",
"2" = "second_generation_algorithm")
# Parsing spectra ---------------------------------------------------------
spectralData <- x %>%
magrittr::extract2("Path") %>%
lapply(parse_fluorimeter_output, adjust = TRUE)
# Read out data -----------------------------------------------------------
# What spectral time windows to extract?
## Just checking ...
maxAcquisitionTime <- x %>%
magrittr::extract2("Timepoint of Measurement (s)") %>%
max(na.rm = TRUE)
# Average over 10 values before dithionite addition for activity baseline
x[["Baseline Fluorescence"]] <- spectralData %>%
vapply(
function(spectral_data_i){
spectral_data_i %>%
magrittr::extract(.$Time.in.sec < 0, ) %>%
utils::tail(n_averaging) %>%
magrittr::extract2("Fluorescence.Intensity") %>%
stats::median(na.rm = TRUE)
},
1)
# Average over last 10 values (in common time range) for activity
x[["Minimum Fluorescence"]] <- mapply(
function(spectral_data_i, timepoint_of_measurement_s)
{
# TODO: this behaviour matches the previous implementation, but using <=
# seems slightly more intuitive
spectral_data_i %>%
magrittr::extract(.$Time.in.sec < timepoint_of_measurement_s, ) %>%
utils::tail(n_averaging) %>%
magrittr::extract2("Fluorescence.Intensity") %>%
stats::median(na.rm = TRUE)
},
spectralData,
x$"Timepoint of Measurement (s)"
)
# Apply volume correction factors as needed
volumeCorrectionFactor <- x[["Fluorescence Assay Vol. with DT (ul)"]] /
x[["Fluorescence Assay Vol. w/o DT (ul)"]]
x[["Minimum Fluorescence, Volume Corrected"]] <- x[["Minimum Fluorescence"]] *
volumeCorrectionFactor
## Calculate activity reduction
x[["Fluorescence Reduction"]] <- 1 -
x[["Minimum Fluorescence, Volume Corrected"]] /
x[["Baseline Fluorescence"]]
# Generate PPR v.s P(>=1 Scramblase/Liposome) -------------------------------
# Split by Experiment
x[["CombinedId"]] <- paste(
x[["Experimental Series"]],
x[["Experiment"]],
sep = "_")
inputListFromX <- split(
x,
x[["CombinedId"]])
# Calculations A
processedListFromX <- inputListFromX %>%
lapply(
function(z){
## Ensure that there's a data point with liposomes ONLY as a unique
## reference point
indexOfLiposomesOnlyData <- z %>%
magrittr::extract2("Protein Reconstituted (mg)") %>%
magrittr::equals(0) %>%
which()
if(length(indexOfLiposomesOnlyData) == 0){
stop("Experimental series '",unique(z[["Experimental Series"]]),"' does
not have the required liposomes-ONLY ('Extract Volume (ul)' of '0')
data point. Aborting.")
}
if(length(indexOfLiposomesOnlyData) > 1){
stop("Experimental series '",unique(z[["Experimental Series"]]),"' has
more than one liposomes-ONLY ('Extract Volume (ul)' of '0') data
point.")
}
##> End-point fluorescence reduction data from flippase activity assays
##> were obtained for proteoliposomes generated over a range of PPR
##> values.
##> The data were transformed according to the formula
##> P(>=1 flippase) = (y - yo)/(yMax - yo),
##> where yo is the percent reduction obtained with liposomes, yMax is the
##> maximum percentage reduction observed and p is the probability that a
##> particular vesicle in the ensemble is 'flippase-active', i.e it
##> possesses at least one flippase.
## Calculate the relative fluorescence reduction
z[["Relative Fluorescence Reduction"]]<- z[["Fluorescence Reduction"]] -
z[["Fluorescence Reduction"]][indexOfLiposomesOnlyData]
## Calculate PPR
z[["Protein per Phospholipid (mg/mmol)"]] <- z %>%
calculate_ppr(ppr_scale_factor = ppr_scale_factor)
## Calculate p>=1Scramblase/Liposome
y <- z %>%
magrittr::extract2("Relative Fluorescence Reduction")
y0 <- y %>%
magrittr::extract(indexOfLiposomesOnlyData)
if(scale_to == "model"){
fit_prep <- z %>%
fit_prep(y = "Relative Fluorescence Reduction")
fitStart <- if(force_through_origin){
list(
a = fit_prep[["estimated_a"]],
b = fit_prep[["estimated_b"]])
} else {
start = list(
a = fit_prep[["estimated_a"]],
b = fit_prep[["estimated_b"]],
c = fit_prep[["estimated_b"]])
}
rMod <- minpack.lm::nlsLM(
formula = formulae[["pre_fit"]][[as.character(force_through_origin)]],
data = fit_prep[["fit_set"]],
start = fitStart,
control = nlsControl)
yMax <- rMod %>%
stats::coef() %>%
magrittr::extract2("b")
} else {
yMax <- z %>%
magrittr::extract2("Relative Fluorescence Reduction") %>%
max(na.rm = TRUE)
}
z[["Probability >= 1 Scramblase in Vesicle"]] <- (y-y0)/(yMax-y0)
return(z)
})
# Is separate handling of experiments required?
if(!split_by_experiment){
processedListFromX <- processedListFromX %>%
plyr::rbind.fill()
processedListFromX <- processedListFromX %>%
split(processedListFromX$"Experimental Series")
}
# Calculatons B
processedListFromX %<>%
lapply(
function(z){
fit_prep <- z %>%
fit_prep(y = "Probability >= 1 Scramblase in Vesicle")
##> Reference: Menon 2011, "Opsin is a phospholipid flippase",
##> supplementary material 01, p2-3, section
##> "Protein-dependence of lipid flipping; analysis of Figure 1f."
##> https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3057128/#SD1
##>
##> The dependence of P(>=1 flippase) on PPR was analyzed as follows.
##> Definitions:
##> f, number of flippases used for reconstitution
##> v, number of vesicles(
##> m, number of flippases per vesicle (=f/v)
##> PPR, mg protein per mmol phospholipid
##> To calculate P(>=1 flippase) as a function of the PPR, we assume that
##> reconstitution of opsin/rhodopsin molecules into vesicles occurs
##> independently and that the vesicles are identical and may have more
##> than one flippase. The probability that a flippase will be
##> reconstituted into a particular vesicle is therefore 1/v. On
##> reconstituting f flippases into an ensemble of v vesicles, the
##> probability P(k) that a particular vesicle contains k flippases is
##> given by the binomial formula:
##>
##> P(k) = C(f,k)*(1/v)^k*(1-1/v)^(f-k)
##>
##> where C(f,k) is the number of combinations of k elements chosen
##> from a set of size f.
##>
##> Because f and v are both large, it is convenient to use the Poisson
##> approximation:
##>
##> P(k) = (m^k/k!)e^(-m), where m = f/v is the average number of
##> flippases per vesicle
##>
##> The probability of a particular vesicle having no flippases is
##> P(0) = e^(-m); therefore, the probability that a vesicle has one or
##> more flippases, i.e, is active in the flippase assay, is
##>
##> P(>=1) = 1-P(0) = 1 - e^(-m)
##>
##> The average number of flippases per vesicle, m, is proportional to PPR
##> and can be written as m = PPR/alpha, where alpha is a constant with units of
##> mg/mmol.
##> Thus,
##>
##> P(>=1) = 1 - e^(-m) = 1 - exp(-PPR/alpha)
##>
##> The mono-exponential fit constant for a graph of P(>=1) vs PPR is
##> alpha mg/mmol; at this PPR value, m = 1 and ~63% (= 1 - e^(-1)) of the
##> vesicles in the population possess >=1 flippase.
## Fit a monoexponential curve to the data
rMod <- minpack.lm::nlsLM(
formula = formulae[[generation_of_algorithm]][[as.character(force_through_origin)]],
data = fit_prep[["fit_set"]],
start = if(force_through_origin){
list(a = fit_prep[["estimated_a"]], b = 1)
} else {
list(a = fit_prep[["estimated_a"]], b = 1, c = 1)
},
control = nlsControl)
gc()
z[["Fit Constant (a)"]] <- stats::coef(rMod)[["a"]]
output <- list(
Raw = z,
FitObject = rMod)
## Generate data to plot the results of the fit
xPredictedFromFit <- seq(
from = min(
z[["Protein per Phospholipid (mg/mmol)"]],
na.rm=TRUE),
to = max(
z[["Protein per Phospholipid (mg/mmol)"]],
na.rm=TRUE),
length.out=200)
yPredictedFromFit <- stats::predict(
rMod,
newdata = list(x = xPredictedFromFit))
output[["Fit"]] <- data.frame(
"Protein per Phospholipid (mg/mmol)"=xPredictedFromFit,
"Probability >= 1 Scramblase in Vesicle"=yPredictedFromFit,
"Experimental Series"=unique(z[["Experimental Series"]]),
"Experiment"=unique(z[["Experiment"]]),
check.names=FALSE,
stringsAsFactors=FALSE)
## Return
return(output)
})
# Output ------------------------------------------------------------------
return(processedListFromX)
}
calculate_ppr <- function(x, ppr_scale_factor = 0.65){
ppr <- x[["Protein Reconstituted (mg)"]]/x[["Lipid in Reconstitution (mmol)"]]
if(!is.null(ppr_scale_factor)){
ppr <- ppr/ppr_scale_factor
}
return(ppr)
}
# Helper Functions --------------------------------------------------------
fit_prep <- function(z, y = c("Relative Fluorescence Reduction", "Probability >= 1 Scramblase in Vesicle")){
y %<>%
match.arg(
choices = c(
"Relative Fluorescence Reduction",
"Probability >= 1 Scramblase in Vesicle"),
several.ok = FALSE)
fit_set <- z %>%
magrittr::extract(
c("Protein per Phospholipid (mg/mmol)",
y)) %>%
magrittr::set_names(c("x", "y"))
### Determine a sensible start point for 'a'
pointSixY <- fit_set %>%
magrittr::extract2("y") %>%
max(na.rm = TRUE) %>%
magrittr::multiply_by(0.6)
estimated_a <- fit_set %>%
magrittr::extract2("x") %>%
magrittr::extract(
fit_set %>%
magrittr::extract2("y") %>%
magrittr::subtract(pointSixY) %>%
abs() %>%
which.min())
estimated_b <- z %>%
magrittr::extract2(y) %>%
max(na.rm = TRUE)
list(
fit_set = fit_set,
estimated_a = estimated_a,
estimated_b = estimated_b) %>%
return()
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.