Title: | Simulation of the Circadian Clock Gene Network |
---|---|
Description: | A preconfigured simulation workflow for the circadian clock gene network. |
Authors: | Ye Yuan [aut, cre] |
Maintainer: | Ye Yuan <[email protected]> |
License: | GPL (>= 3) |
Version: | 0.0.0.9002 |
Built: | 2025-02-22 16:28:56 UTC |
Source: | https://github.com/yeyuan98/clockSim |
Verify that the two matrices share the same time values and the same states.
.compute_validatePair(res_mat1, res_mat2)
.compute_validatePair(res_mat1, res_mat2)
res_mat1 |
Trajectory matrix 1, first column must be time while others are states. |
res_mat2 |
Trajectory matrix 2, first column must be time while others are states. |
# Internal only; error if check fails.
# Internal only; error if check fails.
This package provides preconfigured circadian clock gene network simulation models based on the Odin simulation engine.
Circadian clock is a paradigm for studying negative-positive feedback mechanisms in Biology. Molecular clock, a cell-autonomous transcriptional-translational feedback, is one of the best known gene network oscillators.
Despite widespread interest on molecular clock research, numerical simulation of the clock remains elusive to many if not most wet-lab biologists, which hindered generation of novel hypotheses and efficient exploration of the clock parameter space.
This package aims to provide a low-friction workflow for molecular clock simulation, thereby unleashing the power of in silico exploration to more.
The mypackage functions ...
Cosine Similarity is used to characterize how state direction of two
trajectory agree with each other.
It computes cosine of the angle between the two states at each time.
compute_cosine(res_mat1, res_mat2)
compute_cosine(res_mat1, res_mat2)
res_mat1 |
Trajectory matrix 1, first column must be time while others are states. |
res_mat2 |
Trajectory matrix 2, first column must be time while others are states. |
Formula:
Cosine similarity of the two trajectory, vector with length = number of time steps.
# Perfect alignment (cos=1) mat1 <- cbind(time = 1:3, state1 = 1:3, state2 = 4:6) mat2 <- cbind(time = 1:3, state1 = 2:4, state2 = 5:7) compute_cosine(mat1, mat2) # c(1, 1, 1) # Orthogonal vectors (cos=0) mat3 <- cbind(time = 1:2, state1 = c(1, 0), state2 = c(0, 1)) mat4 <- cbind(time = 1:2, state1 = c(0, 1), state2 = c(1, 0)) compute_cosine(mat3, mat4) # c(0, 0) # Opposite direction (cos=-1) mat5 <- cbind(time = 1:3, state1 = 1:3) mat6 <- cbind(time = 1:3, state1 = -1:-3) compute_cosine(mat5, mat6) # c(-1, -1, -1)
# Perfect alignment (cos=1) mat1 <- cbind(time = 1:3, state1 = 1:3, state2 = 4:6) mat2 <- cbind(time = 1:3, state1 = 2:4, state2 = 5:7) compute_cosine(mat1, mat2) # c(1, 1, 1) # Orthogonal vectors (cos=0) mat3 <- cbind(time = 1:2, state1 = c(1, 0), state2 = c(0, 1)) mat4 <- cbind(time = 1:2, state1 = c(0, 1), state2 = c(1, 0)) compute_cosine(mat3, mat4) # c(0, 0) # Opposite direction (cos=-1) mat5 <- cbind(time = 1:3, state1 = 1:3) mat6 <- cbind(time = 1:3, state1 = -1:-3) compute_cosine(mat5, mat6) # c(-1, -1, -1)
Performs per-state normalization by dividing state values by: range (max-min) or sd (standard deviation). none means no normalization is performed.
compute_normalize(res_mat, method = c("none", "range", "sd"), ref_stats = NULL)
compute_normalize(res_mat, method = c("none", "range", "sd"), ref_stats = NULL)
res_mat |
Trajectory matrix, first column must be time while others are states. |
method |
State normalization method, one of "none", "range", "sd". |
ref_stats |
Optional list containing pre-computed normalization statistics
(e.g., from a previous call to |
Normalized trajectory matrix with attribute "ref_stats"
containing
normalization parameters.
mat <- cbind(time = 1:10, matrix(runif(30, 0, 10), ncol = 3)) # Range normalization with automatic stat computation norm_mat1 <- compute_normalize(mat, "range") # Reuse stats for another matrix mat2 <- cbind(time = 1:10, matrix(runif(30, 5, 15), ncol = 3)) norm_mat2 <- compute_normalize(mat2, "range", attr(norm_mat1, "ref_stats"))
mat <- cbind(time = 1:10, matrix(runif(30, 0, 10), ncol = 3)) # Range normalization with automatic stat computation norm_mat1 <- compute_normalize(mat, "range") # Reuse stats for another matrix mat2 <- cbind(time = 1:10, matrix(runif(30, 5, 15), ncol = 3)) norm_mat2 <- compute_normalize(mat2, "range", attr(norm_mat1, "ref_stats"))
Supports the following methods: fft
and lomb
. lomb
uses the lomb
package for computation. Please note that lomb
can be resource intensive
and significantly slower than fft
.
compute_period(ts_vector, ts_time = NULL, method = "fft")
compute_period(ts_vector, ts_time = NULL, method = "fft")
ts_vector |
Numeric vector of time series |
ts_time |
Numeric vector of time points (if |
method |
Period calculation method. |
If ts_time
is provided, it is passed to the Lomb-Scargle algorithm for
unevenly sampled data computation. In this case, lomb
method returns
period in unit of ts_time
. ts_time
has no effect on the fft
method as
it requires even time spacing.
If ts_time
is NULL
, assume even time spacing and period unit will be
in the (implicitly provided) time spacing of the time series vector.
Power, SNR, and p-value of the period detection are method-specific:
fft
: power = Spectrum power of the peak. \
SNR = (power of peak)/(median power). p-value is not available (NA).
lomb
: power = LS power of the peak. \
SNR is not available (NA). p-value = LS p-value.
Named vector of length 4 (period, power, snr, p.value).
# Generate a period = 50 sine wave data with some noise (even spacing) n <- 1000 time <- seq(1, n, by = 1) ts_data <- sin(2 * pi * time / 50) + rnorm(n, sd = 0.5) compute_period(ts_data) compute_period(ts_data, method = "lomb") # Uneven sampling of the previous data and run lomb method again s <- sample(1:n, n/3) compute_period(ts_data[s], time[s], method = "lomb")
# Generate a period = 50 sine wave data with some noise (even spacing) n <- 1000 time <- seq(1, n, by = 1) ts_data <- sin(2 * pi * time / 50) + rnorm(n, sd = 0.5) compute_period(ts_data) compute_period(ts_data, method = "lomb") # Uneven sampling of the previous data and run lomb method again s <- sample(1:n, n/3) compute_period(ts_data[s], time[s], method = "lomb")
RMSE is used to characterize how absolute values of two trajectory
agree with each other. There are
normalized and default versions:
compute_rmse(res_mat1, res_mat2, normalize = "none")
compute_rmse(res_mat1, res_mat2, normalize = "none")
res_mat1 |
Trajectory matrix 1, first column must be time while others are states. |
res_mat2 |
Trajectory matrix 2, first column must be time while others are states. |
normalize |
Normalization method, one of "none", "range", "sd". |
Default (value scales with variable ):
Normalized (value normalized by range for each variable ):
(spread is customizable)
For normalization, spread is always computed ONLY from res_mat1
. res_mat2
is normalized using spread computed from mat1
to still retain the property
that RMSE allows comparing absolute values of two trajectories.
In this case, normalization is used to prevent numerically large states from dominating RMSE results, giving a more thorough comparison of trajectories.
RMSE of the two trajectory, vector with length = number of input states.
# Perfect agreement (zero error) mat1 <- cbind(time = 1:3, state1 = c(1, 2, 3), state2 = c(4, 5, 6)) mat2 <- mat1 compute_rmse(mat1, mat2) # Simple error case # state1 = sqrt(mean(c(0,0,0.5)^2)) = 0.29 # state2 = sqrt(mean(c(0,0.2,0)^2)) = 0.12 mat3 <- cbind(time = 1:3, state1 = c(1, 2, 3.5), state2 = c(4, 5.2, 6)) compute_rmse(mat1, mat3) # Normalized example (NRMSE = RMSE / spread) # state1 = state2 = 1/20 = 0.05 mat4 <- cbind(time = 1:3, state1 = c(10, 20, 30), state2 = c(40, 50, 60)) mat5 <- cbind(time = 1:3, state1 = c(11, 21, 31), state2 = c(41, 51, 61)) compute_rmse(mat4, mat5, "range") # c(1/20, 1/20) = c(0.05, 0.05)
# Perfect agreement (zero error) mat1 <- cbind(time = 1:3, state1 = c(1, 2, 3), state2 = c(4, 5, 6)) mat2 <- mat1 compute_rmse(mat1, mat2) # Simple error case # state1 = sqrt(mean(c(0,0,0.5)^2)) = 0.29 # state2 = sqrt(mean(c(0,0.2,0)^2)) = 0.12 mat3 <- cbind(time = 1:3, state1 = c(1, 2, 3.5), state2 = c(4, 5.2, 6)) compute_rmse(mat1, mat3) # Normalized example (NRMSE = RMSE / spread) # state1 = state2 = 1/20 = 0.05 mat4 <- cbind(time = 1:3, state1 = c(10, 20, 30), state2 = c(40, 50, 60)) mat5 <- cbind(time = 1:3, state1 = c(11, 21, 31), state2 = c(41, 51, 61)) compute_rmse(mat4, mat5, "range") # c(1/20, 1/20) = c(0.05, 0.05)
Refer to package documentation on preconfigured models in this package. This function is intended for advanced usage only; normally calling other helper functions will suffice for simulation.
getOdinGen()
getOdinGen()
Named list of Odin generator R6 objects.
# TODO
# TODO
This function is useful for running a large scan of parameter combinations for the same model. Typical use cases are probing stability of an attractor, effect of certain parameters on the system, etc.
grid_scan( model_gen, grid, apply.fn = identity, n.core = 2, custom.export = NULL, ... )
grid_scan( model_gen, grid, apply.fn = identity, n.core = 2, custom.export = NULL, ... )
model_gen |
Odin model generator, see |
grid |
Data frame of the parameter grid. |
apply.fn |
Function to apply before return (e.g., some summary). |
n.core |
Number of cores for parallel computing. |
custom.export |
Names of additional variables used by apply.fn. |
... |
Model |
Grid is a data frame whose columns are model parameters.
See model$contents()
for tunable parameters.
List of model run results.
# TODO
# TODO
Provides convenience function to plot simulation trajectory of result from one run in 2-D phase space.
plot_phase(df_result, x, y, time = NULL)
plot_phase(df_result, x, y, time = NULL)
df_result |
Data frame containing results to plot |
x |
Column to plot as X-axis |
y |
Column to plot as Y-axis |
time |
Column to color the trajectory |
ggplot object of phase portrait
# TODO
# TODO
Provides convenience function to plot simulation trajectory of result from one run as time series. Supports plotting multiple states in faceted plot.
plot_timeSeries(df_result, start_time, end_time, sample_time, tick_time, ...)
plot_timeSeries(df_result, start_time, end_time, sample_time, tick_time, ...)
df_result |
Data frame containing results to plot. Must have column $time. |
start_time |
Plot start time. |
end_time |
Plot end time. |
sample_time |
Time interval to subset data. |
tick_time |
X-axis label break interval time (default 2*sample_time). |
... |
... < |
Support the following features:
flexible start-end time of the series by start_time
and end_time
.
row subset of data by a fixed sample_time
(useful when time step is small for numerical precision).
flexible time tick label spacing by tick_time
.
states to plot are selected by tidy-select
.
Please note that this function does not attempt to do exhaustive sanity check of parameters Make sure yourself that parameters make sense (e.g., end-start time must be longer than sample_time and tick_time, etc.)
ggplot object of faceted time series
# TODO
# TODO
Perform test runs of an Odin model with the specified run parameters and return measured run time for planning large-scale simulation repeats and/or model parameter scans.
run_eta(odin_model, ...)
run_eta(odin_model, ...)
odin_model |
An Odin model |
... |
Passed to |
Estimated run time and resource measured by test run.
# TODO
# TODO