1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
|
#$Id$
# Generates fitted curves so the GUI can plot them
# Based on code in IRLSkinfit
# Modifications developed by Hybrid Intelligence (formerly Tessella), part of
# Capgemini Engineering, for Syngenta, Copyright (C) 2011-2022 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
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
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)
if (length(metabolites) > 0) {
for (i in 1:length(metabolites)) {
metabolite <- metabolites[i]
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
}
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]
metabolite0Parameter <- metabolite
if (!(metabolite %in% names(odeini))) {
metabolite0Parameter <- paste0(metabolite, "_0")
}
# 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="_")
}
# 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.
fit$fixed[decay_var,"value"] <- log(2)/1000
}
new_col_name <- paste(metabolite, "DT50_1000_FFM_FITTED", sep="_")
final[, new_col_name] <- CakeSinglePlot(fit, xlim)[metabolite]
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]
# 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
}
}
cat("<PLOT MODEL START>\n")
write.table(final, quote=FALSE)
cat("<PLOT MODEL END>\n")
# View(final)
}
# Get the immediate parents of a metabolite in a topology
GetParents <- function(topology, metabolite)
{
parents <- NULL
for (topo_metab in names(topology)) {
if (metabolite %in% topology[[topo_metab]]$to) {
parents <- append(parents, topo_metab)
}
}
return (parents)
}
# 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)
}
# 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
}
}
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 (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"
}
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(odeResult_transformed[, compartment.name][!is.nan(odeResult_transformed[, compartment.name])]) > 0) {
odeResult_transformed[, compartment.name][is.nan(odeResult_transformed[, compartment.name])] <- 0
}
odeResult_transformed[, compartment.name][odeResult_transformed[, compartment.name] < 0] <- 0
}
return (odeResult_transformed)
}
|
