diff options
Diffstat (limited to 'CakePlot.R')
| -rw-r--r--[-rwxr-xr-x] | CakePlot.R | 290 |
1 files changed, 147 insertions, 143 deletions
diff --git a/CakePlot.R b/CakePlot.R index 1b80a4b..da69ea8 100755..100644 --- a/CakePlot.R +++ b/CakePlot.R @@ -1,8 +1,8 @@ #$Id$ # Generates fitted curves so the GUI can plot them # Based on code in IRLSkinfit -# Modifications developed by Tessella for Syngenta: Copyright (C) 2011-2016 Syngenta -# Tessella Project Reference: 6245, 7247, 8361, 7414 +# Modifications developed by Tessella for Syngenta: Copyright (C) 2011-2020 Syngenta +# Tessella Project Reference: 6245, 7247, 8361, 7414, 10091 # # The CAKE R modules are free software: you can redistribute it and/or modify # it under the terms of the GNU General Public License as published by @@ -19,21 +19,22 @@ CakeDoPlot <- function(fit, xlim = range(fit$data$time), ...) { + original_fit <- fit t.map.names <- names(fit$map) metabolites <- grep("[A-Z]\\d", t.map.names, value = TRUE) t.map.rest <- setdiff(t.map.names, metabolites) # Generate the normal graphs. final <- CakeSinglePlot(fit, xlim) - final_scaled <- final - if(length(metabolites) > 0){ - for(i in 1:length(metabolites)) - { + if (length(metabolites) > 0) { + for (i in 1:length(metabolites)) { metabolite <- metabolites[i] - if (names(fit$map[[metabolite]])[1] != "SFO"){ - # We do not support these curves for models other than SFO (for now; see CU-79). + fit_method = names(fit$map[[metabolite]])[1] + if (fit_method != "SFO" && fit_method != "DFOP") { + # Only SFO and DFOP are current compatible with the process below + # so skip this metabolite. next } @@ -47,72 +48,49 @@ CakeDoPlot <- function(fit, xlim = range(fit$data$time), ...) metabolite0Parameter <- metabolite - if (!(metabolite %in% names(odeini))){ + if (!(metabolite %in% names(odeini))) { metabolite0Parameter <- paste0(metabolite, "_0") } - if (odeini[[metabolite0Parameter]] != 0){ - # We do not support these curves where the initial concentration is non-zero (for now; see CU-79). - next + # If the metabolite is DFOP, we should find a k2_ decay variable, + # otherwise we can fall back to the SFO k_ variant + decay_var <- paste("k2", metabolite, sep="_") + if (!(decay_var %in% names(fit$par)) && + !(decay_var %in% rownames(fit$fixed))) { + decay_var <- paste("k", metabolite, sep="_") } - decay_var <- paste("k", metabolite, sep="_") - - # calculate the new ffm (formation factor) and generate the two ffm scale charts - regex <- paste("f_(.+)_to", metabolite, sep="_") - decays = grep(regex, names(fit$par), value = TRUE) - - if (length(decays) != 1){ - # We do not support these curves where there is formation from more than 1 compartment (for now; see CU-79). - next - } - - ffm_fitted <- sum(fit$par[decays]) - normal <- final - ffm_scale <- NULL - - numeric_DT50 <- as.numeric(fit$distimes[metabolite,][["DT50"]]) - - if (is.na(numeric_DT50)){ - # We can't get anywhere without a numeric DT50. - next - } - - # Generate the curve for DT50=1000d and ffm as fitted. - if (decay_var %in% names(fit$par)){ - k_new <- fit$par[decay_var]*numeric_DT50/1000; - fit$par[decay_var]<- k_new - } - else{ + # Generate the curve for DT50=1000 and ffm as fitted. + if (decay_var %in% names(fit$par)) { + fit$par[decay_var] <- log(2)/1000 + } else { # If decay_var was fixed, need to modify it there. - k_new <- fit$fixed[decay_var,]["value"]*numeric_DT50/1000; - fit$fixed[decay_var,]["value"] <- k_new[decay_var,] + fit$fixed[decay_var,"value"] <- log(2)/1000 } - dt50_1000_ffm_fitted <- CakeSinglePlot(fit, xlim)[metabolite] - - naming <- c(names(final), paste(metabolite, "DT50_1000_FFM_FITTED", sep="_")) - normal <- c(final, dt50_1000_ffm_fitted) - names(normal) <- naming - final <- normal - - # Generate the scaled FFM - if(ffm_fitted != 0) - { - ffm_scale <- 1 / ffm_fitted - final_scaled <- final[c("time", metabolite, paste(metabolite, "DT50_1000_FFM_FITTED", sep="_"))] - final_scaled[t.map.rest] <- NULL; - final_frame <- as.data.frame(final_scaled) - new_names <- c(paste(metabolite, "DT50_FITTED_FFM_1", sep="_"), paste(metabolite, "DT50_1000_FFM_1", sep="_")) - names(final_frame) <- c("time", new_names) - final_frame[new_names]<-final_frame[new_names]*ffm_scale; - - cat("<PLOT MODEL START>\n") + new_col_name <- paste(metabolite, "DT50_1000_FFM_FITTED", sep="_") + final[, new_col_name] <- CakeSinglePlot(fit, xlim)[metabolite] - write.table(final_frame, quote=FALSE) + ffm1_fit <- SetFormationFractionsToOne(fit, metabolite) + if (!is.null(ffm1_fit)) { + # It was possible to set the formation fractions to 1, so generate the curves + new_col_name <- paste(metabolite, "DT50_1000_FFM_1", sep="_") + final[, new_col_name] <- CakeSinglePlot(ffm1_fit, xlim)[metabolite] - cat("<PLOT MODEL END>\n") + # Reset the k values back to their originals + # Then we can generate the curve for DT50 as fitted and ffm=1 + if (decay_var %in% names(ffm1_fit$par)) { + ffm1_fit$par[decay_var] <- original_fit$par[decay_var] + } else { + ffm1_fit$fixed[decay_var,"value"] <- original_fit$fixed[decay_var,"value"] + } + + new_col_name <- paste(metabolite, "DT50_FITTED_FFM_1", sep="_") + final[, new_col_name] <- CakeSinglePlot(ffm1_fit, xlim)[metabolite] } + + # Finally reset the fit for the next metabolite + fit <- original_fit } } @@ -125,98 +103,124 @@ CakeDoPlot <- function(fit, xlim = range(fit$data$time), ...) # View(final) } -CakeSinglePlot <- function(fit, xlim = range(fit$data$time), scale_x = 1.0, ...) +# Get the immediate parents of a metabolite in a topology +GetParents <- function(topology, metabolite) { - solution = fit$solution - if ( is.null(solution) ) { - solution <- "deSolve" - } - atol = fit$atol - if ( is.null(atol) ) { - atol = 1.0e-6 - } - - fixed <- fit$fixed$value - names(fixed) <- rownames(fit$fixed) - parms.all <- c(fit$par, fixed) - ininames <- c( - rownames(subset(fit$start, type == "state")), - rownames(subset(fit$fixed, type == "state"))) - odeini <- parms.all[ininames] - names(odeini) <- gsub("_0$", "", names(odeini)) - odenames <- c( - rownames(subset(fit$start, type == "deparm")), - rownames(subset(fit$fixed, type == "deparm"))) - odeparms <- parms.all[odenames] - odeini <- AdjustOdeInitialValues(odeini, fit, odeparms) + parents <- NULL + + for (topo_metab in names(topology)) { + if (metabolite %in% topology[[topo_metab]]$to) { + parents <- append(parents, topo_metab) + } + } + + return (parents) +} - outtimes <- seq(0, xlim[2], length.out=CAKE.plots.resolution) * scale_x +# Get all the ancestors of a metabolite in a topology +GetAncestors <- function(topology, metabolite, ancestors = NULL) +{ + for (parent in GetParents(topology, metabolite)) { + if (!(parent %in% ancestors)) { + ancestors <- GetAncestors(topology, parent, append(ancestors, parent)) + } + } + + return (ancestors) +} - # Solve the system - evalparse <- function(string) - { - eval(parse(text=string), as.list(c(odeparms, odeini))) +# Set the formation leading to the metabolite to 1 from each parent +# Returns the updated fit or NULL if the metabolite can't be updated for FFM=1 +SetFormationFractionsToOne <- function(fit, metabolite) +{ + parents <- GetParents(fit$topology, metabolite) + + for (parent in parents) { + if (any(GetAncestors(fit$topology, parent) %in% parents)) { + # Reject metabolites where any parent is also an ancestor of another parent + return (NULL) + } + + # Try to find the formation fraction in the fitted and fixed parameters + ffm_var <- paste("f", parent, "to", metabolite, sep="_") + updated_fit <- SetFfmVariable(fit, ffm_var, 1) + if (!is.null(updated_fit)) { + # Formation fraction was found and set to 1 + fit <- updated_fit + next + } + + # The formation fraction for this metabolite is implicitly 1 minus all the others from this parent + # We need to iterate other formation fractions from the parent and set them to zero + for (other_metabolite in fit$topology[[parent]]$to) { + if (other_metabolite == metabolite) { + # Skip the metabolite itself because we know there isn't a formation fraction for it + next + } + + other_ffm_var <- paste("f", parent, "to", other_metabolite, sep="_") + updated_fit <- SetFfmVariable(fit, other_ffm_var, 0) + + # The updated fit shouldn't be able to be null because otherwise the topology was invalid + # or didn't match the parameters in the fit. If it does, we'll throw an error. + if (is.null(updated_fit)) { + stop(paste("The updated fit was NULL. The formation fraction '", other_ffm_var, "' doesn't exist in the fit.", sep="")) + } + + fit <- updated_fit + } } - if (solution == "analytical") { - parent.type = names(fit$map[[1]])[1] - parent.name = names(fit$diffs)[[1]] - o <- switch(parent.type, - SFO = SFO.solution(outtimes, - evalparse(parent.name), - evalparse(paste("k", parent.name, sep="_"))), - FOMC = FOMC.solution(outtimes, - evalparse(parent.name), - evalparse("alpha"), evalparse("beta")), - DFOP = DFOP.solution(outtimes, - evalparse(parent.name), - evalparse(paste("k1", parent.name, sep="_")), - evalparse(paste("k2", parent.name, sep="_")), - evalparse(paste("g", parent.name, sep="_"))), - HS = HS.solution(outtimes, - evalparse(parent.name), - evalparse("k1"), evalparse("k2"), - evalparse("tb")), - IORE = IORE.solution(outtimes, - evalparse(parent.name), - evalparse(paste("k", parent.name, sep="_")), - evalparse("N")) - ) - out <- cbind(outtimes, o) - dimnames(out) <- list(outtimes, c("time", parent.name)) + + return (fit) +} + +# Set the value of ffm_var in the fit, looking in both fitted and fixed parameters +# Returns the updated fit or NULL if the ffm_var couldn't be found +SetFfmVariable <- function(fit, ffm_var, value) +{ + if (ffm_var %in% names(fit$par)) { + # The ffm_var was found as a fitted parameter + fit$par[ffm_var] <- value + return (fit) } - if (solution == "eigen") { - coefmat.num <- matrix(sapply(as.vector(fit$coefmat), evalparse), - nrow = length(odeini)) - e <- eigen(coefmat.num) - c <- solve(e$vectors, odeini) - f.out <- function(t) { - e$vectors %*% diag(exp(e$values * t), nrow=length(odeini)) %*% c - } - o <- matrix(mapply(f.out, outtimes), - nrow = length(odeini), ncol = length(outtimes)) - dimnames(o) <- list(names(odeini), NULL) - out <- cbind(time = outtimes, t(o)) - } - if (solution == "deSolve") { - out <- ode( - y = odeini, - times = outtimes, - func = fit$mkindiff, - parms = odeparms, - atol = atol - ) + + if (ffm_var %in% rownames(fit$fixed)) { + # The ffm_var was found as a fixed parameter + fit$fixed[ffm_var,"value"] <- value + return (fit) + } + + # ffm_var was not found in the fit + return (NULL) +} + +CakeSinglePlot <- function(fit, xlim = range(fit$data$time), scale_x = 1.0, ...) +{ + solution <- fit$solution + if ( is.null(solution) ) { + solution <- "deSolve" } - out_transformed <- PostProcessOdeOutput(out, fit, atol) + atol <- fit$atol + if ( is.null(atol) ) { + atol <- 1.0e-5 + } + + outtimes <- seq(0, xlim[2], length.out=CAKE.plots.resolution) * scale_x + + # Solve the system + odeResult <- CakeOdeSolve(fit, outtimes, solution, atol) + + odeResult_transformed <- PostProcessOdeOutput(odeResult, fit, atol) # Replace values that are incalculably small with 0. for (compartment.name in names(fit$map)) { - if (length(out_transformed[, compartment.name][!is.nan(out_transformed[, compartment.name])]) > 0) { - out_transformed[, compartment.name][is.nan(out_transformed[, compartment.name])] <- 0 + if (length(odeResult_transformed[, compartment.name][!is.nan(odeResult_transformed[, compartment.name])]) > 0) { + odeResult_transformed[, compartment.name][is.nan(odeResult_transformed[, compartment.name])] <- 0 } - out_transformed[, compartment.name][out_transformed[, compartment.name] < 0] <- 0 + odeResult_transformed[, compartment.name][odeResult_transformed[, compartment.name] < 0] <- 0 } - return(out_transformed) + return (odeResult_transformed) } |