From fa0144fe0dcbe1582bc2a9c1c947b6c5916751b2 Mon Sep 17 00:00:00 2001 From: sichao Date: Wed, 23 Aug 2023 15:46:45 -0400 Subject: [PATCH 01/36] add typing for twostep --- dynamo/estimation/tsc/twostep.py | 41 +++++++++++++++++++++++++++----- 1 file changed, 35 insertions(+), 6 deletions(-) diff --git a/dynamo/estimation/tsc/twostep.py b/dynamo/estimation/tsc/twostep.py index dd8beca7d..3513d7a9d 100644 --- a/dynamo/estimation/tsc/twostep.py +++ b/dynamo/estimation/tsc/twostep.py @@ -1,12 +1,21 @@ +from typing import List, Optional, Tuple, Union + import numpy as np -from scipy.sparse import issparse +from scipy.sparse import csc_matrix, csr_matrix, issparse from tqdm import tqdm from ...tools.utils import calc_norm_loglikelihood, calc_R2, elem_prod, find_extreme from ..csc.utils_velocity import fit_linreg, fit_stochastic_linreg -def fit_slope_stochastic(S, U, US, S2, perc_left=None, perc_right=5): +def fit_slope_stochastic( + S: Union[csc_matrix, csr_matrix, np.ndarray], + U: Union[csc_matrix, csr_matrix, np.ndarray], + US: Union[csc_matrix, csr_matrix, np.ndarray], + S2: Union[csc_matrix, csr_matrix, np.ndarray], + perc_left: Optional[int] = None, + perc_right: int = 5, +) -> Tuple: n_var = S.shape[0] k, all_r2, all_logLL = np.zeros(n_var), np.zeros(n_var), np.zeros(n_var) @@ -28,7 +37,14 @@ def fit_slope_stochastic(S, U, US, S2, perc_left=None, perc_right=5): return k, 0, all_r2, all_logLL -def fit_labeling_synthesis(new, total, t, intercept=False, perc_left=None, perc_right=None): +def fit_labeling_synthesis( + new: Union[csc_matrix, csr_matrix, np.ndarray], + total: Union[csc_matrix, csr_matrix, np.ndarray], + t: Union[List, csc_matrix, csr_matrix, np.ndarray], + intercept: bool = False, + perc_left: Optional[int] = None, + perc_right: Optional[int] = None, +): T = np.unique(t) K = np.zeros(len(T)) R2 = np.zeros(len(T)) @@ -40,18 +56,31 @@ def fit_labeling_synthesis(new, total, t, intercept=False, perc_left=None, perc_ return K, R2 -def compute_gamma_synthesis(K, T): +def compute_gamma_synthesis( + K: Union[csc_matrix, csr_matrix, np.ndarray], + T: Union[csc_matrix, csr_matrix, np.ndarray], +) -> Tuple: gamma, _, r2, _ = fit_linreg(T, -np.log(1 - K)) return gamma, r2 -def compute_velocity_synthesis(N, R, gamma, t): +def compute_velocity_synthesis( + N: Union[csc_matrix, csr_matrix, np.ndarray], + R: Union[csc_matrix, csr_matrix, np.ndarray], + gamma: Union[csc_matrix, csr_matrix, np.ndarray], + t: Union[List, csc_matrix, csr_matrix, np.ndarray], +) -> Union[csc_matrix, csr_matrix, np.ndarray]: k = 1 - np.exp(-np.einsum("i,j->ij", t, gamma)) V = elem_prod(gamma, N) / k - elem_prod(gamma, R) return V -def lin_reg_gamma_synthesis(R, N, time, perc_right=100): +def lin_reg_gamma_synthesis( + R: Union[csc_matrix, csr_matrix, np.ndarray], + N: Union[csc_matrix, csr_matrix, np.ndarray], + time: Union[List, csc_matrix, csr_matrix, np.ndarray], + perc_right: int = 100, +) -> Tuple: n_var = R.shape[0] mean_R2, gamma, r2 = np.zeros(n_var), np.zeros(n_var), np.zeros(n_var) K_list, K_fit_list = [None] * n_var, [None] * n_var From 625755feec4efb98897b15235c8b69fb0bf48cc8 Mon Sep 17 00:00:00 2001 From: sichao Date: Wed, 23 Aug 2023 16:02:14 -0400 Subject: [PATCH 02/36] add docstr for twostep gamma fitting --- dynamo/estimation/tsc/twostep.py | 13 +++++++++++++ 1 file changed, 13 insertions(+) diff --git a/dynamo/estimation/tsc/twostep.py b/dynamo/estimation/tsc/twostep.py index 3513d7a9d..fd1512f37 100644 --- a/dynamo/estimation/tsc/twostep.py +++ b/dynamo/estimation/tsc/twostep.py @@ -81,6 +81,19 @@ def lin_reg_gamma_synthesis( time: Union[List, csc_matrix, csr_matrix, np.ndarray], perc_right: int = 100, ) -> Tuple: + """Under the steady-state assumption, alpha / gamma equals the total RNA. Gamma can be calculated from the slope of + total RNA and new RNA. The relationship can be expressed as: + l(t) = (1 - exp(- gamma * t)) alpha / gamma + + Args: + R: a matrix representing total RNA. + N: a matrix representing new RNA. + time: a matrix with time information. + perc_right: the percentile limitation to find extreme data points. + + Returns: + Gamma, R squared, the slope k, the mean of R squared and the fitted k by time and gamma. + """ n_var = R.shape[0] mean_R2, gamma, r2 = np.zeros(n_var), np.zeros(n_var), np.zeros(n_var) K_list, K_fit_list = [None] * n_var, [None] * n_var From 0354d34a05070494c5e8003effed70e87d225c90 Mon Sep 17 00:00:00 2001 From: sichao Date: Wed, 23 Aug 2023 16:48:10 -0400 Subject: [PATCH 03/36] add docstr for twostep helper funcs --- dynamo/estimation/tsc/twostep.py | 47 ++++++++++++++++++++++++++++++++ 1 file changed, 47 insertions(+) diff --git a/dynamo/estimation/tsc/twostep.py b/dynamo/estimation/tsc/twostep.py index fd1512f37..68aa0182c 100644 --- a/dynamo/estimation/tsc/twostep.py +++ b/dynamo/estimation/tsc/twostep.py @@ -16,6 +16,20 @@ def fit_slope_stochastic( perc_left: Optional[int] = None, perc_right: int = 5, ) -> Tuple: + """Estimate the slope of unspliced and spliced RNA with the GMM method: [u, 2*us + u] = gamma * [s, 2*ss - s]. The + answer is denoted as gamma_k most of the time, which equals gamma/beta under steady state. + + Args: + S: a matrix of the first moments of the spliced RNA. + U: a matrix of the first moments of the unspliced RNA. + US: a matrix of the cross moments of unspliced/spliced RNA. + S2: a matrix of the second moments of spliced RNA. + perc_left: the left percentile limitation to find extreme data points. + perc_right: the right percentile limitation to find extreme data points. + + Returns: + The slope, intercept, R squared and log likelihood. + """ n_var = S.shape[0] k, all_r2, all_logLL = np.zeros(n_var), np.zeros(n_var), np.zeros(n_var) @@ -45,6 +59,19 @@ def fit_labeling_synthesis( perc_left: Optional[int] = None, perc_right: Optional[int] = None, ): + """Calculate the slope of total and new RNA under steady-state assumption. + + Args: + new: a matrix representing new RNA. Can be expression or the first moments. + total: a matrix representing total RNA. Can be expression or the first moments. + t: a matrix of time information. + intercept: whether to perform the linear regression with intercept. + perc_left: the left percentile limitation to find extreme data points. + perc_right: the right percentile limitation to find extreme data points. + + Returns: + The slope K and R squared of linear regression. + """ T = np.unique(t) K = np.zeros(len(T)) R2 = np.zeros(len(T)) @@ -60,6 +87,15 @@ def compute_gamma_synthesis( K: Union[csc_matrix, csr_matrix, np.ndarray], T: Union[csc_matrix, csr_matrix, np.ndarray], ) -> Tuple: + """Calculate gamma as the linear regression results of given time and log(1 - slope k). + + Args: + K: a matrix of the slope k. + T: a matrix of time information. + + Returns: + The gamma and R squared of linear regression. + """ gamma, _, r2, _ = fit_linreg(T, -np.log(1 - K)) return gamma, r2 @@ -70,6 +106,17 @@ def compute_velocity_synthesis( gamma: Union[csc_matrix, csr_matrix, np.ndarray], t: Union[List, csc_matrix, csr_matrix, np.ndarray], ) -> Union[csc_matrix, csr_matrix, np.ndarray]: + """Calculate the velocity of total RNA with a physical time unit: velocity = (gamma / k) N - gamma * R. + + Args: + N: a matrix representing new RNA. + R: a matrix representing total RNA. + gamma: a matrix of degradation rate. + t: a matrix of time information. + + Returns: + The velocity. + """ k = 1 - np.exp(-np.einsum("i,j->ij", t, gamma)) V = elem_prod(gamma, N) / k - elem_prod(gamma, R) return V From 421ddc79625d281f3a1490c8def389d732837f1e Mon Sep 17 00:00:00 2001 From: sichao Date: Wed, 23 Aug 2023 16:49:08 -0400 Subject: [PATCH 04/36] update twostelin_reg_gamma_synthesis docstr --- dynamo/estimation/tsc/twostep.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/dynamo/estimation/tsc/twostep.py b/dynamo/estimation/tsc/twostep.py index 68aa0182c..426ee0b29 100644 --- a/dynamo/estimation/tsc/twostep.py +++ b/dynamo/estimation/tsc/twostep.py @@ -133,8 +133,8 @@ def lin_reg_gamma_synthesis( l(t) = (1 - exp(- gamma * t)) alpha / gamma Args: - R: a matrix representing total RNA. - N: a matrix representing new RNA. + R: a matrix representing total RNA. Can be expression or the first moments. + N: a matrix representing new RNA. Can be expression or the first moments. time: a matrix with time information. perc_right: the percentile limitation to find extreme data points. From 0f2236f7d9e7b2bdadd213a01e6dd56087b62bec Mon Sep 17 00:00:00 2001 From: sichao Date: Thu, 24 Aug 2023 14:45:34 -0400 Subject: [PATCH 05/36] add typing for estimation_kinetic funcs --- dynamo/estimation/tsc/estimation_kinetic.py | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/dynamo/estimation/tsc/estimation_kinetic.py b/dynamo/estimation/tsc/estimation_kinetic.py index 7034bf652..e900cb3b5 100755 --- a/dynamo/estimation/tsc/estimation_kinetic.py +++ b/dynamo/estimation/tsc/estimation_kinetic.py @@ -1,5 +1,7 @@ +from typing import List, Optional, Tuple, Union import warnings +import numpy as np from scipy.optimize import least_squares from scipy.stats import chi2 @@ -9,26 +11,26 @@ from .utils_kinetic import * -def guestimate_alpha(x_data, time): +def guestimate_alpha(x_data: np.ndarray, time: np.ndarray) -> np.ndarray: """Roughly estimate p0 for kinetics data.""" imax = np.argmax(x_data) alpha = x_data[imax] / time[imax] return alpha -def guestimate_gamma(x_data, time): +def guestimate_gamma(x_data: np.ndarray, time: np.ndarray) -> np.ndarray: """Roughly estimate gamma0 with the assumption that time starts at 0 for degradation data.""" ga0 = np.clip(np.log(max(x_data[0], 0) / (x_data[-1] + 1e-6)) / time[-1], 1e-3, 1e3) return ga0 -def guestimate_init_cond(x_data): +def guestimate_init_cond(x_data: np.ndarray) -> np.ndarray: """Roughly estimate x0 for degradation data.""" x0 = np.clip(np.max(x_data, 1), 1e-4, np.inf) return x0 -def guestimate_p0_kinetic_chase(x_data, time): +def guestimate_p0_kinetic_chase(x_data: np.ndarray, time: np.ndarray) -> Tuple: t0 = np.min(time) x0 = np.mean(x_data[time == t0]) idx = time != t0 From 33b1867483e153df677cab3dbf54ede49f98e32a Mon Sep 17 00:00:00 2001 From: sichao Date: Thu, 24 Aug 2023 15:27:16 -0400 Subject: [PATCH 06/36] add docstr for estimation_kinetic func --- dynamo/estimation/tsc/estimation_kinetic.py | 39 +++++++++++++++++++-- 1 file changed, 36 insertions(+), 3 deletions(-) diff --git a/dynamo/estimation/tsc/estimation_kinetic.py b/dynamo/estimation/tsc/estimation_kinetic.py index e900cb3b5..0d1df4a75 100755 --- a/dynamo/estimation/tsc/estimation_kinetic.py +++ b/dynamo/estimation/tsc/estimation_kinetic.py @@ -12,25 +12,58 @@ def guestimate_alpha(x_data: np.ndarray, time: np.ndarray) -> np.ndarray: - """Roughly estimate p0 for kinetics data.""" + """Roughly estimate alpha for kinetics data, which is simply the averaged ratio of new RNA and labeling time. + + Args: + x_data: a matrix representing RNA data. + time: a matrix of labeling time. + + Returns: + The estimated alpha. + """ imax = np.argmax(x_data) alpha = x_data[imax] / time[imax] return alpha def guestimate_gamma(x_data: np.ndarray, time: np.ndarray) -> np.ndarray: - """Roughly estimate gamma0 with the assumption that time starts at 0 for degradation data.""" + """Roughly estimate gamma0 with the assumption that time starts at 0. + + Args: + x_data: a matrix representing RNA data. + time: a matrix of labeling time. + + Returns: + The estimated gamma. + """ ga0 = np.clip(np.log(max(x_data[0], 0) / (x_data[-1] + 1e-6)) / time[-1], 1e-3, 1e3) return ga0 def guestimate_init_cond(x_data: np.ndarray) -> np.ndarray: - """Roughly estimate x0 for degradation data.""" + """Roughly estimate x0 for degradation data. + + Args: + x_data: a matrix representing RNA data. + + Returns: + The estimated x0. + """ x0 = np.clip(np.max(x_data, 1), 1e-4, np.inf) return x0 def guestimate_p0_kinetic_chase(x_data: np.ndarray, time: np.ndarray) -> Tuple: + """Roughly estimated alpha, gamma and initial conditions for the kinetic chase experiment. Initail conditions are + the average abundance of labeled RNAs across all cells belonging to the initial labeling time point. + + Args: + x_data: a matrix representing RNA data. + time: a matrix of labeling time. + + Returns: + The estimated alpha0, gamma0 and x0. + """ t0 = np.min(time) x0 = np.mean(x_data[time == t0]) idx = time != t0 From c00ec7db86e0906f38224ecba948d163071daddc Mon Sep 17 00:00:00 2001 From: sichao Date: Thu, 24 Aug 2023 16:05:08 -0400 Subject: [PATCH 07/36] add typing for class kinetic_estimation --- dynamo/estimation/tsc/estimation_kinetic.py | 64 +++++++++++++-------- 1 file changed, 39 insertions(+), 25 deletions(-) diff --git a/dynamo/estimation/tsc/estimation_kinetic.py b/dynamo/estimation/tsc/estimation_kinetic.py index 0d1df4a75..ce90e840a 100755 --- a/dynamo/estimation/tsc/estimation_kinetic.py +++ b/dynamo/estimation/tsc/estimation_kinetic.py @@ -91,7 +91,7 @@ class kinetic_estimation: as well as two functions: integrate (numerical integration), solve (analytical method). """ - def __init__(self, param_ranges, x0_ranges, simulator): + def __init__(self, param_ranges: np.ndarray, x0_ranges: np.ndarray, simulator: LinearODE) -> None: self.simulator = simulator self.ranges = [] @@ -114,7 +114,7 @@ def __init__(self, param_ranges, x0_ranges, simulator): self.popt = None self.cost = None - def sample_p0(self, samples=1, method="lhs"): + def sample_p0(self, samples: int = 1, method: str = "lhs") -> np.ndarray: ret = np.zeros((samples, self.n_params)) if method == "lhs": ret = self._lhsclassic(samples) @@ -127,7 +127,7 @@ def sample_p0(self, samples=1, method="lhs"): ret[n, i] = r * (self.ranges[i][1] - self.ranges[i][0]) + self.ranges[i][0] return ret - def _lhsclassic(self, samples): + def _lhsclassic(self, samples: int) -> np.ndarray: # From PyDOE # Generate the intervals # from .utils import lhsclassic @@ -135,48 +135,55 @@ def _lhsclassic(self, samples): return H - def get_bound(self, axis): + def get_bound(self, axis: int) -> np.ndarray: ret = np.zeros(self.n_params) for i in range(self.n_params): ret[i] = self.ranges[i][axis] return ret - def normalize_data(self, X): + def normalize_data(self, X: np.ndarray) -> np.ndarray: return np.log1p(X) - def extract_data_from_simulator(self, t=None, **kwargs): + def extract_data_from_simulator(self, t: Optional[np.ndarray] = None, **kwargs) -> np.ndarray: if t is None: return self.simulator.x.T else: x = self.simulator.integrate(t=t, **kwargs) return x.T - def assemble_kin_params(self, unfixed_params): + def assemble_kin_params(self, unfixed_params: np.ndarray) -> np.ndarray: p = np.array(self.fixed_parameters[: self.n_tot_kin_params], copy=True) p[np.isnan(p)] = unfixed_params[: self.n_kin_params] return p - def assemble_x0(self, unfixed_params): + def assemble_x0(self, unfixed_params: np.ndarray) -> np.ndarray: p = np.array(self.fixed_parameters[self.n_tot_kin_params :], copy=True) p[np.isnan(p)] = unfixed_params[self.n_kin_params :] return p - def set_params(self, params): + def set_params(self, params: np.ndarray) -> None: self.simulator.set_params(*self.assemble_kin_params(params)) - def get_opt_kin_params(self): + def get_opt_kin_params(self) -> Optional[np.ndarray]: if self.popt is not None: return self.assemble_kin_params(self.popt) else: return None - def get_opt_x0_params(self): + def get_opt_x0_params(self) -> Optional[np.ndarray]: if self.popt is not None: return self.assemble_x0(self.popt) else: return None - def f_lsq(self, params, t, x_data, method=None, normalize=True): + def f_lsq( + self, + params: np.ndarray, + t: np.ndarray, + x_data: np.ndarray, + method: Optional[str] = None, + normalize: bool = True, + ) -> np.ndarray: self.set_params(params) x0 = self.assemble_x0(params) self.simulator.integrate(t, x0, method) @@ -187,15 +194,15 @@ def f_lsq(self, params, t, x_data, method=None, normalize=True): def fit_lsq( self, - t, - x_data, - p0=None, - n_p0=1, - bounds=None, - sample_method="lhs", - method=None, - normalize=True, - ): + t: np.ndarray, + x_data: np.ndarray, + p0: Optional[np.ndarray] = None, + n_p0: int = 1, + bounds: Optional[Tuple] = None, + sample_method: str = "lhs", + method: Optional[str] = None, + normalize: bool = True, + ) -> Tuple: """Fit time-seris data using least squares Arguments @@ -255,19 +262,26 @@ def fit_lsq( self.cost = costs[i_min] return self.popt, self.cost - def export_parameters(self): + def export_parameters(self) -> Optional[np.ndarray]: return self.get_opt_kin_params() - def export_model(self, reinstantiate=True): + def export_model(self, reinstantiate: bool = True) -> LinearODE: if reinstantiate: return self.simulator.__class__() else: return self.simulator - def get_SSE(self): + def get_SSE(self) -> float: return self.cost - def test_chi2(self, t, x_data, species=None, method="matrix", normalize=True): + def test_chi2( + self, + t: np.ndarray, + x_data: np.ndarray, + species: Optional[Union[List, np.ndarray]] = None, + method: str = "matrix", + normalize: bool = True, + ) -> Tuple: """perform a Pearson's chi-square test. The statistics is computed as: sum_i (O_i - E_i)^2 / E_i, where O_i is the data and E_i is the model predication. The data can be either 1. stratified moments: 't' is an array of k distinct time points, 'x_data' is a m-by-k matrix of data, where m is the number of species. From ddea05b975a797456c4158af8194b100b9d3aeb6 Mon Sep 17 00:00:00 2001 From: sichao Date: Thu, 24 Aug 2023 16:55:38 -0400 Subject: [PATCH 08/36] add typing to linearODE and MixtureModels --- dynamo/estimation/tsc/utils_kinetic.py | 28 ++++++++++++++------------ 1 file changed, 15 insertions(+), 13 deletions(-) diff --git a/dynamo/estimation/tsc/utils_kinetic.py b/dynamo/estimation/tsc/utils_kinetic.py index 8f3d2be9a..f7a76e272 100755 --- a/dynamo/estimation/tsc/utils_kinetic.py +++ b/dynamo/estimation/tsc/utils_kinetic.py @@ -1,3 +1,5 @@ +from typing import List, Optional, Tuple + import numpy as np from scipy.integrate import odeint @@ -6,7 +8,7 @@ class LinearODE: - def __init__(self, n_species, x0=None): + def __init__(self, n_species: int, x0: Optional[np.ndarray] = None) -> None: """A general class for linear odes""" self.n_species = n_species # solution @@ -19,12 +21,12 @@ def __init__(self, n_species, x0=None): self.methods = ["numerical", "matrix"] self.default_method = "matrix" - def ode_func(self, x, t): + def ode_func(self, x: np.ndarray, t: np.ndarray) -> np.ndarray: """Implement your own ODE functions here such that dx=f(x, t)""" dx = np.zeros(len(x)) return dx - def integrate(self, t, x0=None, method=None): + def integrate(self, t: np.ndarray, x0: Optional[np.ndarray] = None, method: Optional[str] = None) -> np.ndarray: method = self.default_method if method is None else method if method == "matrix": sol = self.integrate_matrix(t, x0) @@ -36,7 +38,7 @@ def integrate(self, t, x0=None, method=None): self.t = t return sol - def integrate_numerical(self, t, x0=None): + def integrate_numerical(self, t: np.ndarray, x0: Optional[np.ndarray] = None) -> np.ndarray: if x0 is None: x0 = self.x0 else: @@ -44,20 +46,20 @@ def integrate_numerical(self, t, x0=None): sol = odeint(self.ode_func, x0, t) return sol - def reset(self): + def reset(self) -> None: # reset solutions self.t = None self.x = None self.K = None self.p = None - def computeKnp(self): + def computeKnp(self) -> Tuple[np.ndarray, np.ndarray]: """Implement your own vectorized ODE functions here such that dx = Kx + p""" K = np.zeros((self.n_species, self.n_species)) p = np.zeros(self.n_species) return K, p - def integrate_matrix(self, t, x0=None): + def integrate_matrix(self, t: np.ndarray, x0: Optional[np.ndarray] = None) -> np.ndarray: # t0 = t[0] t0 = 0 if x0 is None: @@ -88,7 +90,7 @@ def integrate_matrix(self, t, x0=None): class MixtureModels: - def __init__(self, models, param_distributor): + def __init__(self, models: LinearODE, param_distributor: List) -> None: """A general class for linear odes""" self.n_models = len(models) self.models = models @@ -101,7 +103,7 @@ def __init__(self, models, param_distributor): self.methods = ["numerical", "matrix"] self.default_method = "matrix" - def integrate(self, t, x0=None, method=None): + def integrate(self, t: np.ndarray, x0: Optional[np.ndarray] = None, method: Optional[str] = None) -> None: self.x = np.zeros((len(t), np.sum(self.n_species))) for i, mdl in enumerate(self.models): x0_ = None if x0 is None else x0[self.get_model_species(i)] @@ -110,22 +112,22 @@ def integrate(self, t, x0=None, method=None): self.x[:, self.get_model_species(i)] = mdl.x self.t = np.array(self.models[0].t, copy=True) - def get_model_species(self, model_index): + def get_model_species(self, model_index: int) -> int: id = np.hstack((0, np.cumsum(self.n_species))) idx = np.arange(id[-1] + 1) return idx[id[model_index] : id[model_index + 1]] - def reset(self): + def reset(self) -> None: # reset solutions self.t = None self.x = None for mdl in self.models: mdl.reset() - def param_mixer(self, *params): + def param_mixer(self, *params: Tuple) -> Tuple: return params - def set_params(self, *params): + def set_params(self, *params: Tuple) -> None: params = self.param_mixer(*params) for i, mdl in enumerate(self.models): idx = self.distributor[i] From 35e490d21d61432ebdfdeb54d99bbbcaa9134630 Mon Sep 17 00:00:00 2001 From: sichao Date: Thu, 24 Aug 2023 18:12:21 -0400 Subject: [PATCH 09/36] add docstring for linearODE and MixtureModels --- dynamo/estimation/tsc/utils_kinetic.py | 86 ++++++++++++++++++++++++-- 1 file changed, 82 insertions(+), 4 deletions(-) diff --git a/dynamo/estimation/tsc/utils_kinetic.py b/dynamo/estimation/tsc/utils_kinetic.py index f7a76e272..145d185c7 100755 --- a/dynamo/estimation/tsc/utils_kinetic.py +++ b/dynamo/estimation/tsc/utils_kinetic.py @@ -8,8 +8,14 @@ class LinearODE: + """A general class for linear odes.""" def __init__(self, n_species: int, x0: Optional[np.ndarray] = None) -> None: - """A general class for linear odes""" + """Initialize the LinearODE object. + + Args: + n_species: the number of species. + x0: the initial condition of variable x. + """ self.n_species = n_species # solution self.t = None @@ -22,11 +28,29 @@ def __init__(self, n_species: int, x0: Optional[np.ndarray] = None) -> None: self.default_method = "matrix" def ode_func(self, x: np.ndarray, t: np.ndarray) -> np.ndarray: - """Implement your own ODE functions here such that dx=f(x, t)""" + """ODE functions to be implemented in the derived class such that dx=f(x, t). + + Args: + x: the variable. + t: the aray of time. + + Returns: + The derivatives dx. + """ dx = np.zeros(len(x)) return dx def integrate(self, t: np.ndarray, x0: Optional[np.ndarray] = None, method: Optional[str] = None) -> np.ndarray: + """Integrate the ODE using the given time values. + + Args: + t: array of time values. + x0: array of initial conditions. + method: the method to integrate, including "matrix" and "numerical". + + Returns: + Array containing the integrated solution over the specified time values. + """ method = self.default_method if method is None else method if method == "matrix": sol = self.integrate_matrix(t, x0) @@ -39,6 +63,15 @@ def integrate(self, t: np.ndarray, x0: Optional[np.ndarray] = None, method: Opti return sol def integrate_numerical(self, t: np.ndarray, x0: Optional[np.ndarray] = None) -> np.ndarray: + """Numerically integrate the ODE using the given time values. + + Args: + t: array of time values. + x0: array of initial conditions. + + Returns: + Array containing the integrated solution over the specified time values. + """ if x0 is None: x0 = self.x0 else: @@ -47,6 +80,7 @@ def integrate_numerical(self, t: np.ndarray, x0: Optional[np.ndarray] = None) -> return sol def reset(self) -> None: + """Reset the ODE to initial state.""" # reset solutions self.t = None self.x = None @@ -54,12 +88,21 @@ def reset(self) -> None: self.p = None def computeKnp(self) -> Tuple[np.ndarray, np.ndarray]: - """Implement your own vectorized ODE functions here such that dx = Kx + p""" + """The vectorized ODE functions to be implemented in the derived class such that dx = Kx + p.""" K = np.zeros((self.n_species, self.n_species)) p = np.zeros(self.n_species) return K, p def integrate_matrix(self, t: np.ndarray, x0: Optional[np.ndarray] = None) -> np.ndarray: + """Integrate the system of ordinary differential equations (ODEs) using matrix exponential. + + Args: + t: array of time values. + x0: array of initial conditions. + + Returns: + Array containing the integrated solution over the specified time values. + """ # t0 = t[0] t0 = 0 if x0 is None: @@ -90,8 +133,14 @@ def integrate_matrix(self, t: np.ndarray, x0: Optional[np.ndarray] = None) -> np class MixtureModels: + """The base class for mixture models.""" def __init__(self, models: LinearODE, param_distributor: List) -> None: - """A general class for linear odes""" + """Initialize the MixtureModels class. + + Args: + models: the models to mix. + param_distributor: the index to assign parameters. + """ self.n_models = len(models) self.models = models self.n_species = np.array([mdl.n_species for mdl in self.models]) @@ -104,6 +153,13 @@ def __init__(self, models: LinearODE, param_distributor: List) -> None: self.default_method = "matrix" def integrate(self, t: np.ndarray, x0: Optional[np.ndarray] = None, method: Optional[str] = None) -> None: + """Integrate with time values for all models. + + Args: + t: array of time values. + x0: array of initial conditions. + method: the method or methods to integrate, including "matrix" and "numerical". + """ self.x = np.zeros((len(t), np.sum(self.n_species))) for i, mdl in enumerate(self.models): x0_ = None if x0 is None else x0[self.get_model_species(i)] @@ -113,11 +169,20 @@ def integrate(self, t: np.ndarray, x0: Optional[np.ndarray] = None, method: Opti self.t = np.array(self.models[0].t, copy=True) def get_model_species(self, model_index: int) -> int: + """Get the indices of species associated with the specified model. + + Args: + model_index: index of the model. + + Returns: + Array containing the indices of species associated with the specified model. + """ id = np.hstack((0, np.cumsum(self.n_species))) idx = np.arange(id[-1] + 1) return idx[id[model_index] : id[model_index + 1]] def reset(self) -> None: + """Reset all models.""" # reset solutions self.t = None self.x = None @@ -125,9 +190,22 @@ def reset(self) -> None: mdl.reset() def param_mixer(self, *params: Tuple) -> Tuple: + """Unpack the given parameters. + + Args: + params: tuple of parameters. + + Returns: + The unpacked tuple. + """ return params def set_params(self, *params: Tuple) -> None: + """Set parameters for all models. + + Args: + params: tuple of parameters. + """ params = self.param_mixer(*params) for i, mdl in enumerate(self.models): idx = self.distributor[i] From 0857310d36396219656d851f011f0ed06648389d Mon Sep 17 00:00:00 2001 From: sichao Date: Thu, 24 Aug 2023 18:20:38 -0400 Subject: [PATCH 10/36] debug typing in integrate --- dynamo/estimation/tsc/utils_kinetic.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/dynamo/estimation/tsc/utils_kinetic.py b/dynamo/estimation/tsc/utils_kinetic.py index 145d185c7..aedb21ae2 100755 --- a/dynamo/estimation/tsc/utils_kinetic.py +++ b/dynamo/estimation/tsc/utils_kinetic.py @@ -1,4 +1,4 @@ -from typing import List, Optional, Tuple +from typing import List, Optional, Tuple, Union import numpy as np from scipy.integrate import odeint @@ -152,7 +152,7 @@ def __init__(self, models: LinearODE, param_distributor: List) -> None: self.methods = ["numerical", "matrix"] self.default_method = "matrix" - def integrate(self, t: np.ndarray, x0: Optional[np.ndarray] = None, method: Optional[str] = None) -> None: + def integrate(self, t: np.ndarray, x0: Optional[np.ndarray] = None, method: Optional[Union[str, List]] = None) -> None: """Integrate with time values for all models. Args: From 0e9a6bffce55f0bdeb03a39d85b5df0c55274f9f Mon Sep 17 00:00:00 2001 From: sichao Date: Thu, 24 Aug 2023 18:53:58 -0400 Subject: [PATCH 11/36] add docstr for kinetic_estimation --- dynamo/estimation/tsc/estimation_kinetic.py | 200 ++++++++++++++------ 1 file changed, 142 insertions(+), 58 deletions(-) diff --git a/dynamo/estimation/tsc/estimation_kinetic.py b/dynamo/estimation/tsc/estimation_kinetic.py index ce90e840a..22afa72dd 100755 --- a/dynamo/estimation/tsc/estimation_kinetic.py +++ b/dynamo/estimation/tsc/estimation_kinetic.py @@ -75,23 +75,20 @@ def guestimate_p0_kinetic_chase(x_data: np.ndarray, time: np.ndarray) -> Tuple: class kinetic_estimation: - """A general parameter estimation framework for all types of time-seris data - - Arguments - --------- - param_ranges: :class:`~numpy.ndarray` - A n-by-2 numpy array containing the lower and upper ranges of n parameters - (and initial conditions if not fixed). - x0_ranges: :class:`~numpy.ndarray` - Lower and upper bounds for initial conditions for the integrators. - To fix a parameter, set its lower and upper bounds to the same value. - simulator: :class:`utils_kinetic.Linear_ODE` - An instance of python class which solves ODEs. It should have properties 't' (k time points, 1d numpy array), - 'x0' (initial conditions for m species, 1d numpy array), and 'x' (solution, k-by-m array), - as well as two functions: integrate (numerical integration), solve (analytical method). - """ + """A general parameter estimation framework for all types of time-seris data.""" def __init__(self, param_ranges: np.ndarray, x0_ranges: np.ndarray, simulator: LinearODE) -> None: + """Initialize the kinetic_estimation class. + + Args: + param_ranges: A n-by-2 numpy array containing the lower and upper ranges of n parameters (and initial + conditions if not fixed). + x0_ranges: Lower and upper bounds for initial conditions for the integrators. To fix a parameter, set its + lower and upper bounds to the same value. + simulator: An instance of python class which solves ODEs. It should have properties 't' (k time points, + 1d numpy array), 'x0' (initial conditions for m species, 1d numpy array), and 'x' (solution, k-by-m + array), as well as two functions: integrate (numerical integration), solve (analytical method). + """ self.simulator = simulator self.ranges = [] @@ -115,6 +112,15 @@ def __init__(self, param_ranges: np.ndarray, x0_ranges: np.ndarray, simulator: L self.cost = None def sample_p0(self, samples: int = 1, method: str = "lhs") -> np.ndarray: + """Sample the initial parameters with either Latin Hypercube Sampling or random method. + + Args: + samples: the number of samples. + method: the sampling method. Only support "lhs" and random sampling. + + Returns: + The sampled array. + """ ret = np.zeros((samples, self.n_params)) if method == "lhs": ret = self._lhsclassic(samples) @@ -128,6 +134,14 @@ def sample_p0(self, samples: int = 1, method: str = "lhs") -> np.ndarray: return ret def _lhsclassic(self, samples: int) -> np.ndarray: + """Run the Latin Hypercube Sampling function. + + Args: + samples: the number of samples. + + Returns: + The sampled data array. + """ # From PyDOE # Generate the intervals # from .utils import lhsclassic @@ -136,15 +150,40 @@ def _lhsclassic(self, samples: int) -> np.ndarray: return H def get_bound(self, axis: int) -> np.ndarray: + """Get the bounds of the specified axis for all parameters. + + Args: + axis: the index of axis. + + Returns: + An array containing the bounds of the specified axis for all parameters. + """ ret = np.zeros(self.n_params) for i in range(self.n_params): ret[i] = self.ranges[i][axis] return ret def normalize_data(self, X: np.ndarray) -> np.ndarray: + """Perform log1p normalization on the data. + + Args: + X: target data to normalize. + + Returns: + The normalized data. + """ return np.log1p(X) def extract_data_from_simulator(self, t: Optional[np.ndarray] = None, **kwargs) -> np.ndarray: + """Extract data from the ODE simulator. + + Args: + t: the time information. If provided, the data will be integrated with time information. + kwargs: additional keyword arguments. + + Returns: + The variable from ODE simulator. + """ if t is None: return self.simulator.x.T else: @@ -152,25 +191,56 @@ def extract_data_from_simulator(self, t: Optional[np.ndarray] = None, **kwargs) return x.T def assemble_kin_params(self, unfixed_params: np.ndarray) -> np.ndarray: + """Assemble the kinetic parameters array. + + Args: + unfixed_params: array of unfixed parameters. + + Returns: + The assembled kinetic parameters. + """ p = np.array(self.fixed_parameters[: self.n_tot_kin_params], copy=True) p[np.isnan(p)] = unfixed_params[: self.n_kin_params] return p def assemble_x0(self, unfixed_params: np.ndarray) -> np.ndarray: + """Assemble the initial conditions array. + + Args: + unfixed_params: array of unfixed parameters. + + Returns: + The assembled initial conditions. + """ p = np.array(self.fixed_parameters[self.n_tot_kin_params :], copy=True) p[np.isnan(p)] = unfixed_params[self.n_kin_params :] return p def set_params(self, params: np.ndarray) -> None: + """Set the parameters of the simulator using assembled kinetic parameters. + + Args: + params: array of assembled kinetic parameters. + """ self.simulator.set_params(*self.assemble_kin_params(params)) def get_opt_kin_params(self) -> Optional[np.ndarray]: + """Get the optimized kinetic parameters. + + Returns: + Array containing the optimized kinetic parameters, or None if not available. + """ if self.popt is not None: return self.assemble_kin_params(self.popt) else: return None def get_opt_x0_params(self) -> Optional[np.ndarray]: + """Get the optimized initial conditions. + + Returns: + Array containing the optimized initial conditions, or None if not available. + """ if self.popt is not None: return self.assemble_x0(self.popt) else: @@ -184,6 +254,18 @@ def f_lsq( method: Optional[str] = None, normalize: bool = True, ) -> np.ndarray: + """Calculate the difference between simulated and observed data for least squares fitting. + + Args: + params: array of parameters for the simulation. + t: array of time values. + x_data: the input array. + method: method for integration. + normalize: whether to normalize data. + + Returns: + Residuals representing the differences between simulated and observed data (flattened). + """ self.set_params(params) x0 = self.assemble_x0(params) self.simulator.integrate(t, x0, method) @@ -203,37 +285,26 @@ def fit_lsq( method: Optional[str] = None, normalize: bool = True, ) -> Tuple: - """Fit time-seris data using least squares - - Arguments - --------- - t: :class:`~numpy.ndarray` - A numpy array of n time points. - x_data: :class:`~numpy.ndarray` - A m-by-n numpy a array of m species, each having n values for the n time points. - p0: :class:`numpy.ndarray`, optional, default: None - Initial guesses of parameters. If None, a random number is generated within the bounds. - n_p0: int, optional, default: 1 - Number of initial guesses. - bounds: tuple, optional, default: None - Lower and upper bounds for parameters. - sample_method: str, optional, default: `lhs` - Method used for sampling initial guesses of parameters: + """Fit time-seris data using the least squares method. + + This method iteratively optimizes the parameters for different initial conditions (p0) and returns + the optimized parameters and associated cost. + + Args: + t: a numpy array of n time points. + x_data: an m-by-n numpy array of m species, each having n values for the n time points. + p0: initial guesses of parameters. If None, a random number is generated within the bounds. + n_p0: number of initial guesses. + bounds: lower and upper bounds for parameters. + sample_method: method used for sampling initial guesses of parameters: `lhs`: latin hypercube sampling; `uniform`: uniform random sampling. - method: str or None, optional, default: None - Method used for solving ODEs. See options in simulator classes. - If None, default method is used. - normalize: bool, optional, default: True - Whether or not normalize values in x_data across species, so that large values do + method: method used for solving ODEs. See options in simulator classes. + normalize: whether to normalize values in x_data across species, so that large values do not dominate the optimizer. - Returns - --------- - popt: :class:`~numpy.ndarray` - Optimal parameters. - cost: float - The cost function evaluated at the optimum. + Returns: + Optimal parameters and the cost function evaluated at the optimum. """ if p0 is None: p0 = self.sample_p0(n_p0, sample_method) @@ -263,15 +334,33 @@ def fit_lsq( return self.popt, self.cost def export_parameters(self) -> Optional[np.ndarray]: + """Export the optimized kinetic parameters. + + Returns: + Array containing the optimized kinetic parameters, or None if not available. + """ return self.get_opt_kin_params() def export_model(self, reinstantiate: bool = True) -> LinearODE: + """Export the simulator model. + + Args: + reinstantiate: whether to reinstantiate the model class (default: True). + + Returns: + Exported simulator model. + """ if reinstantiate: return self.simulator.__class__() else: return self.simulator def get_SSE(self) -> float: + """Get the sum of squared errors (SSE) from the least squares fitting. + + Returns: + Sum of squared errors (SSE). + """ return self.cost def test_chi2( @@ -282,24 +371,19 @@ def test_chi2( method: str = "matrix", normalize: bool = True, ) -> Tuple: - """perform a Pearson's chi-square test. The statistics is computed as: sum_i (O_i - E_i)^2 / E_i, where O_i is the data and E_i is the model predication. - - The data can be either 1. stratified moments: 't' is an array of k distinct time points, 'x_data' is a m-by-k matrix of data, where m is the number of species. - or 2. raw data: 't' is an array of k time points for k cells, 'x_data' is a m-by-k matrix of data, where m is the number of species. - Note that if the method is 'numerical', t has to monotonically increasing. - + """Perform a Pearson's chi-square test. The statistics is computed as: sum_i (O_i - E_i)^2 / E_i, where O_i is + the data and E_i is the model predication. + + The data can be either: + 1. stratified moments: 't' is an array of k distinct time points, 'x_data' is an m-by-k matrix of data, + where m is the number of species. + Or + 2. raw data: 't' is an array of k time points for k cells, 'x_data' is an m-by-k matrix of data, where m is + the number of species. Note that if the method is 'numerical', t has to monotonically increasing. If not all species are included in the data, use 'species' to specify the species of interest. - Returns - ------- - p: float - The p-value of a one-tailed chi-square test. - - c2: float - The chi-square statistics. - - df: int - Degree of freedom. + Returns: + The p-value of a one-tailed chi-square test, the chi-square statistics and degree of freedom. """ if x_data.ndim == 1: x_data = x_data[None] From c157abe3eecbe4e06146f33eef3e00d0cd531231 Mon Sep 17 00:00:00 2001 From: sichao Date: Thu, 24 Aug 2023 19:03:28 -0400 Subject: [PATCH 12/36] debug tuple typing in estimation_kinetic --- dynamo/estimation/tsc/estimation_kinetic.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/dynamo/estimation/tsc/estimation_kinetic.py b/dynamo/estimation/tsc/estimation_kinetic.py index 22afa72dd..1a21eacdd 100755 --- a/dynamo/estimation/tsc/estimation_kinetic.py +++ b/dynamo/estimation/tsc/estimation_kinetic.py @@ -53,7 +53,7 @@ def guestimate_init_cond(x_data: np.ndarray) -> np.ndarray: return x0 -def guestimate_p0_kinetic_chase(x_data: np.ndarray, time: np.ndarray) -> Tuple: +def guestimate_p0_kinetic_chase(x_data: np.ndarray, time: np.ndarray) -> Tuple[np.ndarray, np.ndarray, np.ndarray]: """Roughly estimated alpha, gamma and initial conditions for the kinetic chase experiment. Initail conditions are the average abundance of labeled RNAs across all cells belonging to the initial labeling time point. @@ -280,11 +280,11 @@ def fit_lsq( x_data: np.ndarray, p0: Optional[np.ndarray] = None, n_p0: int = 1, - bounds: Optional[Tuple] = None, + bounds: Optional[Tuple[Union[float, int], Union[float, int]]] = None, sample_method: str = "lhs", method: Optional[str] = None, normalize: bool = True, - ) -> Tuple: + ) -> Tuple[np.ndarray, float]: """Fit time-seris data using the least squares method. This method iteratively optimizes the parameters for different initial conditions (p0) and returns @@ -370,7 +370,7 @@ def test_chi2( species: Optional[Union[List, np.ndarray]] = None, method: str = "matrix", normalize: bool = True, - ) -> Tuple: + ) -> Tuple[float, float, int]: """Perform a Pearson's chi-square test. The statistics is computed as: sum_i (O_i - E_i)^2 / E_i, where O_i is the data and E_i is the model predication. From 42ea44a2a4a35f5085c5c5fb11eb34091ad75dbd Mon Sep 17 00:00:00 2001 From: sichao Date: Fri, 25 Aug 2023 16:40:40 -0400 Subject: [PATCH 13/36] add docstr for Estimation_Degradation --- dynamo/estimation/tsc/estimation_kinetic.py | 46 +++++++++++++++++++++ 1 file changed, 46 insertions(+) diff --git a/dynamo/estimation/tsc/estimation_kinetic.py b/dynamo/estimation/tsc/estimation_kinetic.py index 1a21eacdd..d46b32d08 100755 --- a/dynamo/estimation/tsc/estimation_kinetic.py +++ b/dynamo/estimation/tsc/estimation_kinetic.py @@ -408,23 +408,69 @@ def test_chi2( class Estimation_Degradation(kinetic_estimation): + """The base parameters, estimation class for degradation experiments.""" def __init__(self, ranges, x0, simulator): + """Initialize the Estimation_Degradation object. + + Args: + ranges: the lower and upper ranges of parameters. + x0: initial conditions. + simulator: instance of the Python class to solve ODEs. + """ self.kin_param_keys = np.array(["alpha", "gamma"]) super().__init__(np.vstack((np.zeros(2), ranges)), x0, simulator) def guestimate_init_cond(self, x_data): + """Roughly estimate initial conditions for parameter estimation. + + Args: + x_data: a matrix representing RNA data. + + Returns: + Estimated initial conditions. + """ return guestimate_init_cond(x_data) def guestimate_gamma(self, x_data, time): + """Roughly estimate initial conditions for parameter estimation. + + Args: + x_data: a matrix representing RNA data. + time: a matrix of time information. + + Returns: + Estimated gamma. + """ return guestimate_gamma(x_data, time) def get_param(self, key): + """Get the estimated parameter value by key. + + Args: + key: the key of parameter. + + Returns: + The estimated parameter value. + """ return self.popt[np.where(self.kin_param_keys == key)[0][0]] def calc_half_life(self, key): + """Calculate half-life of a parameter. + + Args: + key: the key of parameter. + + Returns: + The half-life value. + """ return np.log(2) / self.get_param(key) def export_dictionary(self): + """Export parameter estimation results as a dictionary. + + Returns: + Dictionary containing model name, kinetic parameters, and initial conditions. + """ mdl_name = type(self.simulator).__name__ params = self.export_parameters() param_dict = {self.kin_param_keys[i]: params[i] for i in range(len(params))} From cc6f011750147461ea0dcbb8704ab8ca943c79ac Mon Sep 17 00:00:00 2001 From: sichao Date: Fri, 25 Aug 2023 17:36:48 -0400 Subject: [PATCH 14/36] add docstr for the deg estimation --- dynamo/estimation/tsc/estimation_kinetic.py | 95 ++++++++++++++++++++- 1 file changed, 92 insertions(+), 3 deletions(-) diff --git a/dynamo/estimation/tsc/estimation_kinetic.py b/dynamo/estimation/tsc/estimation_kinetic.py index d46b32d08..879940e4a 100755 --- a/dynamo/estimation/tsc/estimation_kinetic.py +++ b/dynamo/estimation/tsc/estimation_kinetic.py @@ -489,17 +489,41 @@ class Estimation_DeterministicDeg(Estimation_Degradation): """ def __init__(self, beta=None, gamma=None, x0=None): + """Initialize the Estimation_DeterministicDeg object. + + Args: + beta: the splicing rate. + gamma: the degradation rate. + x0: the initial conditions. + """ self.kin_param_keys = np.array(["alpha", "beta", "gamma"]) if beta is not None and gamma is not None and x0 is not None: self._initialize(beta, gamma, x0) def _initialize(self, beta, gamma, x0): + """Initialize the parameters to the default value. + + Args: + beta: the splicing rate. + gamma: the degradation rate. + x0: the initial conditions. + """ ranges = np.zeros((2, 2)) ranges[0] = beta * np.ones(2) if np.isscalar(beta) else beta ranges[1] = gamma * np.ones(2) if np.isscalar(gamma) else gamma super().__init__(ranges, x0, Deterministic()) def auto_fit(self, time, x_data, **kwargs): + """Estimate the parameters. + + Args: + time: the time information. + x_data: a matrix representing RNA data. + kwargs: the additional keyword arguments. + + Returns: + The optimized parameters and the cost. + """ be0 = self.guestimate_gamma(x_data[0, :], time) ga0 = self.guestimate_gamma(x_data[0, :] + x_data[1, :], time) x0 = self.guestimate_init_cond(x_data) @@ -516,10 +540,22 @@ class Estimation_DeterministicDegNosp(Estimation_Degradation): """An estimation class for degradation (without splicing) experiments.""" def __init__(self, gamma=None, x0=None): + """Initialize the Estimation_DeterministicDegNosp object. + + Args: + gamma: the degradation rate. + x0: the initial conditions. + """ if gamma is not None and x0 is not None: self._initialize(gamma, x0) def _initialize(self, gamma, x0): + """Initialize the parameters to the default value. + + Args: + gamma: the degradation rate. + x0: the initial conditions. + """ ranges = gamma * np.ones(2) if np.isscalar(gamma) else gamma if np.isscalar(x0) or x0.ndim > 1: x0_ = x0 @@ -528,6 +564,20 @@ def _initialize(self, gamma, x0): super().__init__(ranges, x0_, Deterministic_NoSplicing()) def auto_fit(self, time, x_data, sample_method="lhs", method=None, normalize=False): + """Estimate the parameters. + + Args: + time: the time information. + x_data: a matrix representing RNA data. + sample_method: method used for sampling initial guesses of parameters: + `lhs`: latin hypercube sampling; + `uniform`: uniform random sampling. + method: method used for solving ODEs. See options in simulator classes. + normalize: whether to normalize the + + Returns: + The optimized parameters and the cost. + """ ga0 = self.guestimate_gamma(x_data, time) x0 = self.guestimate_init_cond(x_data[None]) gamma_bound = np.array([0, 1e2 * ga0]) @@ -552,18 +602,34 @@ class Estimation_MomentDeg(Estimation_DeterministicDeg): """ def __init__(self, beta=None, gamma=None, x0=None, include_cov=True): + """Initialize the Estimation_MomentDeg object. + + Args: + beta: the splicing rate. + gamma: the degradation rate. + x0: the initial conditions. + include_cov: whether to consider covariance when estimating. + """ self.kin_param_keys = np.array(["alpha", "beta", "gamma"]) self.include_cov = include_cov if beta is not None and gamma is not None and x0 is not None: self._initialize(beta, gamma, x0) def _initialize(self, beta, gamma, x0): + """Initialize the parameters to the default value. + + Args: + beta: the splicing rate. + gamma: the degradation rate. + x0: the initial conditions. + """ ranges = np.zeros((2, 2)) ranges[0] = beta * np.ones(2) if np.isscalar(beta) else beta ranges[1] = gamma * np.ones(2) if np.isscalar(gamma) else gamma super(Estimation_DeterministicDeg, self).__init__(ranges, x0, Moments_NoSwitching()) def extract_data_from_simulator(self): + """Get corresponding data from the LinearODE class.""" if self.include_cov: ret = np.zeros((5, len(self.simulator.t))) ret[0] = self.simulator.x[:, self.simulator.u] @@ -581,20 +647,43 @@ def extract_data_from_simulator(self): class Estimation_MomentDegNosp(Estimation_Degradation): - """An estimation class for degradation (without splicing) experiments.""" + """An estimation class for degradation (without splicing) experiments. Order of species: , .""" def __init__(self, gamma=None, x0=None): - """An estimation class for degradation (without splicing) experiments. - Order of species: , + """Initialize the Estimation_MomentDeg object. + + Args: + gamma: the degradation rate. + x0: the initial conditions. """ if gamma is not None and x0 is not None: self._initialize(gamma, x0) def _initialize(self, gamma, x0): + """Initialize the parameters to the default value. + + Args: + gamma: the degradation rate. + x0: the initial conditions. + """ ranges = gamma * np.ones(2) if np.isscalar(gamma) else gamma super().__init__(ranges, x0, Moments_NoSwitchingNoSplicing()) def auto_fit(self, time, x_data, sample_method="lhs", method=None, normalize=False): + """Estimate the parameters. + + Args: + time: the time information. + x_data: a matrix representing RNA data. + sample_method: method used for sampling initial guesses of parameters: + `lhs`: latin hypercube sampling; + `uniform`: uniform random sampling. + method: method used for solving ODEs. See options in simulator classes. + normalize: whether to normalize the + + Returns: + The optimized parameters and the cost. + """ ga0 = self.guestimate_gamma(x_data[0, :], time) x0 = self.guestimate_init_cond(x_data) gamma_bound = np.array([0, 1e2 * ga0]) From 2ef169080b543cdafaae6bdf025cb61f81cbf171 Mon Sep 17 00:00:00 2001 From: sichao Date: Fri, 25 Aug 2023 18:04:24 -0400 Subject: [PATCH 15/36] add docstr for kin estimation class --- dynamo/estimation/tsc/estimation_kinetic.py | 77 +++++++++++++++++++++ 1 file changed, 77 insertions(+) diff --git a/dynamo/estimation/tsc/estimation_kinetic.py b/dynamo/estimation/tsc/estimation_kinetic.py index 879940e4a..893c6e7e2 100755 --- a/dynamo/estimation/tsc/estimation_kinetic.py +++ b/dynamo/estimation/tsc/estimation_kinetic.py @@ -707,6 +707,17 @@ class Estimation_MomentKin(kinetic_estimation): """ def __init__(self, a, b, alpha_a, alpha_i, beta, gamma, include_cov=True): + """Initialize the Estimation_MomentKin object. + + Args: + a: switching rate from active promoter state to inactive promoter state. + b: switching rate from inactive promoter state to active promoter state. + alpha_a: transcription rate for active promoter. + alpha_i: transcription rate for inactive promoter. + beta: splicing rate. + gamma: degradation rate. + include_cov: whether to include the covariance when estimating. + """ self.param_keys = np.array(["a", "b", "alpha_a", "alpha_i", "beta", "gamma"]) ranges = np.zeros((6, 2)) ranges[0] = a * np.ones(2) if np.isscalar(a) else a @@ -719,6 +730,7 @@ def __init__(self, a, b, alpha_a, alpha_i, beta, gamma, include_cov=True): self.include_cov = include_cov def extract_data_from_simulator(self): + """Get corresponding data from the LinearODE class.""" if self.include_cov: ret = np.zeros((5, len(self.simulator.t))) ret[0] = self.simulator.get_nu() @@ -735,28 +747,40 @@ def extract_data_from_simulator(self): return ret def get_alpha_a(self): + """Get the transcription rate for active promoter.""" return self.popt[2] def get_alpha_i(self): + """Get the transcription rate for inactive promoter.""" return self.popt[3] def get_alpha(self): + """Get all transcription rates.""" alpha = self.simulator.fbar(self.get_alpha_a(), self.get_alpha_i()) return alpha def get_beta(self): + """Get the splicing rate.""" return self.popt[4] def get_gamma(self): + """Get the degradation rate.""" return self.popt[5] def calc_spl_half_life(self): + """Calculate the half life of splicing.""" return np.log(2) / self.get_beta() def calc_deg_half_life(self): + """Calculate the half life of degradation.""" return np.log(2) / self.get_gamma() def export_dictionary(self): + """Export parameter estimation results as a dictionary. + + Returns: + Dictionary containing model name, kinetic parameters, and initial conditions. + """ mdl_name = type(self.simulator).__name__ params = self.export_parameters() param_dict = {self.param_keys[i]: params[i] for i in range(len(params))} @@ -775,6 +799,15 @@ class Estimation_MomentKinNosp(kinetic_estimation): """ def __init__(self, a, b, alpha_a, alpha_i, gamma): + """Initialize the Estimation_MomentKinNosp object. + + Args: + a: switching rate from active promoter state to inactive promoter state. + b: switching rate from inactive promoter state to active promoter state. + alpha_a: transcription rate for active promoter. + alpha_i: transcription rate for inactive promoter. + gamma: degradation rate. + """ self.param_keys = np.array(["a", "b", "alpha_a", "alpha_i", "gamma"]) ranges = np.zeros((5, 2)) ranges[0] = a * np.ones(2) if np.isscalar(a) else a @@ -785,28 +818,39 @@ def __init__(self, a, b, alpha_a, alpha_i, gamma): super().__init__(ranges, np.zeros((3, 2)), Moments_Nosplicing()) def extract_data_from_simulator(self): + """Get corresponding data from the LinearODE class.""" ret = np.zeros((2, len(self.simulator.t))) ret[0] = self.simulator.get_nu() ret[1] = self.simulator.x[:, self.simulator.uu] return ret def get_alpha_a(self): + """Get the transcription rate for active promoter.""" return self.popt[2] def get_alpha_i(self): + """Get the transcription rate for inactive promoter.""" return self.popt[3] def get_alpha(self): + """Get all transcription rates.""" alpha = self.simulator.fbar(self.get_alpha_a().self.get_alpha_i()) return alpha def get_gamma(self): + """Get the degradation rate.""" return self.popt[4] def calc_deg_half_life(self): + """Calculate the half life of degradation.""" return np.log(2) / self.get_gamma() def export_dictionary(self): + """Export parameter estimation results as a dictionary. + + Returns: + Dictionary containing model name, kinetic parameters, and initial conditions. + """ mdl_name = type(self.simulator).__name__ params = self.export_parameters() param_dict = {self.param_keys[i]: params[i] for i in range(len(params))} @@ -825,6 +869,13 @@ class Estimation_DeterministicKinNosp(kinetic_estimation): """ def __init__(self, alpha, gamma, x0=0): + """Initialize the Estimation_DeterministicKinNosp object. + + Args: + alpha: transcription rate. + gamma: degradation rate. + x0: the initial condition. + """ self.param_keys = np.array(["alpha", "gamma"]) ranges = np.zeros((2, 2)) ranges[0] = alpha * np.ones(2) if np.isscalar(alpha) else alpha @@ -834,15 +885,23 @@ def __init__(self, alpha, gamma, x0=0): super().__init__(ranges, x0, Deterministic_NoSplicing()) def get_alpha(self): + """Get the transcription rate.""" return self.popt[0] def get_gamma(self): + """Get the degradation rate.""" return self.popt[1] def calc_half_life(self, key): + """Calculate the half life.""" return np.log(2) / self.get_param(key) def export_dictionary(self): + """Export parameter estimation results as a dictionary. + + Returns: + Dictionary containing model name, kinetic parameters, and initial conditions. + """ mdl_name = type(self.simulator).__name__ params = self.export_parameters() param_dict = {self.param_keys[i]: params[i] for i in range(len(params))} @@ -861,6 +920,14 @@ class Estimation_DeterministicKin(kinetic_estimation): """ def __init__(self, alpha, beta, gamma, x0=np.zeros(2)): + """Initialize the Estimation_DeterministicKin object. + + Args: + alpha: transcription rate. + beta: splicing rate. + gamma: degradation rate. + x0: the initial condition. + """ self.param_keys = np.array(["alpha", "beta", "gamma"]) ranges = np.zeros((3, 2)) ranges[0] = alpha * np.ones(2) if np.isscalar(alpha) else alpha @@ -873,21 +940,31 @@ def __init__(self, alpha, beta, gamma, x0=np.zeros(2)): super().__init__(ranges, x0, Deterministic()) def get_alpha(self): + """Get the transcription rate.""" return self.popt[0] def get_beta(self): + """Get the splicing rate.""" return self.popt[1] def get_gamma(self): + """Get the degradation rate.""" return self.popt[2] def calc_spl_half_life(self): + """Calculate the half life of splicing.""" return np.log(2) / self.get_beta() def calc_deg_half_life(self): + """Calculate the half life of degradation.""" return np.log(2) / self.get_gamma() def export_dictionary(self): + """Export parameter estimation results as a dictionary. + + Returns: + Dictionary containing model name, kinetic parameters, and initial conditions. + """ mdl_name = type(self.simulator).__name__ params = self.export_parameters() param_dict = {self.param_keys[i]: params[i] for i in range(len(params))} From 9c4877d1fff1e79c6a4f4cfcc12b10179001368f Mon Sep 17 00:00:00 2001 From: sichao Date: Fri, 25 Aug 2023 18:24:20 -0400 Subject: [PATCH 16/36] add docstr for mixture estimation class --- dynamo/estimation/tsc/estimation_kinetic.py | 133 +++++++++++++++++++- 1 file changed, 131 insertions(+), 2 deletions(-) diff --git a/dynamo/estimation/tsc/estimation_kinetic.py b/dynamo/estimation/tsc/estimation_kinetic.py index 893c6e7e2..d97f039aa 100755 --- a/dynamo/estimation/tsc/estimation_kinetic.py +++ b/dynamo/estimation/tsc/estimation_kinetic.py @@ -983,6 +983,16 @@ class Mixture_KinDeg_NoSwitching(kinetic_estimation): """ def __init__(self, model1, model2, alpha=None, gamma=None, x0=None, beta=None): + """Initialize the Mixture_KinDeg_NoSwitching object. + + Args: + model1: the first model to mix. + model2: the second model to mix. + alpha: transcription rate. + gamma: degradation rate. + x0: the initial condition. + beta: splicing rate. + """ self.model1 = model1 self.model2 = model2 self.scale = 1 @@ -990,6 +1000,14 @@ def __init__(self, model1, model2, alpha=None, gamma=None, x0=None, beta=None): self._initialize(alpha, gamma, x0, beta) def _initialize(self, alpha, gamma, x0, beta=None): + """Initialize the parameters to the default value. + + Args: + alpha: transcription rate. + gamma: degradation rate. + x0: the initial condition. + beta: splicing rate. + """ if type(self.model1) in nosplicing_models: self.param_distributor = [[0, 2], [1, 2]] self.param_keys = ["alpha", "alpha_2", "gamma"] @@ -1010,6 +1028,16 @@ def _initialize(self, alpha, gamma, x0, beta=None): super().__init__(ranges, x0_, model) def normalize_deg_data(self, x_data, weight): + """Normalize the degradation data while preserving the relative proportions between species. It calculates + scaling factors to ensure the data's range remains within a certain limit. + + Args: + x_data: a matrix representing RNA data. + weight: weight for scaling. + + Returns: + A tuple containing the normalized degradation data and the scaling factor. + """ x_data_norm = np.array(x_data, copy=True) x_data_kin = x_data_norm[: self.model1.n_species, :] @@ -1022,6 +1050,21 @@ def normalize_deg_data(self, x_data, weight): return x_data_norm, scale def auto_fit(self, time, x_data, alpha_min=0.1, beta_min=50, gamma_min=10, kin_weight=2, use_p0=True, **kwargs): + """Estimate the parameters. + + Args: + time: the time information. + x_data: a matrix representing RNA data. + alpha_min: the minimum limitation on transcription rate. + beta_min: the minimum limitation on splicing rate. + gamma_min: the minimum limitation on degradation rate. + kin_weight: weight for scaling during normalization. + use_p0: whether to use initial parameters when estimating. + kwargs: the additional keyword arguments. + + Returns: + The optimized parameters and the cost. + """ if kin_weight is not None: x_data_norm, self.scale = self.normalize_deg_data(x_data, kin_weight) else: @@ -1059,17 +1102,35 @@ def auto_fit(self, time, x_data, alpha_min=0.1, beta_min=50, gamma_min=10, kin_w return popt, cost def export_model(self, reinstantiate=True): + """Export the mixture model. + + Args: + reinstantiate: whether to reinstantiate the model. + + Returns: + MixtureModels or LinearODE. + """ if reinstantiate: return MixtureModels([self.model1, self.model2], self.param_distributor) else: return self.simulator def export_x0(self): + """Export optimized initial conditions for the mixture of models analysis. + + Returns: + Exported initial conditions. + """ x = self.get_opt_x0_params() x[self.model1.n_species :] *= self.scale return x def export_dictionary(self): + """Export parameter estimation results as a dictionary. + + Returns: + Dictionary containing model nameS, kinetic parameters, and initial conditions. + """ mdl1_name = type(self.model1).__name__ mdl2_name = type(self.model2).__name__ params = self.export_parameters() @@ -1099,6 +1160,17 @@ def __init__( x0=None, beta=None, ): + """Initialize the Lambda_NoSwitching object. + + Args: + model1: the first model to mix. + model2: the second model to mix. + alpha: transcription rate. + lambd: the lambd value. + gamma: degradation rate. + x0: the initial condition. + beta: splicing rate. + """ self.model1 = model1 self.model2 = model2 self.scale = 1 @@ -1106,8 +1178,13 @@ def __init__( self._initialize(alpha, gamma, x0, beta) def _initialize(self, alpha, gamma, x0, beta=None): - """ - parameter order: alpha, lambda, (beta), gamma + """Initialize the parameters to the default value. + + Args: + alpha: transcription rate. + gamma: degradation rate. + x0: the initial condition. + beta: splicing rate. """ if type(self.model1) in nosplicing_models and type(self.model2) in nosplicing_models: self.param_keys = ["alpha", "lambda", "gamma"] @@ -1127,9 +1204,27 @@ def _initialize(self, alpha, gamma, x0, beta=None): super(Mixture_KinDeg_NoSwitching, self).__init__(ranges, x0_, model) def auto_fit(self, time, x_data, **kwargs): + """Estimate the parameters. + + Args: + time: the time information. + x_data: a matrix representing RNA data. + kwargs: the additional keyword arguments. + + Returns: + The optimized parameters and the cost. + """ return super().auto_fit(time, x_data, kin_weight=None, use_p0=False, **kwargs) def export_model(self, reinstantiate=True): + """Export the mixture model. + + Args: + reinstantiate: whether to reinstantiate the model. + + Returns: + MixtureModels or LinearODE. + """ if reinstantiate: return LambdaModels_NoSwitching(self.model1, self.model2) else: @@ -1137,18 +1232,43 @@ def export_model(self, reinstantiate=True): class Estimation_KineticChase(kinetic_estimation): + """An estimation class for kinetic chase experiment.""" def __init__(self, alpha=None, gamma=None, x0=None): + """Initialize the Estimation_KineticChase object. + + Args: + alpha: transcription rate. + gamma: degradation rate. + x0: the initial condition. + """ self.kin_param_keys = np.array(["alpha", "gamma"]) if alpha is not None and gamma is not None and x0 is not None: self._initialize(alpha, gamma, x0) def _initialize(self, alpha, gamma, x0): + """Initialize the parameters to the default value. + + Args: + alpha: transcription rate. + gamma: degradation rate. + x0: the initial condition. + """ ranges = np.zeros((2, 2)) ranges[0] = alpha * np.ones(2) if np.isscalar(alpha) else alpha ranges[1] = gamma * np.ones(2) if np.isscalar(gamma) else gamma super().__init__(ranges, np.atleast_2d(x0), KineticChase()) def auto_fit(self, time, x_data, **kwargs): + """Estimate the parameters. + + Args: + time: the time information. + x_data: a matrix representing RNA data. + kwargs: the additional keyword arguments. + + Returns: + The optimized parameters and the cost. + """ if len(time) != len(np.unique(time)): t = np.unique(time) x = strat_mom(x_data, time, np.mean) @@ -1163,18 +1283,27 @@ def auto_fit(self, time, x_data, **kwargs): return popt, cost def get_param(self, key): + """Get corresponding parameter according to the key.""" return self.popt[np.where(self.kin_param_keys == key)[0][0]] def get_alpha(self): + """Get the transcription rate.""" return self.popt[0] def get_gamma(self): + """Get the degradation rate.""" return self.popt[1] def calc_half_life(self, key): + """Calculate the half life.""" return np.log(2) / self.get_param(key) def export_dictionary(self): + """Export parameter estimation results as a dictionary. + + Returns: + Dictionary containing model nameS, kinetic parameters, and initial conditions. + """ mdl_name = type(self.simulator).__name__ params = self.export_parameters() param_dict = {self.kin_param_keys[i]: params[i] for i in range(len(params))} From d65ffcd35d98384a7571d2301a69ee29ff5ae6cc Mon Sep 17 00:00:00 2001 From: sichao Date: Fri, 25 Aug 2023 18:39:51 -0400 Subject: [PATCH 17/36] add docstr for goodness of fit --- dynamo/estimation/tsc/estimation_kinetic.py | 58 +++++++++++++++++++++ 1 file changed, 58 insertions(+) diff --git a/dynamo/estimation/tsc/estimation_kinetic.py b/dynamo/estimation/tsc/estimation_kinetic.py index d97f039aa..6307fee22 100755 --- a/dynamo/estimation/tsc/estimation_kinetic.py +++ b/dynamo/estimation/tsc/estimation_kinetic.py @@ -1317,7 +1317,19 @@ def export_dictionary(self): class GoodnessOfFit: + """Evaluate goodness of fitting. + + This class provides methods for assessing the quality of predictions, using various metrics including Gaussian + likelihood, Gaussian log-likelihood, and mean squared deviation. + """ def __init__(self, simulator, params=None, x0=None): + """Initialize the GoodnessOfFit object. + + Args: + simulator: the linearODE class. + params: the parameters. + x0: the initial conditions. + """ self.simulator = simulator if params is not None: self.simulator.set_params(*params) @@ -1329,6 +1341,14 @@ def __init__(self, simulator, params=None, x0=None): self.pred = None def extract_data_from_simulator(self, species=None): + """Extract data from the simulator's results. + + Args: + species: index of the species to extract. + + Returns: + Extracted data from the simulator's results. + """ ret = self.simulator.x.T if species is not None: ret = ret[species] @@ -1343,6 +1363,16 @@ def prepare_data( normalize=True, reintegrate=True, ): + """Prepare data for evaluation. + + Args: + t: the time information. + x_data: the RNA data. + species: index of the species to consider. + method: integration method. + normalize: whether to normalize data. + reintegrate: whether to reintegrate the model. + """ if reintegrate: self.simulator.integrate(t, method=method) x_model = self.extract_data_from_simulator(species=species) @@ -1361,12 +1391,27 @@ def prepare_data( self.pred = strat_mom(x_model_norm.T, t, np.mean) def normalize(self, x_data, x_model, scale=None): + """Normalize data and model predictions. + + Args: + x_data: the RNA data. + x_model: predictions from model. + scale: Scaling factors for normalization. + + Returns: + Normalized observations and model predictions. + """ scale = np.max(x_data, 1) if scale is None else scale x_data_norm = (x_data.T / scale).T x_model_norm = (x_model.T / scale).T return x_data_norm, x_model_norm def calc_gaussian_likelihood(self): + """ Calculate the Gaussian likelihood between model predictions and observations. + + Returns: + Gaussian likelihood value. + """ sig = np.array(self.sigm, copy=True) if np.any(sig == 0): main_warning("Some standard deviations are 0; Set to 1 instead.") @@ -1376,6 +1421,11 @@ def calc_gaussian_likelihood(self): return ret def calc_gaussian_loglikelihood(self): + """Calculate the Gaussian log-likelihood between model predictions and observations. + + Returns: + Gaussian log-likelihood value. + """ sig = np.array(self.sigm, copy=True) if np.any(sig == 0): main_warning("Some standard deviations are 0; Set to 1 instead.") @@ -1385,6 +1435,14 @@ def calc_gaussian_loglikelihood(self): return ret def calc_mean_squared_deviation(self, weighted=True): + """Calculate the mean squared deviation between model predictions and observations. + + Args: + weighted: whether to weight the output. + + Returns: + Mean squared deviation + """ sig = np.array(self.sigm, copy=True) if np.any(sig == 0): main_warning("Some standard deviations are 0; Set to 1 instead.") From 4114133e4ecb6d4f98c5ad6c19e1e7f335809b53 Mon Sep 17 00:00:00 2001 From: Sichao25 Date: Sat, 26 Aug 2023 14:47:15 -0400 Subject: [PATCH 18/36] add typing for Deg experiments --- dynamo/estimation/tsc/estimation_kinetic.py | 63 ++++++++++++++------- 1 file changed, 44 insertions(+), 19 deletions(-) diff --git a/dynamo/estimation/tsc/estimation_kinetic.py b/dynamo/estimation/tsc/estimation_kinetic.py index 6307fee22..4aeaa17a6 100755 --- a/dynamo/estimation/tsc/estimation_kinetic.py +++ b/dynamo/estimation/tsc/estimation_kinetic.py @@ -1,4 +1,4 @@ -from typing import List, Optional, Tuple, Union +from typing import Dict, List, Optional, Tuple, Union import warnings import numpy as np @@ -409,7 +409,7 @@ def test_chi2( class Estimation_Degradation(kinetic_estimation): """The base parameters, estimation class for degradation experiments.""" - def __init__(self, ranges, x0, simulator): + def __init__(self, ranges: np.ndarray, x0: np.ndarray, simulator: LinearODE) -> None: """Initialize the Estimation_Degradation object. Args: @@ -420,7 +420,7 @@ def __init__(self, ranges, x0, simulator): self.kin_param_keys = np.array(["alpha", "gamma"]) super().__init__(np.vstack((np.zeros(2), ranges)), x0, simulator) - def guestimate_init_cond(self, x_data): + def guestimate_init_cond(self, x_data: np.ndarray) -> np.ndarray: """Roughly estimate initial conditions for parameter estimation. Args: @@ -431,7 +431,7 @@ def guestimate_init_cond(self, x_data): """ return guestimate_init_cond(x_data) - def guestimate_gamma(self, x_data, time): + def guestimate_gamma(self, x_data: np.ndarray, time: np.ndarray) -> np.ndarray: """Roughly estimate initial conditions for parameter estimation. Args: @@ -443,7 +443,7 @@ def guestimate_gamma(self, x_data, time): """ return guestimate_gamma(x_data, time) - def get_param(self, key): + def get_param(self, key: str) -> np.ndarray: """Get the estimated parameter value by key. Args: @@ -454,7 +454,7 @@ def get_param(self, key): """ return self.popt[np.where(self.kin_param_keys == key)[0][0]] - def calc_half_life(self, key): + def calc_half_life(self, key: str) -> np.ndarray: """Calculate half-life of a parameter. Args: @@ -465,7 +465,7 @@ def calc_half_life(self, key): """ return np.log(2) / self.get_param(key) - def export_dictionary(self): + def export_dictionary(self) -> Dict: """Export parameter estimation results as a dictionary. Returns: @@ -488,7 +488,12 @@ class Estimation_DeterministicDeg(Estimation_Degradation): Order of species: , """ - def __init__(self, beta=None, gamma=None, x0=None): + def __init__( + self, + beta: Optional[np.ndarray] = None, + gamma: Optional[np.ndarray] = None, + x0: Optional[np.ndarray] = None, + ) -> None: """Initialize the Estimation_DeterministicDeg object. Args: @@ -500,7 +505,7 @@ def __init__(self, beta=None, gamma=None, x0=None): if beta is not None and gamma is not None and x0 is not None: self._initialize(beta, gamma, x0) - def _initialize(self, beta, gamma, x0): + def _initialize(self, beta: np.ndarray, gamma: np.ndarray, x0: np.ndarray) -> None: """Initialize the parameters to the default value. Args: @@ -513,7 +518,7 @@ def _initialize(self, beta, gamma, x0): ranges[1] = gamma * np.ones(2) if np.isscalar(gamma) else gamma super().__init__(ranges, x0, Deterministic()) - def auto_fit(self, time, x_data, **kwargs): + def auto_fit(self, time: np.ndarray, x_data: np.ndarray, **kwargs) -> Tuple[np.ndarray, float]: """Estimate the parameters. Args: @@ -539,7 +544,7 @@ def auto_fit(self, time, x_data, **kwargs): class Estimation_DeterministicDegNosp(Estimation_Degradation): """An estimation class for degradation (without splicing) experiments.""" - def __init__(self, gamma=None, x0=None): + def __init__(self, gamma: Optional[np.ndarray] = None, x0: Optional[np.ndarray] = None) -> None: """Initialize the Estimation_DeterministicDegNosp object. Args: @@ -549,7 +554,7 @@ def __init__(self, gamma=None, x0=None): if gamma is not None and x0 is not None: self._initialize(gamma, x0) - def _initialize(self, gamma, x0): + def _initialize(self, gamma: np.ndarray, x0: np.ndarray) -> None: """Initialize the parameters to the default value. Args: @@ -563,7 +568,14 @@ def _initialize(self, gamma, x0): x0_ = np.array([x0]) super().__init__(ranges, x0_, Deterministic_NoSplicing()) - def auto_fit(self, time, x_data, sample_method="lhs", method=None, normalize=False): + def auto_fit( + self, + time: np.ndarray, + x_data: np.ndarray, + sample_method: str = "lhs", + method: Optional[str] = None, + normalize: bool = False, + ) -> Tuple[np.ndarray, float]: """Estimate the parameters. Args: @@ -601,7 +613,13 @@ class Estimation_MomentDeg(Estimation_DeterministicDeg): Order of parameters: beta, gamma """ - def __init__(self, beta=None, gamma=None, x0=None, include_cov=True): + def __init__( + self, + beta: Optional[np.ndarray] = None, + gamma: Optional[np.ndarray] = None, + x0: Optional[np.ndarray] = None, + include_cov: bool = True, + ) -> None: """Initialize the Estimation_MomentDeg object. Args: @@ -615,7 +633,7 @@ def __init__(self, beta=None, gamma=None, x0=None, include_cov=True): if beta is not None and gamma is not None and x0 is not None: self._initialize(beta, gamma, x0) - def _initialize(self, beta, gamma, x0): + def _initialize(self, beta: np.ndarray, gamma: np.ndarray, x0: np.ndarray) -> None: """Initialize the parameters to the default value. Args: @@ -628,7 +646,7 @@ def _initialize(self, beta, gamma, x0): ranges[1] = gamma * np.ones(2) if np.isscalar(gamma) else gamma super(Estimation_DeterministicDeg, self).__init__(ranges, x0, Moments_NoSwitching()) - def extract_data_from_simulator(self): + def extract_data_from_simulator(self) -> np.ndarray: """Get corresponding data from the LinearODE class.""" if self.include_cov: ret = np.zeros((5, len(self.simulator.t))) @@ -649,7 +667,7 @@ def extract_data_from_simulator(self): class Estimation_MomentDegNosp(Estimation_Degradation): """An estimation class for degradation (without splicing) experiments. Order of species: , .""" - def __init__(self, gamma=None, x0=None): + def __init__(self, gamma: Optional[np.ndarray] = None, x0: Optional[np.ndarray] = None) -> None: """Initialize the Estimation_MomentDeg object. Args: @@ -659,7 +677,7 @@ def __init__(self, gamma=None, x0=None): if gamma is not None and x0 is not None: self._initialize(gamma, x0) - def _initialize(self, gamma, x0): + def _initialize(self, gamma: np.ndarray, x0: np.ndarray) -> None: """Initialize the parameters to the default value. Args: @@ -669,7 +687,14 @@ def _initialize(self, gamma, x0): ranges = gamma * np.ones(2) if np.isscalar(gamma) else gamma super().__init__(ranges, x0, Moments_NoSwitchingNoSplicing()) - def auto_fit(self, time, x_data, sample_method="lhs", method=None, normalize=False): + def auto_fit( + self, + time: np.ndarray, + x_data: np.ndarray, + sample_method: str = "lhs", + method: Optional[str] = None, + normalize: bool = False, + ) -> Tuple[np.ndarray, float]: """Estimate the parameters. Args: From 279f67aa605d4353cd5189cbd14834db7bb48e23 Mon Sep 17 00:00:00 2001 From: sichao Date: Mon, 28 Aug 2023 10:28:23 -0400 Subject: [PATCH 19/36] debug estimation_kinetic docstr --- dynamo/estimation/tsc/estimation_kinetic.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/dynamo/estimation/tsc/estimation_kinetic.py b/dynamo/estimation/tsc/estimation_kinetic.py index 4aeaa17a6..c420a5095 100755 --- a/dynamo/estimation/tsc/estimation_kinetic.py +++ b/dynamo/estimation/tsc/estimation_kinetic.py @@ -585,7 +585,7 @@ def auto_fit( `lhs`: latin hypercube sampling; `uniform`: uniform random sampling. method: method used for solving ODEs. See options in simulator classes. - normalize: whether to normalize the + normalize: whether to normalize the data. Returns: The optimized parameters and the cost. @@ -704,7 +704,7 @@ def auto_fit( `lhs`: latin hypercube sampling; `uniform`: uniform random sampling. method: method used for solving ODEs. See options in simulator classes. - normalize: whether to normalize the + normalize: whether to normalize the data. Returns: The optimized parameters and the cost. @@ -1466,7 +1466,7 @@ def calc_mean_squared_deviation(self, weighted=True): weighted: whether to weight the output. Returns: - Mean squared deviation + Mean squared deviation. """ sig = np.array(self.sigm, copy=True) if np.any(sig == 0): From be91ce2b324b45714a7ae3295cb128e69c45ede9 Mon Sep 17 00:00:00 2001 From: sichao Date: Mon, 28 Aug 2023 11:04:20 -0400 Subject: [PATCH 20/36] add docstr and typing for lambda models_noswitching --- dynamo/estimation/tsc/utils_kinetic.py | 20 +++++++++++++++----- 1 file changed, 15 insertions(+), 5 deletions(-) diff --git a/dynamo/estimation/tsc/utils_kinetic.py b/dynamo/estimation/tsc/utils_kinetic.py index aedb21ae2..f2ac5d33f 100755 --- a/dynamo/estimation/tsc/utils_kinetic.py +++ b/dynamo/estimation/tsc/utils_kinetic.py @@ -217,10 +217,15 @@ def set_params(self, *params: Tuple) -> None: class LambdaModels_NoSwitching(MixtureModels): - def __init__(self, model1, model2): - """ - parameter order: alpha, lambda, (beta), gamma - distributor order: alpha_1, alpha_2, (beta), gamma + """Linear ODEs for the lambda mixture model. The order of params is: + parameter order: alpha, lambda, (beta), gamma + distributor order: alpha_1, alpha_2, (beta), gamma""" + def __init__(self, model1: LinearODE, model2: LinearODE) -> None: + """Initialize the LambdaModels_NoSwitching class. + + Args: + model1: the first model to mix. + model2: the second model to mix. """ models = [model1, model2] if type(model1) in nosplicing_models and type(model2) in nosplicing_models: @@ -231,7 +236,12 @@ def __init__(self, model1, model2): param_distributor = [dist1, dist2] super().__init__(models, param_distributor) - def param_mixer(self, *params): + def param_mixer(self, *params) -> np.ndarray: + """Set parameters for all models. + + Args: + params: tuple of parameters. + """ lam = params[1] alp_1 = params[0] * lam alp_2 = params[0] * (1 - lam) From bf70c3c7bdd0148ff4af47616f16d5039467dc08 Mon Sep 17 00:00:00 2001 From: sichao Date: Mon, 28 Aug 2023 12:59:43 -0400 Subject: [PATCH 21/36] add docstr for determinstic model --- dynamo/estimation/tsc/utils_kinetic.py | 97 ++++++++++++++++++++++++-- 1 file changed, 93 insertions(+), 4 deletions(-) diff --git a/dynamo/estimation/tsc/utils_kinetic.py b/dynamo/estimation/tsc/utils_kinetic.py index f2ac5d33f..26f278fa9 100755 --- a/dynamo/estimation/tsc/utils_kinetic.py +++ b/dynamo/estimation/tsc/utils_kinetic.py @@ -691,9 +691,16 @@ def computeKnp(self): class Deterministic(LinearODE): + """This class simulates the deterministic dynamics of a transcription-splicing system.""" def __init__(self, alpha=None, beta=None, gamma=None, x0=None): - """This class simulates the deterministic dynamics of - a transcription-splicing system.""" + """Initialize the Deterministic object. + + Args: + alpha: the transcription rate. + beta: the splicing rate. + gamma: the degradation rate. + x0: the inital conditions. + """ # species self.u = 0 self.s = 1 @@ -711,6 +718,17 @@ def __init__(self, alpha=None, beta=None, gamma=None, x0=None): self.set_params(alpha, beta, gamma) def ode_func(self, x, t): + """The ODE functions to solve: + dx[u] = alpha - beta * u + dx[v] = beta * u - gamma * s + + Args: + x: the x variable. + t: the time information. + + Returns: + An array containing the ODEs' output. + """ dx = np.zeros(len(x)) # parameters al = self.al @@ -724,6 +742,13 @@ def ode_func(self, x, t): return dx def set_params(self, alpha, beta, gamma): + """Set the parameters. + + Args: + alpha: the transcription rate. + beta: the splicing rate. + gamma: the degradation rate. + """ self.al = alpha self.be = beta self.ga = gamma @@ -732,6 +757,13 @@ def set_params(self, alpha, beta, gamma): super().reset() def integrate(self, t, x0=None, method="analytical"): + """Integrate the ODE using the given time values. + + Args: + t: array of time values. + x0: array of initial conditions. + method: the method to integrate, including "matrix" and "numerical". + """ method = self.default_method if method is None else method if method == "analytical": sol = self.integrate_analytical(t, x0) @@ -741,6 +773,11 @@ def integrate(self, t, x0=None, method="analytical"): self.t = t def computeKnp(self): + """Calculate the K and p from ODE function such that dx = Kx + p. + + Returns: + A tuple containing K and p. + """ # parameters al = self.al be = self.be @@ -760,6 +797,15 @@ def computeKnp(self): return K, p def integrate_analytical(self, t, x0=None): + """Integrate the odes with the analytical solution. + + Args: + t: the time information. + x0: the initial conditions. + + Returns: + The solution of unspliced and splcied mRNA wrapped in an array. + """ if x0 is None: x0 = self.x0 else: @@ -770,9 +816,15 @@ def integrate_analytical(self, t, x0=None): class Deterministic_NoSplicing(LinearODE): + """The class simulates the deterministic dynamics of a transcription-splicing system.""" def __init__(self, alpha=None, gamma=None, x0=None): - """This class simulates the deterministic dynamics of - a transcription-splicing system.""" + """Initialize the Deterministic_NoSplicing object. + + Args: + alpha: the transcription rate. + gamma: the degradation rate. + x0: the initial conditions. + """ # species self.u = 0 @@ -789,6 +841,16 @@ def __init__(self, alpha=None, gamma=None, x0=None): self.set_params(alpha, gamma) def ode_func(self, x, t): + """The ODE functions to solve: + dx[u] = alpha - gamma * u + + Args: + x: the x variable. + t: the time information. + + Returns: + An array containing the ODEs' output. + """ dx = np.zeros(len(x)) # parameters al = self.al @@ -800,6 +862,12 @@ def ode_func(self, x, t): return dx def set_params(self, alpha, gamma): + """Set the parameters. + + Args: + alpha: the transcription rate. + gamma: the degradation rate. + """ self.al = alpha self.ga = gamma @@ -807,6 +875,13 @@ def set_params(self, alpha, gamma): super().reset() def integrate(self, t, x0=None, method="analytical"): + """Integrate the ODE using the given time values. + + Args: + t: array of time values. + x0: array of initial conditions. + method: the method to integrate, including "matrix" and "numerical". + """ method = self.default_method if method is None else method if method == "analytical": sol = self.integrate_analytical(t, x0) @@ -816,6 +891,11 @@ def integrate(self, t, x0=None, method="analytical"): self.t = t def computeKnp(self): + """Calculate the K and p from ODE function such that dx = Kx + p. + + Returns: + A tuple containing K and p. + """ # parameters al = self.al ga = self.ga @@ -830,6 +910,15 @@ def computeKnp(self): return K, p def integrate_analytical(self, t, x0=None): + """Integrate the odes with the analytical solution. + + Args: + t: the time information. + x0: the initial conditions. + + Returns: + The solution of unspliced mRNA as an array. + """ x0 = self.x0 if x0 is None else x0 x0 = x0 if np.isscalar(x0) else x0[self.u] if self.x0 is None: From d9b033764b73ecf98bbb1635b11d2e4f911309b3 Mon Sep 17 00:00:00 2001 From: sichao Date: Mon, 28 Aug 2023 13:26:08 -0400 Subject: [PATCH 22/36] add typing for deterministic odes --- dynamo/estimation/tsc/utils_kinetic.py | 35 +++++++++++++++++--------- 1 file changed, 23 insertions(+), 12 deletions(-) diff --git a/dynamo/estimation/tsc/utils_kinetic.py b/dynamo/estimation/tsc/utils_kinetic.py index 26f278fa9..2a55e0ceb 100755 --- a/dynamo/estimation/tsc/utils_kinetic.py +++ b/dynamo/estimation/tsc/utils_kinetic.py @@ -692,7 +692,13 @@ def computeKnp(self): class Deterministic(LinearODE): """This class simulates the deterministic dynamics of a transcription-splicing system.""" - def __init__(self, alpha=None, beta=None, gamma=None, x0=None): + def __init__( + self, + alpha: Optional[np.ndarray] = None, + beta: Optional[np.ndarray] = None, + gamma: Optional[np.ndarray] = None, + x0: Optional[np.ndarray] = None, + ) -> None: """Initialize the Deterministic object. Args: @@ -717,7 +723,7 @@ def __init__(self, alpha=None, beta=None, gamma=None, x0=None): if not (alpha is None or beta is None or gamma is None): self.set_params(alpha, beta, gamma) - def ode_func(self, x, t): + def ode_func(self, x: np.ndarray, t: np.ndarray) -> np.ndarray: """The ODE functions to solve: dx[u] = alpha - beta * u dx[v] = beta * u - gamma * s @@ -741,7 +747,7 @@ def ode_func(self, x, t): return dx - def set_params(self, alpha, beta, gamma): + def set_params(self, alpha: np.ndarray, beta: np.ndarray, gamma: np.ndarray) -> None: """Set the parameters. Args: @@ -756,7 +762,7 @@ def set_params(self, alpha, beta, gamma): # reset solutions super().reset() - def integrate(self, t, x0=None, method="analytical"): + def integrate(self, t: np.ndarray, x0: Optional[np.ndarray] = None, method: str = "analytical") -> None: """Integrate the ODE using the given time values. Args: @@ -772,7 +778,7 @@ def integrate(self, t, x0=None, method="analytical"): self.x = sol self.t = t - def computeKnp(self): + def computeKnp(self) -> Tuple[np.ndarray, np.ndarray]: """Calculate the K and p from ODE function such that dx = Kx + p. Returns: @@ -796,7 +802,7 @@ def computeKnp(self): return K, p - def integrate_analytical(self, t, x0=None): + def integrate_analytical(self, t: np.ndarray, x0: Optional[np.ndarray] = None) -> np.ndarray: """Integrate the odes with the analytical solution. Args: @@ -817,7 +823,12 @@ def integrate_analytical(self, t, x0=None): class Deterministic_NoSplicing(LinearODE): """The class simulates the deterministic dynamics of a transcription-splicing system.""" - def __init__(self, alpha=None, gamma=None, x0=None): + def __init__( + self, + alpha: Optional[np.ndarray] = None, + gamma: Optional[np.ndarray] = None, + x0: Optional[np.ndarray] = None, + ) -> None: """Initialize the Deterministic_NoSplicing object. Args: @@ -840,7 +851,7 @@ def __init__(self, alpha=None, gamma=None, x0=None): if not (alpha is None or gamma is None): self.set_params(alpha, gamma) - def ode_func(self, x, t): + def ode_func(self, x: np.ndarray, t: np.ndarray) -> np.ndarray: """The ODE functions to solve: dx[u] = alpha - gamma * u @@ -861,7 +872,7 @@ def ode_func(self, x, t): return dx - def set_params(self, alpha, gamma): + def set_params(self, alpha: np.ndarray, gamma: np.ndarray) -> np.ndarray: """Set the parameters. Args: @@ -874,7 +885,7 @@ def set_params(self, alpha, gamma): # reset solutions super().reset() - def integrate(self, t, x0=None, method="analytical"): + def integrate(self, t: np.ndarray, x0: Optional[np.ndarray] = None, method: str = "analytical") -> None: """Integrate the ODE using the given time values. Args: @@ -890,7 +901,7 @@ def integrate(self, t, x0=None, method="analytical"): self.x = sol self.t = t - def computeKnp(self): + def computeKnp(self) -> Tuple[np.ndarray, np.ndarray]: """Calculate the K and p from ODE function such that dx = Kx + p. Returns: @@ -909,7 +920,7 @@ def computeKnp(self): return K, p - def integrate_analytical(self, t, x0=None): + def integrate_analytical(self, t: np.ndarray, x0: Optional[np.ndarray] = None) -> np.ndarray: """Integrate the odes with the analytical solution. Args: From c1d844d6ac2101b01cb3997df2dda611ade6d859 Mon Sep 17 00:00:00 2001 From: sichao Date: Mon, 28 Aug 2023 16:15:08 -0400 Subject: [PATCH 23/36] add docstr for Moments linearODEs --- dynamo/estimation/tsc/utils_kinetic.py | 273 ++++++++++++++++++++++++- 1 file changed, 265 insertions(+), 8 deletions(-) diff --git a/dynamo/estimation/tsc/utils_kinetic.py b/dynamo/estimation/tsc/utils_kinetic.py index 2a55e0ceb..1fb46a9de 100755 --- a/dynamo/estimation/tsc/utils_kinetic.py +++ b/dynamo/estimation/tsc/utils_kinetic.py @@ -250,6 +250,8 @@ def param_mixer(self, *params) -> np.ndarray: class Moments(LinearODE): + """The class simulates the dynamics of first and second moments of a transcription-splicing system with promoter + switching.""" def __init__( self, a=None, @@ -260,8 +262,17 @@ def __init__( gamma=None, x0=None, ): - """This class simulates the dynamics of first and second moments of - a transcription-splicing system with promoter switching.""" + """Initialize the Moments object. + + Args: + a: switching rate from active promoter state to inactive promoter state. + b: switching rate from inactive promoter state to active promoter state. + alpha_a: transcription rate for active promoter. + alpha_i: transcription rate for inactive promoter. + beta: splicing rate. + gamma: degradation rate. + x0: the initial conditions. + """ # species self.ua = 0 self.ui = 1 @@ -281,6 +292,20 @@ def __init__( self.set_params(a, b, alpha_a, alpha_i, beta, gamma) def ode_func(self, x, t): + """ODE functions to solve: + dx[u_a] = alpha_a - beta * u_a + a * (u_i - u_a) + dx[u_i] = ai - beta * u_i - b * (u_i - u_a) + dx[x_a] = beta * u_a - ga * x_a + a * (x_i - x_a) + dx[x_i] = beta * u_i - ga * x_i - b * (x_i - x_a) + The second moments is calculated from the variance and covariance of variable u and x. + + Args: + x: the variable. + t: the aray of time. + + Returns: + The derivatives dx. + """ dx = np.zeros(len(x)) # parameters a = self.a @@ -304,9 +329,28 @@ def ode_func(self, x, t): return dx def fbar(self, x_a, x_i): + """Calculate the count of a variable by averaging active and inactive states. + + Args: + x_a: the variable x under the active state. + x_i: the variable x under the inactive state. + + Returns: + The count of variable x. + """ return self.b / (self.a + self.b) * x_a + self.a / (self.a + self.b) * x_i def set_params(self, a, b, alpha_a, alpha_i, beta, gamma): + """Set the parameters. + + Args: + a: switching rate from active promoter state to inactive promoter state. + b: switching rate from inactive promoter state to active promoter state. + alpha_a: transcription rate for active promoter. + alpha_i: transcription rate for inactive promoter. + beta: splicing rate. + gamma: degradation rate. + """ self.a = a self.b = b self.aa = alpha_a @@ -318,6 +362,11 @@ def set_params(self, a, b, alpha_a, alpha_i, beta, gamma): super().reset() def get_all_central_moments(self): + """Get the first and second central moments for all variables. + + Returns: + An array containing all central moments. + """ ret = np.zeros((4, len(self.t))) ret[0] = self.get_nu() ret[1] = self.get_nx() @@ -326,37 +375,82 @@ def get_all_central_moments(self): return ret def get_nosplice_central_moments(self): + """Get the central moments for labeled data. + + Returns: + The central moments. + """ ret = np.zeros((2, len(self.t))) ret[0] = self.get_n_labeled() ret[1] = self.get_var_labeled() return ret def get_nu(self): + """Get the number of the variable u from the mean averaging active and inactive state. + + Returns: + The number of the variable u. + """ return self.fbar(self.x[:, self.ua], self.x[:, self.ui]) def get_nx(self): + """Get the number of the variable x from the mean averaging active and inactive state. + + Returns: + The number of the variable x. + """ return self.fbar(self.x[:, self.xa], self.x[:, self.xi]) def get_n_labeled(self): + """Get the number of the labeled data by combining the count of two variables. + + Returns: + The number of the labeled data. + """ return self.get_nu() + self.get_nx() def get_var_nu(self): + """Get the variance of the variable u. + + Returns: + The variance of the variable u. + """ c = self.get_nu() return self.x[:, self.uu] + c - c**2 def get_var_nx(self): + """Get the variance of the variable x. + + Returns: + The variance of the variable x. + """ c = self.get_nx() return self.x[:, self.xx] + c - c**2 def get_cov_ux(self): + """Get the covariance of the variable x and u. + + Returns: + The covariance. + """ cu = self.get_nu() cx = self.get_nx() return self.x[:, self.ux] - cu * cx def get_var_labeled(self): + """Get the variance of the labeled data by combining the variance of two variables and covariance. + + Returns: + The variance of the labeled data. + """ return self.get_var_nu() + self.get_var_nx() + 2 * self.get_cov_ux() def computeKnp(self): + """Calculate the K and p from ODE function such that dx = Kx + p. + + Returns: + A tuple containing K and p. + """ # parameters a = self.a b = self.b @@ -411,9 +505,19 @@ def computeKnp(self): class Moments_Nosplicing(LinearODE): + """The class simulates the dynamics of first and second moments of a transcription-splicing system with promoter + switching.""" def __init__(self, a=None, b=None, alpha_a=None, alpha_i=None, gamma=None, x0=None): - """This class simulates the dynamics of first and second moments of - a transcription-splicing system with promoter switching.""" + """Initialize the Moments_Nosplicing object. + + Args: + a: switching rate from active promoter state to inactive promoter state. + b: switching rate from inactive promoter state to active promoter state. + alpha_a: transcription rate for active promoter. + alpha_i: transcription rate for inactive promoter. + gamma: degradation rate. + x0: the initial conditions. + """ # species self.ua = 0 self.ui = 1 @@ -429,6 +533,15 @@ def __init__(self, a=None, b=None, alpha_a=None, alpha_i=None, gamma=None, x0=No self.set_params(a, b, alpha_a, alpha_i, gamma) def ode_func(self, x, t): + """ODE functions to solve. Ignore the splicing part in the base class. + + Args: + x: the variable. + t: the aray of time. + + Returns: + The derivatives dx. + """ dx = np.zeros(len(x)) # parameters a = self.a @@ -447,9 +560,27 @@ def ode_func(self, x, t): return dx def fbar(self, x_a, x_i): + """Calculate the count of a variable by averaging active and inactive states. + + Args: + x_a: the variable x under the active state. + x_i: the variable x under the inactive state. + + Returns: + The count of variable x. + """ return self.b / (self.a + self.b) * x_a + self.a / (self.a + self.b) * x_i def set_params(self, a, b, alpha_a, alpha_i, gamma): + """Set the parameters. + + Args: + a: switching rate from active promoter state to inactive promoter state. + b: switching rate from inactive promoter state to active promoter state. + alpha_a: transcription rate for active promoter. + alpha_i: transcription rate for inactive promoter. + gamma: degradation rate. + """ self.a = a self.b = b self.aa = alpha_a @@ -460,19 +591,39 @@ def set_params(self, a, b, alpha_a, alpha_i, gamma): super().reset() def get_all_central_moments(self): + """Get the first and second central moments for all variables. + + Returns: + An array containing all central moments. + """ ret = np.zeros((2, len(self.t))) ret[0] = self.get_nu() ret[1] = self.get_var_nu() return ret def get_nu(self): + """Get the number of the variable u from the mean averaging active and inactive state. + + Returns: + The number of the variable u. + """ return self.fbar(self.x[:, self.ua], self.x[:, self.ui]) def get_var_nu(self): + """Get the variance of the variable u. + + Returns: + The variance of the variable u. + """ c = self.get_nu() return self.x[:, self.uu] + c - c**2 def computeKnp(self): + """Calculate the K and p from ODE function such that dx = Kx + p. + + Returns: + A tuple containing K and p. + """ # parameters a = self.a b = self.b @@ -502,9 +653,17 @@ def computeKnp(self): class Moments_NoSwitching(LinearODE): + """The class simulates the dynamics of first and second moments of a transcription-splicing system without promoter + switching.""" def __init__(self, alpha=None, beta=None, gamma=None, x0=None): - """This class simulates the dynamics of first and second moments of - a transcription-splicing system without promoter switching.""" + """Initialize the Moments_NoSwitching object. + + Args: + alpha: transcription rate. + beta: splicing rate. + gamma: degradation rate. + x0: the initial conditions. + """ # species self.u = 0 self.s = 1 @@ -522,6 +681,15 @@ def __init__(self, alpha=None, beta=None, gamma=None, x0=None): self.set_params(alpha, beta, gamma) def ode_func(self, x, t): + """ODE functions to solve. Ignore the switching part in the base class. + + Args: + x: the variable. + t: the aray of time. + + Returns: + The derivatives dx. + """ dx = np.zeros(len(x)) # parameters al = self.al @@ -540,6 +708,13 @@ def ode_func(self, x, t): return dx def set_params(self, alpha, beta, gamma): + """Set the parameters. + + Args: + alpha: transcription rate. + beta: splicing rate. + gamma: degradation rate. + """ self.al = alpha self.be = beta self.ga = gamma @@ -548,6 +723,11 @@ def set_params(self, alpha, beta, gamma): super().reset() def get_all_central_moments(self): + """Get the first and second central moments for all variables. + + Returns: + An array containing all central moments. + """ ret = np.zeros((5, len(self.t))) ret[0] = self.get_mean_u() ret[1] = self.get_mean_s() @@ -557,31 +737,66 @@ def get_all_central_moments(self): return ret def get_nosplice_central_moments(self): + """Get the central moments for labeled data. + + Returns: + The central moments. + """ ret = np.zeros((2, len(self.t))) ret[0] = self.get_mean_u() + self.get_mean_s() ret[1] = self.x[:, self.uu] + self.x[:, self.ss] + 2 * self.x[:, self.us] return ret def get_mean_u(self): + """Get the mean of the variable u. + + Returns: + The mean of the variable u. + """ return self.x[:, self.u] def get_mean_s(self): + """Get the mean of the variable s. + + Returns: + The mean of the variable s. + """ return self.x[:, self.s] def get_var_u(self): + """Get the variance of the variable u. + + Returns: + The variance of the variable u. + """ c = self.get_mean_u() return self.x[:, self.uu] - c**2 def get_var_s(self): + """Get the variance of the variable s. + + Returns: + The variance of the variable s. + """ c = self.get_mean_s() return self.x[:, self.ss] - c**2 def get_cov_us(self): + """Get the covariance of the variable s and u. + + Returns: + The covariance. + """ cu = self.get_mean_u() cs = self.get_mean_s() return self.x[:, self.us] - cu * cs def computeKnp(self): + """Calculate the K and p from ODE function such that dx = Kx + p. + + Returns: + A tuple containing K and p. + """ # parameters al = self.al be = self.be @@ -618,9 +833,16 @@ def computeKnp(self): class Moments_NoSwitchingNoSplicing(LinearODE): + """The class simulates the dynamics of first and second moments of a transcription system without promoter + switching.""" def __init__(self, alpha=None, gamma=None, x0=None): - """This class simulates the dynamics of first and second moments of - a transcription system without promoter switching.""" + """Initialize the Moments_NoSwitchingNoSplicing object. + + Args: + alpha: transcription rate. + gamma: degradation rate. + x0: the initial conditions. + """ # species self.u = 0 self.uu = 1 @@ -635,6 +857,15 @@ def __init__(self, alpha=None, gamma=None, x0=None): self.set_params(alpha, gamma) def ode_func(self, x, t): + """ODE functions to solve. Both splicing and switching part in the base class are ignored. + + Args: + x: the variable. + t: the aray of time. + + Returns: + The derivatives dx. + """ dx = np.zeros(len(x)) # parameters al = self.al @@ -649,6 +880,12 @@ def ode_func(self, x, t): return dx def set_params(self, alpha, gamma): + """Set the parameters. + + Args: + alpha: transcription rate. + gamma: degradation rate. + """ self.al = alpha self.ga = gamma @@ -656,19 +893,39 @@ def set_params(self, alpha, gamma): super().reset() def get_all_central_moments(self): + """Get the first and second central moments for all variables. + + Returns: + An array containing all central moments. + """ ret = np.zeros((2, len(self.t))) ret[0] = self.get_mean_u() ret[1] = self.get_var_u() return ret def get_mean_u(self): + """Get the mean of the variable u. + + Returns: + The mean of the variable u. + """ return self.x[:, self.u] def get_var_u(self): + """Get the variance of the variable u. + + Returns: + The variance of the variable u. + """ c = self.get_mean_u() return self.x[:, self.uu] - c**2 def computeKnp(self): + """Calculate the K and p from ODE function such that dx = Kx + p. + + Returns: + A tuple containing K and p. + """ # parameters al = self.al ga = self.ga From c36177f2eeccbee0e162b5bff5208114b641b11a Mon Sep 17 00:00:00 2001 From: sichao Date: Mon, 28 Aug 2023 16:33:54 -0400 Subject: [PATCH 24/36] add typing to the Moments linearODEs --- dynamo/estimation/tsc/utils_kinetic.py | 121 +++++++++++++++---------- 1 file changed, 74 insertions(+), 47 deletions(-) diff --git a/dynamo/estimation/tsc/utils_kinetic.py b/dynamo/estimation/tsc/utils_kinetic.py index 1fb46a9de..046afbcff 100755 --- a/dynamo/estimation/tsc/utils_kinetic.py +++ b/dynamo/estimation/tsc/utils_kinetic.py @@ -254,14 +254,14 @@ class Moments(LinearODE): switching.""" def __init__( self, - a=None, - b=None, - alpha_a=None, - alpha_i=None, - beta=None, - gamma=None, - x0=None, - ): + a: Optional[np.ndarray] = None, + b: Optional[np.ndarray] = None, + alpha_a: Optional[np.ndarray] = None, + alpha_i: Optional[np.ndarray] = None, + beta: Optional[np.ndarray] = None, + gamma: Optional[np.ndarray] = None, + x0: Optional[np.ndarray] = None, + ) -> None: """Initialize the Moments object. Args: @@ -291,7 +291,7 @@ def __init__( if not (a is None or b is None or alpha_a is None or alpha_i is None or beta is None or gamma is None): self.set_params(a, b, alpha_a, alpha_i, beta, gamma) - def ode_func(self, x, t): + def ode_func(self, x: np.ndarray, t: np.ndarray) -> np.ndarray: """ODE functions to solve: dx[u_a] = alpha_a - beta * u_a + a * (u_i - u_a) dx[u_i] = ai - beta * u_i - b * (u_i - u_a) @@ -328,7 +328,7 @@ def ode_func(self, x, t): return dx - def fbar(self, x_a, x_i): + def fbar(self, x_a: np.ndarray, x_i: np.ndarray) -> np.ndarray: """Calculate the count of a variable by averaging active and inactive states. Args: @@ -340,7 +340,15 @@ def fbar(self, x_a, x_i): """ return self.b / (self.a + self.b) * x_a + self.a / (self.a + self.b) * x_i - def set_params(self, a, b, alpha_a, alpha_i, beta, gamma): + def set_params( + self, + a: np.ndarray, + b: np.ndarray, + alpha_a: np.ndarray, + alpha_i: np.ndarray, + beta: np.ndarray, + gamma: np.ndarray, + ) -> None: """Set the parameters. Args: @@ -361,7 +369,7 @@ def set_params(self, a, b, alpha_a, alpha_i, beta, gamma): # reset solutions super().reset() - def get_all_central_moments(self): + def get_all_central_moments(self) -> np.ndarray: """Get the first and second central moments for all variables. Returns: @@ -374,7 +382,7 @@ def get_all_central_moments(self): ret[3] = self.get_var_nx() return ret - def get_nosplice_central_moments(self): + def get_nosplice_central_moments(self) -> np.ndarray: """Get the central moments for labeled data. Returns: @@ -385,7 +393,7 @@ def get_nosplice_central_moments(self): ret[1] = self.get_var_labeled() return ret - def get_nu(self): + def get_nu(self) -> np.ndarray: """Get the number of the variable u from the mean averaging active and inactive state. Returns: @@ -393,7 +401,7 @@ def get_nu(self): """ return self.fbar(self.x[:, self.ua], self.x[:, self.ui]) - def get_nx(self): + def get_nx(self) -> np.ndarray: """Get the number of the variable x from the mean averaging active and inactive state. Returns: @@ -401,7 +409,7 @@ def get_nx(self): """ return self.fbar(self.x[:, self.xa], self.x[:, self.xi]) - def get_n_labeled(self): + def get_n_labeled(self) -> np.ndarray: """Get the number of the labeled data by combining the count of two variables. Returns: @@ -409,7 +417,7 @@ def get_n_labeled(self): """ return self.get_nu() + self.get_nx() - def get_var_nu(self): + def get_var_nu(self) -> np.ndarray: """Get the variance of the variable u. Returns: @@ -418,7 +426,7 @@ def get_var_nu(self): c = self.get_nu() return self.x[:, self.uu] + c - c**2 - def get_var_nx(self): + def get_var_nx(self) -> np.ndarray: """Get the variance of the variable x. Returns: @@ -427,7 +435,7 @@ def get_var_nx(self): c = self.get_nx() return self.x[:, self.xx] + c - c**2 - def get_cov_ux(self): + def get_cov_ux(self) -> np.ndarray: """Get the covariance of the variable x and u. Returns: @@ -437,7 +445,7 @@ def get_cov_ux(self): cx = self.get_nx() return self.x[:, self.ux] - cu * cx - def get_var_labeled(self): + def get_var_labeled(self) -> np.ndarray: """Get the variance of the labeled data by combining the variance of two variables and covariance. Returns: @@ -445,7 +453,7 @@ def get_var_labeled(self): """ return self.get_var_nu() + self.get_var_nx() + 2 * self.get_cov_ux() - def computeKnp(self): + def computeKnp(self) -> Tuple[np.ndarray, np.ndarray]: """Calculate the K and p from ODE function such that dx = Kx + p. Returns: @@ -507,7 +515,15 @@ def computeKnp(self): class Moments_Nosplicing(LinearODE): """The class simulates the dynamics of first and second moments of a transcription-splicing system with promoter switching.""" - def __init__(self, a=None, b=None, alpha_a=None, alpha_i=None, gamma=None, x0=None): + def __init__( + self, + a: Optional[np.ndarray] = None, + b: Optional[np.ndarray] = None, + alpha_a: Optional[np.ndarray] = None, + alpha_i: Optional[np.ndarray] = None, + gamma: Optional[np.ndarray] = None, + x0: Optional[np.ndarray] = None + ) -> None: """Initialize the Moments_Nosplicing object. Args: @@ -532,7 +548,7 @@ def __init__(self, a=None, b=None, alpha_a=None, alpha_i=None, gamma=None, x0=No if not (a is None or b is None or alpha_a is None or alpha_i is None or gamma is None): self.set_params(a, b, alpha_a, alpha_i, gamma) - def ode_func(self, x, t): + def ode_func(self, x: np.ndarray, t: np.ndarray) -> np.ndarray: """ODE functions to solve. Ignore the splicing part in the base class. Args: @@ -559,7 +575,7 @@ def ode_func(self, x, t): return dx - def fbar(self, x_a, x_i): + def fbar(self, x_a: np.ndarray, x_i: np.ndarray) -> np.ndarray: """Calculate the count of a variable by averaging active and inactive states. Args: @@ -571,7 +587,7 @@ def fbar(self, x_a, x_i): """ return self.b / (self.a + self.b) * x_a + self.a / (self.a + self.b) * x_i - def set_params(self, a, b, alpha_a, alpha_i, gamma): + def set_params(self, a: np.ndarray, b: np.ndarray, alpha_a: np.ndarray, alpha_i: np.ndarray, gamma: np.ndarray) -> None: """Set the parameters. Args: @@ -590,7 +606,7 @@ def set_params(self, a, b, alpha_a, alpha_i, gamma): # reset solutions super().reset() - def get_all_central_moments(self): + def get_all_central_moments(self) -> np.ndarray: """Get the first and second central moments for all variables. Returns: @@ -601,7 +617,7 @@ def get_all_central_moments(self): ret[1] = self.get_var_nu() return ret - def get_nu(self): + def get_nu(self) -> np.ndarray: """Get the number of the variable u from the mean averaging active and inactive state. Returns: @@ -609,7 +625,7 @@ def get_nu(self): """ return self.fbar(self.x[:, self.ua], self.x[:, self.ui]) - def get_var_nu(self): + def get_var_nu(self) -> np.ndarray: """Get the variance of the variable u. Returns: @@ -618,7 +634,7 @@ def get_var_nu(self): c = self.get_nu() return self.x[:, self.uu] + c - c**2 - def computeKnp(self): + def computeKnp(self) -> Tuple[np.ndarray, np.ndarray]: """Calculate the K and p from ODE function such that dx = Kx + p. Returns: @@ -655,7 +671,13 @@ def computeKnp(self): class Moments_NoSwitching(LinearODE): """The class simulates the dynamics of first and second moments of a transcription-splicing system without promoter switching.""" - def __init__(self, alpha=None, beta=None, gamma=None, x0=None): + def __init__( + self, + alpha: Optional[np.ndarray] = None, + beta: Optional[np.ndarray] = None, + gamma: Optional[np.ndarray] = None, + x0: Optional[np.ndarray] = None, + ) -> None: """Initialize the Moments_NoSwitching object. Args: @@ -680,7 +702,7 @@ def __init__(self, alpha=None, beta=None, gamma=None, x0=None): if not (alpha is None or beta is None or gamma is None): self.set_params(alpha, beta, gamma) - def ode_func(self, x, t): + def ode_func(self, x: np.ndarray, t: np.ndarray) -> np.ndarray: """ODE functions to solve. Ignore the switching part in the base class. Args: @@ -707,7 +729,7 @@ def ode_func(self, x, t): return dx - def set_params(self, alpha, beta, gamma): + def set_params(self, alpha: np.ndarray, beta: np.ndarray, gamma: np.ndarray) -> None: """Set the parameters. Args: @@ -722,7 +744,7 @@ def set_params(self, alpha, beta, gamma): # reset solutions super().reset() - def get_all_central_moments(self): + def get_all_central_moments(self) -> np.ndarray: """Get the first and second central moments for all variables. Returns: @@ -736,7 +758,7 @@ def get_all_central_moments(self): ret[4] = self.get_var_s() return ret - def get_nosplice_central_moments(self): + def get_nosplice_central_moments(self) -> np.ndarray: """Get the central moments for labeled data. Returns: @@ -747,7 +769,7 @@ def get_nosplice_central_moments(self): ret[1] = self.x[:, self.uu] + self.x[:, self.ss] + 2 * self.x[:, self.us] return ret - def get_mean_u(self): + def get_mean_u(self) -> np.ndarray: """Get the mean of the variable u. Returns: @@ -755,7 +777,7 @@ def get_mean_u(self): """ return self.x[:, self.u] - def get_mean_s(self): + def get_mean_s(self) -> np.ndarray: """Get the mean of the variable s. Returns: @@ -763,7 +785,7 @@ def get_mean_s(self): """ return self.x[:, self.s] - def get_var_u(self): + def get_var_u(self) -> np.ndarray: """Get the variance of the variable u. Returns: @@ -772,7 +794,7 @@ def get_var_u(self): c = self.get_mean_u() return self.x[:, self.uu] - c**2 - def get_var_s(self): + def get_var_s(self) -> np.ndarray: """Get the variance of the variable s. Returns: @@ -781,7 +803,7 @@ def get_var_s(self): c = self.get_mean_s() return self.x[:, self.ss] - c**2 - def get_cov_us(self): + def get_cov_us(self) -> np.ndarray: """Get the covariance of the variable s and u. Returns: @@ -791,7 +813,7 @@ def get_cov_us(self): cs = self.get_mean_s() return self.x[:, self.us] - cu * cs - def computeKnp(self): + def computeKnp(self) -> Tuple[np.ndarray, np.ndarray]: """Calculate the K and p from ODE function such that dx = Kx + p. Returns: @@ -835,7 +857,12 @@ def computeKnp(self): class Moments_NoSwitchingNoSplicing(LinearODE): """The class simulates the dynamics of first and second moments of a transcription system without promoter switching.""" - def __init__(self, alpha=None, gamma=None, x0=None): + def __init__( + self, + alpha: Optional[np.ndarray] = None, + gamma: Optional[np.ndarray] = None, + x0: Optional[np.ndarray] = None, + ) -> None: """Initialize the Moments_NoSwitchingNoSplicing object. Args: @@ -856,7 +883,7 @@ def __init__(self, alpha=None, gamma=None, x0=None): if not (alpha is None or gamma is None): self.set_params(alpha, gamma) - def ode_func(self, x, t): + def ode_func(self, x: np.ndarray, t: np.ndarray) -> np.ndarray: """ODE functions to solve. Both splicing and switching part in the base class are ignored. Args: @@ -879,7 +906,7 @@ def ode_func(self, x, t): return dx - def set_params(self, alpha, gamma): + def set_params(self, alpha: np.ndarray, gamma: np.ndarray) -> None: """Set the parameters. Args: @@ -892,7 +919,7 @@ def set_params(self, alpha, gamma): # reset solutions super().reset() - def get_all_central_moments(self): + def get_all_central_moments(self) -> np.ndarray: """Get the first and second central moments for all variables. Returns: @@ -903,7 +930,7 @@ def get_all_central_moments(self): ret[1] = self.get_var_u() return ret - def get_mean_u(self): + def get_mean_u(self) -> np.ndarray: """Get the mean of the variable u. Returns: @@ -911,7 +938,7 @@ def get_mean_u(self): """ return self.x[:, self.u] - def get_var_u(self): + def get_var_u(self) -> np.ndarray: """Get the variance of the variable u. Returns: @@ -920,7 +947,7 @@ def get_var_u(self): c = self.get_mean_u() return self.x[:, self.uu] - c**2 - def computeKnp(self): + def computeKnp(self) -> Tuple[np.ndarray, np.ndarray]: """Calculate the K and p from ODE function such that dx = Kx + p. Returns: From 10ff2e06050f1fbbddbf55e212e796637c689464 Mon Sep 17 00:00:00 2001 From: sichao Date: Tue, 29 Aug 2023 12:48:43 -0400 Subject: [PATCH 25/36] add missing docstr for utils_velocity --- dynamo/estimation/csc/utils_velocity.py | 49 +++++++++++++++++++++++++ 1 file changed, 49 insertions(+) diff --git a/dynamo/estimation/csc/utils_velocity.py b/dynamo/estimation/csc/utils_velocity.py index f19fbc95d..db99b2256 100755 --- a/dynamo/estimation/csc/utils_velocity.py +++ b/dynamo/estimation/csc/utils_velocity.py @@ -776,11 +776,36 @@ def concat_time_series_matrices(mats, t=None): # --------------------------------------------------------------------------------------------------- # negbin method related def compute_dispersion(mX, varX): + """Compute the reciprocal of the dispersion parameter of NB distribution from the following relationship: + variance = mean + phi * mean^2 + + Args: + mX: the mean of variable X. + varX: the variance of variable X. + + Returns: + The dispersion parameter. + """ phi = fit_linreg(mX**2, varX - mX, intercept=False)[0] return phi def fit_k_negative_binomial(n, r, var, phi=None, k0=None, return_k0=False): + """Fit parameter k for the negative binomial distribution of one gene using the least squares optimizer: + k * r = n + k^2 * variance(r) = k * n - phi * n * 2 + + Args: + n: unspliced mRNA or new/labeled mRNA of one gene. + r: spliced mRNA or total mRNA of one gene. + var: the variance of r. + phi: the dispersion parameter of the negative binomial distribution of one gene. + k0: the initial guess for k. + return_k0: whether to return k0. + + Returns: + The resolution from the least squares optimizer. + """ k0 = fit_linreg(r, n, intercept=False)[0] if k0 is None else k0 phi = compute_dispersion(r, var) if phi is None else phi @@ -795,6 +820,19 @@ def fit_k_negative_binomial(n, r, var, phi=None, k0=None, return_k0=False): def fit_K_negbin(N, R, varR, perc_left=None, perc_right=None, return_phi=False): + """Fit parameter K of negative binomial distribution to each gene. + + Args: + N: unspliced mRNA or new/labeled mRNA. + R: spliced mRNA or total mRNA. + varR: the variance of R. + perc_left: left percentile for selecting extreme data points. + perc_right: right percentile for selecting extreme data points. + return_phi: whether to return dispersion parameter Phi. + + Returns: + Fitted parameter K for each gene. + """ n_gene = N.shape[1] K = np.zeros(n_gene) Phi = np.zeros(n_gene) @@ -815,6 +853,17 @@ def fit_K_negbin(N, R, varR, perc_left=None, perc_right=None, return_phi=False): def compute_velocity_labeling(N, R, K, tau): + """Compute velocity for labeling data by: velocity = gamma / k * new - gamma * total. + + Parameters: + N: new or labeled mRNA. + R: total mRNA. + K: fitted slope k. + tau: time information. + + Returns: + Computed velocity. + """ Kc = np.clip(K, 0, 1 - 1e-3) if np.isscalar(tau): Beta_or_gamma = -np.log(1 - Kc) / tau From e2e8dc0e5aa16a252d29da0ba484940df496ac5f Mon Sep 17 00:00:00 2001 From: sichao Date: Tue, 29 Aug 2023 13:14:14 -0400 Subject: [PATCH 26/36] add typing for fit_k_nb and compute_velocity_labeling --- dynamo/estimation/csc/utils_velocity.py | 29 +++++++++++++++++++++---- 1 file changed, 25 insertions(+), 4 deletions(-) diff --git a/dynamo/estimation/csc/utils_velocity.py b/dynamo/estimation/csc/utils_velocity.py index db99b2256..9ef8d1e25 100755 --- a/dynamo/estimation/csc/utils_velocity.py +++ b/dynamo/estimation/csc/utils_velocity.py @@ -1,3 +1,5 @@ +from typing import Dict, List, Optional, Tuple, Union + import numpy as np import statsmodels.api as sm from scipy.optimize import least_squares @@ -775,7 +777,7 @@ def concat_time_series_matrices(mats, t=None): # --------------------------------------------------------------------------------------------------- # negbin method related -def compute_dispersion(mX, varX): +def compute_dispersion(mX: Union[np.ndarray, csr_matrix], varX: Union[np.ndarray, csr_matrix]) -> float: """Compute the reciprocal of the dispersion parameter of NB distribution from the following relationship: variance = mean + phi * mean^2 @@ -790,7 +792,14 @@ def compute_dispersion(mX, varX): return phi -def fit_k_negative_binomial(n, r, var, phi=None, k0=None, return_k0=False): +def fit_k_negative_binomial( + n: Union[np.ndarray, csr_matrix], + r: Union[np.ndarray, csr_matrix], + var: float, + phi: Optional[float] = None, + k0: Optional[float] = None, + return_k0: bool = False, +) -> Union[Tuple[float, float], float]: """Fit parameter k for the negative binomial distribution of one gene using the least squares optimizer: k * r = n k^2 * variance(r) = k * n - phi * n * 2 @@ -819,7 +828,14 @@ def fit_k_negative_binomial(n, r, var, phi=None, k0=None, return_k0=False): return ret.x[0] -def fit_K_negbin(N, R, varR, perc_left=None, perc_right=None, return_phi=False): +def fit_K_negbin( + N: Union[np.ndarray, csr_matrix], + R: Union[np.ndarray, csr_matrix], + varR: Union[np.ndarray, csr_matrix], + perc_left: Optional[Union[float, int]] = None, + perc_right: Optional[Union[float, int]] = None, + return_phi: bool = False, +) -> Union[Tuple[np.ndarray, np.ndarray], np.ndarray]: """Fit parameter K of negative binomial distribution to each gene. Args: @@ -852,7 +868,12 @@ def fit_K_negbin(N, R, varR, perc_left=None, perc_right=None, return_phi=False): return K -def compute_velocity_labeling(N, R, K, tau): +def compute_velocity_labeling( + N: Union[np.ndarray, csr_matrix], + R: Union[np.ndarray, csr_matrix], + K: Union[np.ndarray, csr_matrix], + tau: Union[np.ndarray, csr_matrix], +) -> Union[np.ndarray, csr_matrix]: """Compute velocity for labeling data by: velocity = gamma / k * new - gamma * total. Parameters: From d236f4abf6a3d962299a5c6ee68f12af22eef6d9 Mon Sep 17 00:00:00 2001 From: sichao Date: Tue, 29 Aug 2023 18:40:36 -0400 Subject: [PATCH 27/36] add typing for sol funcs in utils_velocity --- dynamo/estimation/csc/utils_velocity.py | 30 +++++++++++++++++++------ 1 file changed, 23 insertions(+), 7 deletions(-) diff --git a/dynamo/estimation/csc/utils_velocity.py b/dynamo/estimation/csc/utils_velocity.py index 9ef8d1e25..564e93ced 100755 --- a/dynamo/estimation/csc/utils_velocity.py +++ b/dynamo/estimation/csc/utils_velocity.py @@ -10,7 +10,7 @@ from ...tools.utils import elem_prod, find_extreme -def sol_u(t, u0, alpha, beta): +def sol_u(t: np.ndarray, u0: float, alpha: float, beta: float) -> np.ndarray: """The analytical solution of unspliced mRNA kinetics. Arguments @@ -32,7 +32,7 @@ def sol_u(t, u0, alpha, beta): return u0 * np.exp(-beta * t) + alpha / beta * (1 - np.exp(-beta * t)) -def sol_u_2p(t, u0, t1, alpha0, alpha1, beta): +def sol_u_2p(t: np.ndarray, u0: float, t1: np.ndarray, alpha0: float, alpha1: float, beta: float) -> np.ndarray: """The combined 2-piece analytical solution of unspliced mRNA kinetics. Arguments @@ -62,7 +62,7 @@ def sol_u_2p(t, u0, t1, alpha0, alpha1, beta): return np.concatenate((u_pre, u_aft)) -def sol_s(t, s0, u0, alpha, beta, gamma): +def sol_s(t: np.ndarray, s0: float, u0: float, alpha: float, beta: float, gamma: float) -> np.ndarray: """The analytical solution of spliced mRNA kinetics. Arguments @@ -97,7 +97,17 @@ def sol_s(t, s0, u0, alpha, beta, gamma): return s -def sol_p(t, p0, s0, u0, alpha, beta, gamma, eta, delta): +def sol_p( + t: np.ndarray, + p0: float, + s0: float, + u0: float, + alpha: float, + beta: float, + gamma: float, + eta: float, + delta: float, +) -> Tuple[np.ndarray, np.ndarray, np.ndarray]: """The analytical solution of protein kinetics. Arguments @@ -139,7 +149,7 @@ def sol_p(t, p0, s0, u0, alpha, beta, gamma, eta, delta): return p, s, u -def solve_gamma(t, old, total): +def solve_gamma(t: float, old: np.ndarray, total: np.ndarray) -> np.ndarray: """Analytical solution to calculate gamma (degradation rate) using first-order degradation kinetics. Parameters @@ -162,7 +172,7 @@ def solve_gamma(t, old, total): return gamma -def solve_alpha_2p(t0, t1, alpha0, beta, u1): +def solve_alpha_2p(t0: float, t1: float, alpha0: float, beta: float, u1: Union[csr_matrix, np.ndarray]) -> float: """Given known steady state alpha and beta, solve stimulation alpha for a mixed steady state and stimulation labeling experiment. Parameters @@ -195,7 +205,13 @@ def solve_alpha_2p(t0, t1, alpha0, beta, u1): return alpha1 -def solve_alpha_2p_mat(t0, t1, alpha0, beta, u1): +def solve_alpha_2p_mat( + t0: np.ndarray, + t1: np.ndarray, + alpha0: np.ndarray, + beta: np.ndarray, + u1: np.ndarray, +) -> Tuple[csr_matrix, csr_matrix, csr_matrix]: """Given known steady state alpha and beta, solve stimulation alpha for a mixed steady state and stimulation labeling experiment in a matrix form. From 1eeb70ddf804a0bd4681a1027508c0b366270494 Mon Sep 17 00:00:00 2001 From: sichao Date: Wed, 30 Aug 2023 10:48:05 -0400 Subject: [PATCH 28/36] add docstr for fit funcs in utils_velocity --- dynamo/estimation/csc/utils_velocity.py | 78 +++++++++++++++++++++---- 1 file changed, 67 insertions(+), 11 deletions(-) diff --git a/dynamo/estimation/csc/utils_velocity.py b/dynamo/estimation/csc/utils_velocity.py index 564e93ced..049fdee82 100755 --- a/dynamo/estimation/csc/utils_velocity.py +++ b/dynamo/estimation/csc/utils_velocity.py @@ -250,7 +250,13 @@ def solve_alpha_2p_mat( return csr_matrix(u0), csr_matrix(u_new), csr_matrix(alpha1) -def fit_linreg(x, y, mask=None, intercept=False, r2=True): +def fit_linreg( + x: Union[csr_matrix, np.ndarray], + y: Union[csr_matrix, np.ndarray], + mask: Optional[Union[List, np.ndarray]] = None, + intercept: bool = False, + r2: bool = True, +) -> Union[Tuple[float, float, float, float], Tuple[float, float]]: """Simple linear regression: y = kx + b. Arguments @@ -315,7 +321,14 @@ def fit_linreg(x, y, mask=None, intercept=False, r2=True): return k, b -def fit_linreg_robust(x, y, mask=None, intercept=False, r2=True, est_method="rlm"): +def fit_linreg_robust( + x: Union[csr_matrix, np.ndarray], + y: Union[csr_matrix, np.ndarray], + mask: Optional[Union[List, np.ndarray]] = None, + intercept: bool = False, + r2: bool = True, + est_method: str = "rlm", +) -> Union[Tuple[float, float, float, float], Tuple[float, float]]: """Apply robust linear regression of y w.r.t x. Arguments @@ -398,7 +411,14 @@ def fit_linreg_robust(x, y, mask=None, intercept=False, r2=True, est_method="rlm return k, b -def fit_stochastic_linreg(u, s, us, ss, fit_2_gammas=True, err_cov=False): +def fit_stochastic_linreg( + u: Union[csr_matrix, np.ndarray], + s: Union[csr_matrix, np.ndarray], + us: Union[csr_matrix, np.ndarray], + ss: Union[csr_matrix, np.ndarray], + fit_2_gammas: bool = True, + err_cov: bool = False, +) -> float: """generalized method of moments: [u, 2*us + u] = gamma * [s, 2*ss - s]. Arguments @@ -449,7 +469,13 @@ def fit_stochastic_linreg(u, s, us, ss, fit_2_gammas=True, err_cov=False): return gamma -def fit_first_order_deg_lsq(t, l, bounds=(0, np.inf), fix_l0=False, beta_0=1): +def fit_first_order_deg_lsq( + t: np.ndarray, + l: Union[csr_matrix, np.ndarray], + bounds: Tuple = (0, np.inf), + fix_l0: bool = False, + beta_0: Union[float, int] = 1, +) -> Tuple[float, float]: """Estimate beta with degradation data using least squares method. Arguments @@ -490,7 +516,7 @@ def fit_first_order_deg_lsq(t, l, bounds=(0, np.inf), fix_l0=False, beta_0=1): return beta, l0 -def solve_first_order_deg(t, l): +def solve_first_order_deg(t: np.ndarray, l: Union[csr_matrix, np.ndarray]) -> Tuple[float, float, float]: """Solve for the initial amount and the rate constant of a species (for example, labeled mRNA) with time-series data under first-order degration kinetics model. @@ -523,7 +549,14 @@ def solve_first_order_deg(t, l): return l0, k, half_life -def fit_gamma_lsq(t, s, beta, u0, bounds=(0, np.inf), fix_s0=False): +def fit_gamma_lsq( + t: np.ndarray, + s: Union[csr_matrix, np.ndarray], + beta: float, + u0: float, + bounds: Tuple = (0, np.inf), + fix_s0: bool = False, +) -> Tuple[float, float]: """Estimate gamma with degradation data using least squares method. Arguments @@ -570,7 +603,7 @@ def fit_gamma_lsq(t, s, beta, u0, bounds=(0, np.inf), fix_s0=False): return gamma, s0 -def fit_alpha_synthesis(t, u, beta): +def fit_alpha_synthesis(t: np.ndarray, u: Union[csr_matrix, np.ndarray], beta: float) -> float: """Estimate alpha with synthesis data using linear regression with fixed zero intercept. It is assumed that u(0) = 0. @@ -599,7 +632,12 @@ def fit_alpha_synthesis(t, u, beta): return beta * np.mean(u) / np.mean(x) -def fit_alpha_degradation(t, u, beta, intercept=False): +def fit_alpha_degradation( + t: np.ndarray, + u: Union[csr_matrix, np.ndarray], + beta: float, + intercept: bool = False, +) -> Union[float, float, float]: """Estimate alpha with degradation data using linear regression. This is a lsq version of the following function that constrains u0 to be larger than 0 @@ -641,7 +679,12 @@ def fit_alpha_degradation(t, u, beta, intercept=False): return alpha, u0, r2 -def solve_alpha_degradation(t, u, beta, intercept=False): +def solve_alpha_degradation( + t: np.ndarray, + u: Union[csr_matrix, np.ndarray], + beta: float, + intercept: bool = False, +) -> Union[float, float, float]: """Estimate alpha with degradation data using linear regression. Arguments @@ -694,7 +737,13 @@ def solve_alpha_degradation(t, u, beta, intercept=False): return k * beta, b, r2 -def fit_alpha_beta_synthesis(t, l, bounds=(0, np.inf), alpha_0=1, beta_0=1): +def fit_alpha_beta_synthesis( + t: np.ndarray, + l: Union[csr_matrix, np.ndarray], + bounds: Tuple = (0, np.inf), + alpha_0: Union[float, int] = 1, + beta_0: Union[float, int] = 1, +) -> Tuple[float, float]: """Estimate alpha and beta with synthesis data using least square method. It is assumed that u(0) = 0. @@ -728,7 +777,14 @@ def fit_alpha_beta_synthesis(t, l, bounds=(0, np.inf), alpha_0=1, beta_0=1): return ret.x[0], ret.x[1] -def fit_all_synthesis(t, l, bounds=(0, np.inf), alpha_0=1, beta_0=1, gamma_0=1): +def fit_all_synthesis( + t: np.ndarray, + l: Union[csr_matrix, np.ndarray], + bounds: Tuple = (0, np.inf), + alpha_0: Union[float, int] = 1, + beta_0: Union[float, int] = 1, + gamma_0: Union[float, int] = 1, +) -> Tuple[float, float, float]: """Estimate alpha, beta and gamma with synthesis data using least square method. It is assumed that u(0) = 0 and s(0) = 0. From 19ce570675ebdfd78bf229b252aaf46b3c5ed884 Mon Sep 17 00:00:00 2001 From: sichao Date: Wed, 30 Aug 2023 10:58:20 -0400 Subject: [PATCH 29/36] add docstr for concat_time_series_matrices and compute_bursting_properties --- dynamo/estimation/csc/utils_velocity.py | 11 +++++++++-- 1 file changed, 9 insertions(+), 2 deletions(-) diff --git a/dynamo/estimation/csc/utils_velocity.py b/dynamo/estimation/csc/utils_velocity.py index 049fdee82..09c58079c 100755 --- a/dynamo/estimation/csc/utils_velocity.py +++ b/dynamo/estimation/csc/utils_velocity.py @@ -822,7 +822,10 @@ def fit_all_synthesis( return ret.x[0], ret.x[1], ret.x[2] -def concat_time_series_matrices(mats, t=None): +def concat_time_series_matrices( + mats: np.ndarray, + t: Optional[Union[List, np.ndarray]] = None, +) -> Union[Tuple[np.ndarray, np.ndarray], np.ndarray]: """Concatenate a list of gene x cell matrices into a single matrix. Arguments @@ -965,7 +968,11 @@ def compute_velocity_labeling( return elem_prod(Beta_or_gamma, N) / K - elem_prod(Beta_or_gamma, R) -def compute_bursting_properties(M_t, phi, gamma): +def compute_bursting_properties( + M_t: Union[float, np.ndarray], + phi: float, + gamma: float, +) -> Tuple[Union[float, np.ndarray], np.ndarray]: """Compute bursting frequency and size for the negative binomial regression model. The equations come from: Anton J.M. Larsson et al. Genomic encoding of transcriptional burst kinetics, Nature volume 565, pages251–254(2019) From de740a9924d102d3034e03f3bd5a3b938be99348 Mon Sep 17 00:00:00 2001 From: sichao Date: Wed, 30 Aug 2023 11:07:46 -0400 Subject: [PATCH 30/36] update docstr for negbin funcs --- dynamo/estimation/csc/utils_velocity.py | 30 +++++++++---------------- 1 file changed, 11 insertions(+), 19 deletions(-) diff --git a/dynamo/estimation/csc/utils_velocity.py b/dynamo/estimation/csc/utils_velocity.py index 09c58079c..0dfadee40 100755 --- a/dynamo/estimation/csc/utils_velocity.py +++ b/dynamo/estimation/csc/utils_velocity.py @@ -951,7 +951,7 @@ def compute_velocity_labeling( ) -> Union[np.ndarray, csr_matrix]: """Compute velocity for labeling data by: velocity = gamma / k * new - gamma * total. - Parameters: + Args: N: new or labeled mRNA. R: total mRNA. K: fitted slope k. @@ -973,26 +973,18 @@ def compute_bursting_properties( phi: float, gamma: float, ) -> Tuple[Union[float, np.ndarray], np.ndarray]: - """Compute bursting frequency and size for the negative binomial regression model. - The equations come from: - Anton J.M. Larsson et al. Genomic encoding of transcriptional burst kinetics, Nature volume 565, pages251–254(2019) + """Compute bursting frequency and size for the negative binomial regression model. The equations come from: + Anton J.M. Larsson et al. Genomic encoding of transcriptional burst kinetics, Nature volume 565, pages251–254(2019) - Arguments - --------- - M_t: :class:`~numpy.ndarray` or float - The first moment of the number of total mRNA. - If an array is passed, a cell-wise bursting size will be calculated. - phi: float - The reciprocal dispersion parameter. - gamma: float - Degradation rate constant. + Args: + M_t: the first moment of the number of total mRNA. If an array is passed, a cell-wise bursting size will be + calculated. + phi: the reciprocal dispersion parameter. + gamma: degradation rate constant. - Returns - ------- - bs: :class:`~numpy.ndarray` or float - Bursting size - bf: float - Bursting frequency + Returns: + bs: bursting size + bf: bursting frequency """ bs = M_t * phi bf = gamma / phi From 77421cdae5a3e1943cd597e47eaf6fbb90333d1f Mon Sep 17 00:00:00 2001 From: sichao Date: Wed, 30 Aug 2023 12:09:13 -0400 Subject: [PATCH 31/36] update docstr style for utils_velocity --- dynamo/estimation/csc/utils_velocity.py | 544 ++++++++---------------- 1 file changed, 187 insertions(+), 357 deletions(-) diff --git a/dynamo/estimation/csc/utils_velocity.py b/dynamo/estimation/csc/utils_velocity.py index 0dfadee40..dceb0b3d8 100755 --- a/dynamo/estimation/csc/utils_velocity.py +++ b/dynamo/estimation/csc/utils_velocity.py @@ -13,20 +13,13 @@ def sol_u(t: np.ndarray, u0: float, alpha: float, beta: float) -> np.ndarray: """The analytical solution of unspliced mRNA kinetics. - Arguments - --------- - t: :class:`~numpy.ndarray` - A vector of time points. - u0: float - Initial value of u. - alpha: float - Transcription rate. - beta: float - Splicing rate constant. - - Returns - ------- - u: :class:`~numpy.ndarray` + Args: + t: a vector of time points. + u0: initial value of u. + alpha: transcription rate. + beta: splicing rate constant. + + Returns: Unspliced mRNA counts at given time points. """ return u0 * np.exp(-beta * t) + alpha / beta * (1 - np.exp(-beta * t)) @@ -35,24 +28,15 @@ def sol_u(t: np.ndarray, u0: float, alpha: float, beta: float) -> np.ndarray: def sol_u_2p(t: np.ndarray, u0: float, t1: np.ndarray, alpha0: float, alpha1: float, beta: float) -> np.ndarray: """The combined 2-piece analytical solution of unspliced mRNA kinetics. - Arguments - --------- - t: :class:`~numpy.ndarray` - A vector of time points for both steady state and stimulation labeling. - u0: float - Initial value of u. - t1: :class:`~numpy.ndarray` - The time point when the cells switch from steady state to stimulation. - alpha0: float - Transcription rate for steady state labeling. - alpha1: float - Transcription rate for stimulation based labeling. - beta: float - Splicing rate constant. - - Returns - ------- - u: :class:`~numpy.ndarray` + Args: + t: a vector of time points for both steady state and stimulation labeling. + u0: initial value of u. + t1: the time point when the cells switch from steady state to stimulation. + alpha0: transcription rate for steady state labeling. + alpha1: transcription rate for stimulation based labeling. + beta: splicing rate constant. + + Returns: Unspliced mRNA counts at given time points. """ u1 = sol_u(t1, u0, alpha0, beta) @@ -65,24 +49,15 @@ def sol_u_2p(t: np.ndarray, u0: float, t1: np.ndarray, alpha0: float, alpha1: fl def sol_s(t: np.ndarray, s0: float, u0: float, alpha: float, beta: float, gamma: float) -> np.ndarray: """The analytical solution of spliced mRNA kinetics. - Arguments - --------- - t: :class:`~numpy.ndarray` - A vector of time points. - s0: float - Initial value of s. - u0: float - Initial value of u. - alpha: float - Transcription rate. - beta: float - Splicing rate constant. - gamma: float - Degradation rate constant for spliced mRNA. - - Returns - ------- - s: :class:`~numpy.ndarray` + Args: + t: a vector of time points. + s0: initial value of s. + u0: initial value of u. + alpha: transcription rate. + beta: splicing rate constant. + gamma: degradation rate constant for spliced mRNA. + + Returns: Spliced mRNA counts at given time points. """ exp_gt = np.exp(-gamma * t) @@ -110,35 +85,20 @@ def sol_p( ) -> Tuple[np.ndarray, np.ndarray, np.ndarray]: """The analytical solution of protein kinetics. - Arguments - --------- - t: :class:`~numpy.ndarray` - A vector of time points. - p0: float - Initial value of p. - s0: float - Initial value of s. - u0: float - Initial value of u. - alpha: float - Transcription rate. - beta: float - Splicing rate constant. - gamma: float - Degradation rate constant for spliced mRNA. - eta: float - Synthesis rate constant for protein. - delta: float - Degradation rate constant for protein. - - Returns - ------- - p: :class:`~numpy.ndarray` - Protein counts at given time points. - s: :class:`~numpy.ndarray` - Spliced mRNA counts at given time points. - u: :class:`~numpy.ndarray` - Unspliced mRNA counts at given time points. + Args: + t: a vector of time points. + p0: initial value of p. + s0: initial value of s. + u0: initial value of u. + alpha: transcription rate. + beta: splicing rate constant. + gamma: degradation rate constant for spliced mRNA. + eta: synthesis rate constant for protein. + delta: degradation rate constant for protein. + + Returns: + Protein counts at given time points (p), spliced mRNA counts at given time points (s) and unspliced mRNA counts + at given time points (u). """ u = sol_u(t, u0, alpha, beta) s = sol_s(t, s0, u0, alpha, beta, gamma) @@ -152,18 +112,13 @@ def sol_p( def solve_gamma(t: float, old: np.ndarray, total: np.ndarray) -> np.ndarray: """Analytical solution to calculate gamma (degradation rate) using first-order degradation kinetics. - Parameters - ---------- - t: `float` - Metabolic labeling time period. - old: :class:`~numpy.ndarray` or sparse `csr_matrix` - A vector of old RNA amount in each cell - total: :class:`~numpy.ndarray` or sparse `csr_matrix` - A vector of total RNA amount in each cell - - Returns - ------- - Returns the degradation rate (gamma) + Args: + t: metabolic labeling time period. + old: a vector of old RNA amount in each cell + total: a vector of total RNA amount in each cell + + Returns: + The degradation rate (gamma). """ old, total = np.mean(old), np.mean(total) @@ -173,24 +128,18 @@ def solve_gamma(t: float, old: np.ndarray, total: np.ndarray) -> np.ndarray: def solve_alpha_2p(t0: float, t1: float, alpha0: float, beta: float, u1: Union[csr_matrix, np.ndarray]) -> float: - """Given known steady state alpha and beta, solve stimulation alpha for a mixed steady state and stimulation labeling experiment. - - Parameters - ---------- - t0: `float` - Time period for steady state labeling. - t1: `float` - Time period for stimulation labeling. - alpha0: `float` - steady state transcription rate calculated from one-shot experiment mode. - beta: `float` - steady state (and simulation) splicing rate calculated from one-shot experiment mode. - u1: :class:`~numpy.ndarray` or sparse `csr_matrix` - A vector of labeled RNA amount in each cell observed at time t0 + t1. - - Returns - ------- - Returns the transcription rate (alpha1) for the stimulation period in the data. + """Given known steady state alpha and beta, solve stimulation alpha for a mixed steady state and stimulation + labeling experiment. + + Args: + t0: time period for steady state labeling. + t1: time period for stimulation labeling. + alpha0: steady state transcription rate calculated from one-shot experiment mode. + beta: steady state (and simulation) splicing rate calculated from one-shot experiment mode. + u1: a vector of labeled RNA amount in each cell observed at time t0 + t1. + + Returns: + The transcription rate (alpha1) for the stimulation period in the data. """ u1 = np.mean(u1) @@ -215,22 +164,15 @@ def solve_alpha_2p_mat( """Given known steady state alpha and beta, solve stimulation alpha for a mixed steady state and stimulation labeling experiment in a matrix form. - Parameters - ---------- - t0: :class:`~numpy.ndarray` - Time period for steady state labeling for each cell. - t1: :class:`~numpy.ndarray` - Time period for stimulation labeling for each cell. - alpha0: :class:`~numpy.ndarray` - steady state transcription rate calculated from one-shot experiment mode for each gene. - beta: :class:`~numpy.ndarray` - steady state (and simulation) splicing rate calculated from one-shot experiment mode for each gene. - u1: :class:`~numpy.ndarray` or sparse `csr_matrix` - A vector of labeled RNA amount in each cell observed at time t0 + t1. - - Returns - ------- - Returns the transcription rate (alpha1) for the stimulation period in the data. + Args: + t0: time period for steady state labeling for each cell. + t1: time period for stimulation labeling for each cell. + alpha0: steady state transcription rate calculated from one-shot experiment mode for each gene. + beta: steady state (and simulation) splicing rate calculated from one-shot experiment mode for each gene. + u1: a vector of labeled RNA amount in each cell observed at time t0 + t1. + + Returns: + The transcription rate (alpha1) for the stimulation period in the data. """ alpha0 = np.repeat(alpha0.reshape((-1, 1)), u1.shape[1], axis=1) @@ -259,27 +201,19 @@ def fit_linreg( ) -> Union[Tuple[float, float, float, float], Tuple[float, float]]: """Simple linear regression: y = kx + b. - Arguments - --------- - x: :class:`~numpy.ndarray` or sparse `csr_matrix` - A vector of independent variables. - y: :class:`~numpy.ndarray` or sparse `csr_matrix` - A vector of dependent variables. - intercept: bool - If using steady state assumption for fitting, then: - True -- the linear regression is performed with an unfixed intercept; - False -- the linear regresssion is performed with a fixed zero intercept - - Returns - ------- - k: float - The estimated slope. - b: float - The estimated intercept. - r2: float - Coefficient of determination or r square calculated with the extreme data points. - all_r2: float - The r2 calculated using all data points. + Args: + x: a vector of independent variables. + y: a vector of dependent variables. + mask: the index to filter the data before regression. + intercept: if using steady state assumption for fitting, then: + True -- the linear regression is performed with an unfixed intercept; + False -- the linear regression is performed with a fixed zero intercept. + r2: the R squared value calculated from the regression. + + Returns: + The estimated slope (k) , the estimated intercept (b), coefficient of determination or r square calculated with + the extreme data points (r2), the r2 calculated using all data points (all_r2). If argument r2 is False, r2 and + all_r2 will not be returned. """ x = x.A if issparse(x) else x y = y.A if issparse(y) else y @@ -331,29 +265,18 @@ def fit_linreg_robust( ) -> Union[Tuple[float, float, float, float], Tuple[float, float]]: """Apply robust linear regression of y w.r.t x. - Arguments - --------- - x: :class:`~numpy.ndarray` or sparse `csr_matrix` - A vector of independent variables. - y: :class:`~numpy.ndarray` or sparse `csr_matrix` - A vector of dependent variables. - intercept: bool - If using steady state assumption for fitting, then: - True -- the linear regression is performed with an unfixed intercept; - False -- the linear regresssion is performed with a fixed zero intercept. - est_method: str (default: `rlm`) - The linear regression estimation method that will be used. - - Returns - ------- - k: float - The estimated slope. - b: float - The estimated intercept. - r2: float - Coefficient of determination or r square calculated with the extreme data points. - all_r2: float - The r2 calculated using all data points. + Args: + x: a vector of independent variables. + y: a vector of dependent variables. + intercept: if using steady state assumption for fitting, then: + True -- the linear regression is performed with an unfixed intercept; + False -- the linear regression is performed with a fixed zero intercept. + est_method: the linear regression estimation method that will be used. + + Returns: + The estimated slope (k), the estimated intercept (b), coefficient of determination or r square calculated with + the extreme data points (r2), the r2 calculated using all data points (all_r2).If argument r2 is False, r2 and + all_r2 will not be returned. """ x = x.A if issparse(x) else x @@ -419,22 +342,15 @@ def fit_stochastic_linreg( fit_2_gammas: bool = True, err_cov: bool = False, ) -> float: - """generalized method of moments: [u, 2*us + u] = gamma * [s, 2*ss - s]. - - Arguments - --------- - u: :class:`~numpy.ndarray` or sparse `csr_matrix` - A vector of first moments (mean) of unspliced (or new) RNA expression. - s: :class:`~numpy.ndarray` or sparse `csr_matrix` - A vector of first moments (mean) of spliced (or total) RNA expression. - us: :class:`~numpy.ndarray` or sparse `csr_matrix` - A vector of second moments (uncentered co-variance) of unspliced/spliced (or new/total) RNA expression. - ss: :class:`~numpy.ndarray` or sparse `csr_matrix` - A vector of second moments (uncentered variance) of spliced (or total) RNA expression. - - Returns - ------- - gamma: float + """Generalized method of moments: [u, 2*us + u] = gamma * [s, 2*ss - s]. + + Args: + u: a vector of first moments (mean) of unspliced (or new) RNA expression. + s: a vector of first moments (mean) of spliced (or total) RNA expression. + us: a vector of second moments (uncentered co-variance) of unspliced/spliced (or new/total) RNA expression. + ss: a vector of second moments (uncentered variance) of spliced (or total) RNA expression. + + Returns: The estimated gamma. """ y = np.vstack((u.flatten(), (u + 2 * us).flatten())) @@ -478,26 +394,16 @@ def fit_first_order_deg_lsq( ) -> Tuple[float, float]: """Estimate beta with degradation data using least squares method. - Arguments - --------- - t: :class:`~numpy.ndarray` - A vector of time points. - l: :class:`~numpy.ndarray` or sparse `csr_matrix` - A vector of unspliced, labeled mRNA counts for each time point. - bounds: tuple - The bound for beta. The default is beta > 0. - fixed_l0: bool - True: l0 will be calculated by averaging the first column of l; - False: l0 is a parameter that will be estimated all together with beta using lsq. - beta_0: float - Initial guess for beta. - - Returns - ------- - beta: float - The estimated value for beta. - l0: float - The estimated value for the initial spliced, labeled mRNA count. + Args: + t: a vector of time points. + l: a vector of unspliced, labeled mRNA counts for each time point. + bounds: the bound for beta. The default is beta > 0. + fix_l0: whether to used fixed l0. If ture, l0 will be calculated by averaging the first column of l; + If false, l0 is a parameter that will be estimated all together with beta using lsq. + beta_0: initial guess for beta. + + Returns: + The estimated value for beta (beta) and the estimated value for the initial spliced, labeled mRNA count (l0). """ l = l.A.flatten() if issparse(l) else l @@ -518,23 +424,14 @@ def fit_first_order_deg_lsq( def solve_first_order_deg(t: np.ndarray, l: Union[csr_matrix, np.ndarray]) -> Tuple[float, float, float]: """Solve for the initial amount and the rate constant of a species (for example, labeled mRNA) with time-series data - under first-order degration kinetics model. - - Parameters - ---------- - t: :class:`~numpy.ndarray` - A vector of time points. - l: :class:`~numpy.ndarray` or sparse `csr_matrix` - A vector of labeled mRNA counts for each time point. - - Returns - ------- - l0:: `float` - The intial counts of the species (for example, labeled mRNA). - k: `float` - Degradation rate constant. - half_life: `float` - Half-life the species. + under first-order degradation kinetics model. + + Args: + t: a vector of time points. + l: a vector of labeled mRNA counts for each time point. + + Returns: + The initial counts of the species (for example, labeled mRNA), degradation rate constant and half-life the species. """ x = l.A.flatten() if issparse(l) else l @@ -559,28 +456,17 @@ def fit_gamma_lsq( ) -> Tuple[float, float]: """Estimate gamma with degradation data using least squares method. - Arguments - --------- - t: :class:`~numpy.ndarray` - A vector of time points. - s: :class:`~numpy.ndarray` or sparse `csr_matrix` - A vector of spliced, labeled mRNA counts for each time point. - beta: float - The value of beta. - u0: float - Initial number of unspliced mRNA. - bounds: tuple - The bound for gamma. The default is gamma > 0. - fixed_s0: bool - True: s0 will be calculated by averaging the first column of s; - False: s0 is a parameter that will be estimated all together with gamma using lsq. - - Returns - ------- - gamma: float - The estimated value for gamma. - s0: float - The estimated value for the initial spliced mRNA count. + Args: + t: a vector of time points. + s: a vector of spliced, labeled mRNA counts for each time point. + beta: the value of beta. + u0: initial number of unspliced mRNA. + bounds: the bound for gamma. The default is gamma > 0. + fix_s0: whether to use fixed s0. If true, s0 will be calculated by averaging the first column of s; + If false, s0 is a parameter that will be estimated all together with gamma using lsq. + + Returns: + The estimated value for gamma and the estimated value for the initial spliced mRNA count. """ s = s.A.flatten() if issparse(s) else s @@ -607,18 +493,12 @@ def fit_alpha_synthesis(t: np.ndarray, u: Union[csr_matrix, np.ndarray], beta: f """Estimate alpha with synthesis data using linear regression with fixed zero intercept. It is assumed that u(0) = 0. - Arguments - --------- - t: :class:`~numpy.ndarray` - A vector of time points. - u: :class:`~numpy.ndarray` or sparse `csr_matrix` - A matrix of unspliced mRNA counts. Dimension: cells x time points. - beta: float - The value of beta. - - Returns - ------- - alpha: float + Args: + t: a vector of time points. + u: a matrix of unspliced mRNA counts. Dimension: cells x time points. + beta: the value of beta. + + Returns: The estimated value for alpha. """ u = u.A if issparse(u) else u @@ -639,29 +519,19 @@ def fit_alpha_degradation( intercept: bool = False, ) -> Union[float, float, float]: """Estimate alpha with degradation data using linear regression. This is a lsq version of the following function that - constrains u0 to be larger than 0 - - Arguments - --------- - t: :class:`~numpy.ndarray` - A vector of time points. - u: :class:`~numpy.ndarray` or sparse `csr_matrix` - A matrix of unspliced mRNA counts. Dimension: cells x time points. - beta: float - The value of beta. - intercept: bool - If using steady state assumption for fitting, then: - True -- the linear regression is performed with an unfixed intercept; - False -- the linear regresssion is performed with a fixed zero intercept. - - Returns - ------- - alpha: float - The estimated value for alpha. - u0: float - The initial unspliced mRNA count. - r2: float - Coefficient of determination or r square. + constrains u0 to be larger than 0. + + Args: + t: a vector of time points. + u: a matrix of unspliced mRNA counts. Dimension: cells x time points. + beta: the value of beta. + intercept: whether to steady state assumption for fitting, then: + True -- the linear regression is performed with an unfixed intercept; + False -- the linear regression is performed with a fixed zero intercept. + + Returns: + The estimated value for alpha (alpha), the initial unspliced mRNA count (u0), and coefficient of determination + or r square (r2). """ x = u.A if issparse(u) else u @@ -687,27 +557,17 @@ def solve_alpha_degradation( ) -> Union[float, float, float]: """Estimate alpha with degradation data using linear regression. - Arguments - --------- - t: :class:`~numpy.ndarray` - A vector of time points. - u: :class:`~numpy.ndarray` or sparse `csr_matrix` - A matrix of unspliced mRNA counts. Dimension: cells x time points. - beta: float - The value of beta. - intercept: bool - If using steady state assumption for fitting, then: - True -- the linear regression is performed with an unfixed intercept; - False -- the linear regresssion is performed with a fixed zero intercept. - - Returns - ------- - alpha: float - The estimated value for alpha. - b: float - The initial unspliced mRNA count. - r2: float - Coefficient of determination or r square. + Args: + t: a vector of time points. + u: a matrix of unspliced mRNA counts. Dimension: cells x time points. + beta: the value of beta. + intercept: if using steady state assumption for fitting, then: + True -- the linear regression is performed with an unfixed intercept; + False -- the linear regression is performed with a fixed zero intercept. + + Returns: + The estimated value for alpha (alpha), the initial unspliced mRNA count (b), coefficient of determination or r + square. """ u = u.A if issparse(u) else u @@ -744,28 +604,18 @@ def fit_alpha_beta_synthesis( alpha_0: Union[float, int] = 1, beta_0: Union[float, int] = 1, ) -> Tuple[float, float]: - """Estimate alpha and beta with synthesis data using least square method. + """Estimate alpha and beta with synthesis data using the least square method. It is assumed that u(0) = 0. - Arguments - --------- - t: :class:`~numpy.ndarray` - A vector of time points. - l: :class:`~numpy.ndarray` or sparse `csr_matrix` - A matrix of labeled mRNA counts. Dimension: cells x time points. - bounds: tuple - The bound for alpha and beta. The default is alpha / beta > 0. - alpha_0: float - Initial guess for alpha. - beta_0: float - Initial guess for beta. - - Returns - ------- - alpha: float - The estimated value for alpha. - beta: float - The estimated value for beta. + Args: + t: a vector of time points. + l: a matrix of labeled mRNA counts. Dimension: cells x time points. + bounds: the bound for alpha and beta. The default is alpha / beta > 0. + alpha_0: initial guess for alpha. + beta_0: initial guess for beta. + + Returns: + The estimated value for alpha and the estimated value for beta. """ l = l.A if issparse(l) else l @@ -785,32 +635,19 @@ def fit_all_synthesis( beta_0: Union[float, int] = 1, gamma_0: Union[float, int] = 1, ) -> Tuple[float, float, float]: - """Estimate alpha, beta and gamma with synthesis data using least square method. + """Estimate alpha, beta and gamma with synthesis data using the least square method. It is assumed that u(0) = 0 and s(0) = 0. - Arguments - --------- - t: :class:`~numpy.ndarray` - A vector of time points. - l: :class:`~numpy.ndarray` or sparse `csr_matrix` - A matrix of labeled mRNA counts. Dimension: cells x time points. - bounds: tuple - The bound for alpha and beta. The default is alpha / beta > 0. - alpha_0: float - Initial guess for alpha. - beta_0: float - Initial guess for beta. - gamma_0: float - Initial guess for gamma. - - Returns - ------- - alpha: float - The estimated value for alpha. - beta: float - The estimated value for beta. - gamma: float - The estimated value for gamma. + Args: + t: a vector of time points. + l: a matrix of labeled mRNA counts. Dimension: cells x time points. + bounds: the bound for alpha and beta. The default is alpha / beta > 0. + alpha_0: initial guess for alpha. + beta_0: initial guess for beta. + gamma_0: initial guess for gamma. + + Returns: + The estimated value for alpha, beta and gamma. """ l = l.A if issparse(l) else l @@ -828,19 +665,12 @@ def concat_time_series_matrices( ) -> Union[Tuple[np.ndarray, np.ndarray], np.ndarray]: """Concatenate a list of gene x cell matrices into a single matrix. - Arguments - --------- - mats: :class:`~numpy.ndarray` - A list of gene x cell matrices. The length of the list equals the number of time points - t: :class:`~numpy.ndarray` or list - A vector or list of time points - - Returns - ------- - ret_mat: :class:`~numpy.ndarray` - Concatenated gene x cell matrix. - ret_t: :class:`~numpy.ndarray` - A vector of time point for each cell. + Args: + mats: a list of gene x cell matrices. The length of the list equals the number of time points + t: a vector or list of time points. + + Returns: + Concatenated gene x cell matrix (ret_mat) and a vector of time point for each cell (ret_t). """ ret_mat = np.concatenate(mats, axis=1) if t is not None: From 01547c3f121a0dae5b27e162d89246bfa006deb4 Mon Sep 17 00:00:00 2001 From: sichao Date: Mon, 23 Oct 2023 17:20:30 -0400 Subject: [PATCH 32/36] add typing and docstr for fit jacobian --- dynamo/estimation/fit_jacobian.py | 175 ++++++++++++++++++++++++++++-- 1 file changed, 167 insertions(+), 8 deletions(-) diff --git a/dynamo/estimation/fit_jacobian.py b/dynamo/estimation/fit_jacobian.py index f03b936c9..700685036 100644 --- a/dynamo/estimation/fit_jacobian.py +++ b/dynamo/estimation/fit_jacobian.py @@ -1,30 +1,122 @@ +from typing import Callable, Dict, Optional, Tuple, Union + import numpy as np from scipy.optimize import curve_fit from ..dynamo_logger import main_warning -def hill_inh_func(x, A, K, n, g): +def hill_inh_func( + x: np.ndarray, + A: np.ndarray, + K: np.ndarray, + n: np.ndarray, + g: Union[float, int, np.ndarray], +) -> np.ndarray: + """The inhibitory Hill equation function. + + Args: + x: the input transcription factor. + A: the maximum transcription rate. + K: is the ligand concentration producing half occupation. + n: the Hill coefficient. + g: the degradation constant as the offset. + + Returns: + The output of the Hill equation representing the inhibitory effect of the transcription factor. + """ Kd = K**n return A * Kd / (Kd + x**n) - g * x -def hill_inh_grad(x, A, K, n, g): +def hill_inh_grad( + x: np.ndarray, + A: np.ndarray, + K: np.ndarray, + n: np.ndarray, + g: Union[float, int, np.ndarray], +) -> np.ndarray: + """The gradient function of the inhibitory Hill equation. + + Args: + x: the input transcription factor. + A: the maximum transcription rate. + K: is the ligand concentration producing half occupation. + n: the Hill coefficient. + g: the degradation constant as the offset. + + Returns: + The gradient of the Hill equation. + """ Kd = K**n return -A * n * Kd * x ** (n - 1) / (Kd + x**n) ** 2 - g -def hill_act_func(x, A, K, n, g): +def hill_act_func( + x: np.ndarray, + A: np.ndarray, + K: np.ndarray, + n: np.ndarray, + g: Union[float, int, np.ndarray], +) -> np.ndarray: + """The activation Hill equation function. + + Args: + x: the input transcription factor. + A: the maximum transcription rate. + K: is the ligand concentration producing half occupation. + n: the Hill coefficient. + g: the degradation constant as the offset. + + Returns: + The output of the Hill equation representing the activation effect of the transcription factor. + """ Kd = K**n return A * x**n / (Kd + x**n) - g * x -def hill_act_grad(x, A, K, n, g): +def hill_act_grad( + x: np.ndarray, + A: np.ndarray, + K: np.ndarray, + n: np.ndarray, + g: Union[float, int, np.ndarray], +) -> np.ndarray: + """The gradient function of the activation Hill equation. + + Args: + x: the input transcription factor. + A: the maximum transcription rate. + K: is the ligand concentration producing half occupation. + n: the Hill coefficient. + g: the degradation constant as the offset. + + Returns: + The gradient of the Hill equation. + """ Kd = K**n return A * n * Kd * x ** (n - 1) / (Kd + x**n) ** 2 - g -def calc_mean_squared_deviation(func, x_data, y_mean, y_sigm, weighted=True): +def calc_mean_squared_deviation( + func: Callable, + x_data: np.ndarray, + y_mean: np.ndarray, + y_sigm: np.ndarray, + weighted=True, +) -> float: + """Calculate the mean squared deviation of the fit. + + Args: + func: the function to evaluate. + x_data: an array of the x data. + y_mean: an array of the mean of y data. + y_sigm: an array of the standard deviation of the y data. + weighted: whether to use weighted mean squared deviation. + + Returns: + The mean squared deviation. + """ err = func(x_data) - y_mean if weighted: sig = np.array(y_sigm, copy=True) @@ -35,7 +127,34 @@ def calc_mean_squared_deviation(func, x_data, y_mean, y_sigm, weighted=True): return np.sqrt(err.dot(err)) -def fit_hill_grad(x_data, y_mean, type, y_sigm=None, fix_g=None, n_num=5, x_tol=1e-5, x_shift=0): +def fit_hill_grad( + x_data: np.ndarray, + y_mean: np.ndarray, + type: str, + y_sigm: Optional[np.ndarray] = None, + fix_g: Optional[Union[float, int, np.ndarray]] = None, + n_num: int = 5, + x_tol: float = 1e-5, + x_shift: int = 0, +) -> Tuple: + """Fit the Hill equation to the data. + + This can reveal whether a gene is regulated (activated or suppressed) by certain transcription factor through + fitting the derivatives of inhibitory or activation Hill equations to the corresponding Jacobian elements. + + Args: + x_data: an array of the x data. + y_mean: an array of the mean of y data. + type: the type of the Hill equation, either `act` or `inh`. + y_sigm: an array of the standard deviation of the y data. + fix_g: the degradation constant as the offset. + n_num: the number of Hill coefficients to try. + x_tol: the tolerance of the x data. + x_shift: the shift of the x data. + + Returns: + The fitted parameters and the mean squared deviation. + """ assert type in ["act", "inh"], "`type` must be either `act` or `inh`." A0 = -np.min(y_mean) if type == "inh" else np.max(y_mean) @@ -92,7 +211,27 @@ def fit_hill_grad(x_data, y_mean, type, y_sigm=None, fix_g=None, n_num=5, x_tol= return {"A": A, "K": K, "n": n, "g": g}, msd_min -def fit_hill_inh_grad(x_data, y_mean, y_sigm=None, n_num=5, x_tol=1e-5, x_shift=1e-4): +def fit_hill_inh_grad( + x_data: np.ndarray, + y_mean: np.ndarray, + y_sigm: Optional[np.ndarray] = None, + n_num: int = 5, + x_tol: float = 1e-5, + x_shift: float = 1e-4, +) -> Tuple[Dict, float]: + """Fit the inhibitory Hill equation to the data. + + Args: + x_data: an array of the x data. + y_mean: an array of the mean of y data. + y_sigm: an array of the standard deviation of the y data. + n_num: the number of Hill coefficients to try. + x_tol: the tolerance of the x data. + x_shift: the shift of the x data. + + Returns: + The fitted parameters and the mean squared deviation. + """ A0 = -np.min(y_mean) K0 = (x_data[0] + x_data[-1]) / 2 g0 = 0 @@ -132,7 +271,27 @@ def fit_hill_inh_grad(x_data, y_mean, y_sigm=None, n_num=5, x_tol=1e-5, x_shift= return {"A": p_opt_min[0], "K": p_opt_min[1], "n": np.exp(p_opt_min[2]), "g": p_opt_min[3]}, msd_min -def fit_hill_act_grad(x_data, y_mean, y_sigm=None, n_num=5, x_tol=1e-5, x_shift=1e-4): +def fit_hill_act_grad( + x_data: np.ndarray, + y_mean: np.ndarray, + y_sigm: Optional[np.ndarray] = None, + n_num: int = 5, + x_tol: float = 1e-5, + x_shift: float = 1e-4, +) -> Tuple[Dict, float]: + """Fit the activation Hill equation to the data. + + Args: + x_data: an array of the x data. + y_mean: an array of the mean of y data. + y_sigm: an array of the standard deviation of the y data. + n_num: the number of Hill coefficients to try. + x_tol: the tolerance of the x data. + x_shift: the shift of the x data. + + Returns: + The fitted parameters and the mean squared deviation. + """ A0 = np.max(y_mean) K0 = (x_data[0] + x_data[-1]) / 2 g0 = 0 From 0043b469cd10ec8fa018486c1dd8e90137e26fd3 Mon Sep 17 00:00:00 2001 From: sichao Date: Mon, 23 Oct 2023 17:58:43 -0400 Subject: [PATCH 33/36] fix docstr error and format --- dynamo/estimation/csc/velocity.py | 121 +++++++++----------- dynamo/estimation/tsc/estimation_kinetic.py | 58 ++++++++-- dynamo/estimation/tsc/utils_kinetic.py | 55 ++++++--- 3 files changed, 143 insertions(+), 91 deletions(-) diff --git a/dynamo/estimation/csc/velocity.py b/dynamo/estimation/csc/velocity.py index 8588836fa..59cb7c12b 100755 --- a/dynamo/estimation/csc/velocity.py +++ b/dynamo/estimation/csc/velocity.py @@ -22,36 +22,32 @@ class Velocity: - """The class that computes RNA/protein velocity given unknown parameters. - - Arguments - --------- - alpha: :class:`~numpy.ndarray` - A matrix of transcription rate. - beta: :class:`~numpy.ndarray` - A vector of splicing rate constant for each gene. - gamma: :class:`~numpy.ndarray` - A vector of spliced mRNA degradation rate constant for each gene. - eta: :class:`~numpy.ndarray` - A vector of protein synthesis rate constant for each gene. - delta: :class:`~numpy.ndarray` - A vector of protein degradation rate constant for each gene. - t: :class:`~numpy.ndarray` or None (default: None) - A vector of the measured time points for cells - estimation: :class:`~ss_estimation` - An instance of the estimation class. If this not None, the parameters will be taken from this class instead of the input arguments. - """ - + """The class that computes RNA/protein velocity given unknown parameters.""" def __init__( self, - alpha=None, - beta=None, - gamma=None, - eta=None, - delta=None, - t=None, + alpha: Optional[np.ndarray] = None, + beta: Optional[np.ndarray] = None, + gamma: Optional[np.ndarray] = None, + eta: Optional[np.ndarray] = None, + delta: Optional[np.ndarray] = None, + t: Optional[np.ndarray] = None, estimation=None, ): + """Initialize the velocity class. + + Args: + alpha: a matrix of transcription rate. + beta: a vector of splicing rate constant for each gene. + gamma: a vector of spliced mRNA degradation rate constant for each gene. + eta: a vector of protein synthesis rate constant for each gene. + delta: a vector of protein degradation rate constant for each gene. + t: a vector of the measured time points for cells + estimation: a ss_estimation instance of the estimation class. If this not None, the parameters will be + taken from this class instead of the input arguments. + + Returns: + A velocity instance. + """ if estimation is not None: self.parameters = {} self.parameters["alpha"] = estimation.parameters["alpha"] @@ -70,21 +66,21 @@ def __init__( "t": t, } - def vel_u(self, U, repeat=None, update_alpha=True): + def vel_u( + self, + U: Union[csr_matrix, np.ndarray], + repeat: Optional[bool] = None, + update_alpha: bool = True, + ) -> Union[csr_matrix, np.ndarray]: """Calculate the unspliced mRNA velocity. - Arguments - --------- - U: :class:`~numpy.ndarray` or sparse `csr_matrix` - A matrix of unspliced mRNA count. Dimension: genes x cells. - repeat: bool or None - Whether to use average alpha or cell-wise alpha with the formula: + Args: + U: a matrix of unspliced mRNA count. Dimension: genes x cells. + repeat: whether to use average alpha or cell-wise alpha with the formula: $a = \frac{n \gamma}{1 - e^{-\gamma t}}$. - Returns - ------- - V: :class:`~numpy.ndarray` or sparse `csr_matrix` - Each column of V is a velocity vector for the corresponding cell. Dimension: genes x cells. + Returns: + Each column of V is a velocity vector for the corresponding cell. Dimension: genes x cells. """ t = self.parameters["t"] @@ -167,20 +163,15 @@ def vel_u(self, U, repeat=None, update_alpha=True): return V - def vel_s(self, U, S): + def vel_s(self, U: Union[csr_matrix, np.ndarray], S: Union[csr_matrix, np.ndarray]) -> Union[csr_matrix, np.ndarray]: """Calculate the unspliced mRNA velocity. - Arguments - --------- - U: :class:`~numpy.ndarray` or sparse `csr_matrix` - A matrix of unspliced mRNA counts. Dimension: genes x cells. - S: :class:`~numpy.ndarray` or sparse `csr_matrix` - A matrix of spliced mRNA counts. Dimension: genes x cells. + Args: + U: a matrix of unspliced mRNA counts. Dimension: genes x cells. + S: a matrix of spliced mRNA counts. Dimension: genes x cells. - Returns - ------- - V: :class:`~numpy.ndarray` or sparse `csr_matrix` - Each column of V is a velocity vector for the corresponding cell. Dimension: genes x cells. + Returns: + Each column of V is a velocity vector for the corresponding cell. Dimension: genes x cells. """ t = self.parameters["t"] @@ -236,20 +227,15 @@ def vel_s(self, U, S): V = np.nan return V - def vel_p(self, S, P): + def vel_p(self, S: Union[csr_matrix, np.ndarray], P: Union[csr_matrix, np.ndarray]) -> Union[csr_matrix, np.ndarray]: """Calculate the protein velocity. - Arguments - --------- - S: :class:`~numpy.ndarray` or sparse `csr_matrix` - A matrix of spliced mRNA counts. Dimension: genes x cells. - P: :class:`~numpy.ndarray` or sparse `csr_matrix` - A matrix of protein counts. Dimension: genes x cells. + Args: + S: a matrix of spliced mRNA counts. Dimension: genes x cells. + P: a matrix of protein counts. Dimension: genes x cells. - Returns - ------- - V: :class:`~numpy.ndarray` or sparse `csr_matrix` - Each column of V is a velocity vector for the corresponding cell. Dimension: genes x cells. + Returns: + Each column of V is a velocity vector for the corresponding cell. Dimension: genes x cells. """ t = self.parameters["t"] @@ -292,13 +278,11 @@ def vel_p(self, S, P): V = np.nan return V - def get_n_cells(self): + def get_n_cells(self) -> int: """Get the number of cells if the parameter alpha is given. - Returns - ------- - n_cells: int - The second dimension of the alpha matrix, if alpha is given. + Returns: + The second dimension of the alpha matrix, if alpha is given. """ if self.parameters["alpha"] is not None: n_cells = self.parameters["alpha"].shape[1] @@ -306,13 +290,12 @@ def get_n_cells(self): n_cells = np.nan return n_cells - def get_n_genes(self): + def get_n_genes(self) -> int: """Get the number of genes. - Returns - ------- - n_genes: int - The first dimension of the alpha matrix, if alpha is given. Or, the length of beta, gamma, eta, or delta, if they are given. + Returns: + The first dimension of the alpha matrix, if alpha is given. Or, the length of beta, gamma, eta, or delta, + if they are given. """ if self.parameters["alpha"] is not None: n_genes = self.parameters["alpha"].shape[0] diff --git a/dynamo/estimation/tsc/estimation_kinetic.py b/dynamo/estimation/tsc/estimation_kinetic.py index c420a5095..25d9caf24 100755 --- a/dynamo/estimation/tsc/estimation_kinetic.py +++ b/dynamo/estimation/tsc/estimation_kinetic.py @@ -77,7 +77,7 @@ def guestimate_p0_kinetic_chase(x_data: np.ndarray, time: np.ndarray) -> Tuple[n class kinetic_estimation: """A general parameter estimation framework for all types of time-seris data.""" - def __init__(self, param_ranges: np.ndarray, x0_ranges: np.ndarray, simulator: LinearODE) -> None: + def __init__(self, param_ranges: np.ndarray, x0_ranges: np.ndarray, simulator: LinearODE): """Initialize the kinetic_estimation class. Args: @@ -88,6 +88,9 @@ def __init__(self, param_ranges: np.ndarray, x0_ranges: np.ndarray, simulator: L simulator: An instance of python class which solves ODEs. It should have properties 't' (k time points, 1d numpy array), 'x0' (initial conditions for m species, 1d numpy array), and 'x' (solution, k-by-m array), as well as two functions: integrate (numerical integration), solve (analytical method). + + Returns: + An instance of the kinetic_estimation class. """ self.simulator = simulator @@ -409,13 +412,16 @@ def test_chi2( class Estimation_Degradation(kinetic_estimation): """The base parameters, estimation class for degradation experiments.""" - def __init__(self, ranges: np.ndarray, x0: np.ndarray, simulator: LinearODE) -> None: + def __init__(self, ranges: np.ndarray, x0: np.ndarray, simulator: LinearODE): """Initialize the Estimation_Degradation object. Args: ranges: the lower and upper ranges of parameters. x0: initial conditions. simulator: instance of the Python class to solve ODEs. + + Returns: + An instance of the Estimation_Degradation class. """ self.kin_param_keys = np.array(["alpha", "gamma"]) super().__init__(np.vstack((np.zeros(2), ranges)), x0, simulator) @@ -493,13 +499,16 @@ def __init__( beta: Optional[np.ndarray] = None, gamma: Optional[np.ndarray] = None, x0: Optional[np.ndarray] = None, - ) -> None: + ): """Initialize the Estimation_DeterministicDeg object. Args: beta: the splicing rate. gamma: the degradation rate. x0: the initial conditions. + + Returns: + An instance of the Estimation_DeterministicDeg class. """ self.kin_param_keys = np.array(["alpha", "beta", "gamma"]) if beta is not None and gamma is not None and x0 is not None: @@ -544,12 +553,15 @@ def auto_fit(self, time: np.ndarray, x_data: np.ndarray, **kwargs) -> Tuple[np.n class Estimation_DeterministicDegNosp(Estimation_Degradation): """An estimation class for degradation (without splicing) experiments.""" - def __init__(self, gamma: Optional[np.ndarray] = None, x0: Optional[np.ndarray] = None) -> None: + def __init__(self, gamma: Optional[np.ndarray] = None, x0: Optional[np.ndarray] = None): """Initialize the Estimation_DeterministicDegNosp object. Args: gamma: the degradation rate. x0: the initial conditions. + + Returns: + An instance of the Estimation_DeterministicDegNosp class. """ if gamma is not None and x0 is not None: self._initialize(gamma, x0) @@ -619,7 +631,7 @@ def __init__( gamma: Optional[np.ndarray] = None, x0: Optional[np.ndarray] = None, include_cov: bool = True, - ) -> None: + ): """Initialize the Estimation_MomentDeg object. Args: @@ -627,6 +639,9 @@ def __init__( gamma: the degradation rate. x0: the initial conditions. include_cov: whether to consider covariance when estimating. + + Returns: + An instance of the Estimation_MomentDeg class. """ self.kin_param_keys = np.array(["alpha", "beta", "gamma"]) self.include_cov = include_cov @@ -667,12 +682,15 @@ def extract_data_from_simulator(self) -> np.ndarray: class Estimation_MomentDegNosp(Estimation_Degradation): """An estimation class for degradation (without splicing) experiments. Order of species: , .""" - def __init__(self, gamma: Optional[np.ndarray] = None, x0: Optional[np.ndarray] = None) -> None: + def __init__(self, gamma: Optional[np.ndarray] = None, x0: Optional[np.ndarray] = None): """Initialize the Estimation_MomentDeg object. Args: gamma: the degradation rate. x0: the initial conditions. + + Returns: + An instance of the Estimation_MomentDeg class. """ if gamma is not None and x0 is not None: self._initialize(gamma, x0) @@ -742,6 +760,9 @@ def __init__(self, a, b, alpha_a, alpha_i, beta, gamma, include_cov=True): beta: splicing rate. gamma: degradation rate. include_cov: whether to include the covariance when estimating. + + Returns: + An instance of the Estimation_MomentKin class. """ self.param_keys = np.array(["a", "b", "alpha_a", "alpha_i", "beta", "gamma"]) ranges = np.zeros((6, 2)) @@ -832,6 +853,9 @@ def __init__(self, a, b, alpha_a, alpha_i, gamma): alpha_a: transcription rate for active promoter. alpha_i: transcription rate for inactive promoter. gamma: degradation rate. + + Returns: + An instance of the Estimation_MomentKinNosp class. """ self.param_keys = np.array(["a", "b", "alpha_a", "alpha_i", "gamma"]) ranges = np.zeros((5, 2)) @@ -900,6 +924,9 @@ def __init__(self, alpha, gamma, x0=0): alpha: transcription rate. gamma: degradation rate. x0: the initial condition. + + Returns: + An instance of the Estimation_DeterministicKinNosp class. """ self.param_keys = np.array(["alpha", "gamma"]) ranges = np.zeros((2, 2)) @@ -952,6 +979,9 @@ def __init__(self, alpha, beta, gamma, x0=np.zeros(2)): beta: splicing rate. gamma: degradation rate. x0: the initial condition. + + Returns: + An instance of the Estimation_DeterministicKin class. """ self.param_keys = np.array(["alpha", "beta", "gamma"]) ranges = np.zeros((3, 2)) @@ -1017,6 +1047,9 @@ def __init__(self, model1, model2, alpha=None, gamma=None, x0=None, beta=None): gamma: degradation rate. x0: the initial condition. beta: splicing rate. + + Returns: + An instance of the Mixture_KinDeg_NoSwitching class. """ self.model1 = model1 self.model2 = model2 @@ -1171,8 +1204,8 @@ def export_dictionary(self): class Lambda_NoSwitching(Mixture_KinDeg_NoSwitching): - """An estimation class with the mixture model. - If beta is None, it is assumed that the data does not have the splicing process. + """An estimation class with the mixture model. If beta is None, it is assumed that the data does not have the + splicing process. """ def __init__( @@ -1195,6 +1228,9 @@ def __init__( gamma: degradation rate. x0: the initial condition. beta: splicing rate. + + Returns: + An instance of the Lambda_NoSwitching class. """ self.model1 = model1 self.model2 = model2 @@ -1265,6 +1301,9 @@ def __init__(self, alpha=None, gamma=None, x0=None): alpha: transcription rate. gamma: degradation rate. x0: the initial condition. + + Returns: + An instance of the Estimation_KineticChase class. """ self.kin_param_keys = np.array(["alpha", "gamma"]) if alpha is not None and gamma is not None and x0 is not None: @@ -1354,6 +1393,9 @@ def __init__(self, simulator, params=None, x0=None): simulator: the linearODE class. params: the parameters. x0: the initial conditions. + + Returns: + An instance of the GoodnessOfFit class. """ self.simulator = simulator if params is not None: diff --git a/dynamo/estimation/tsc/utils_kinetic.py b/dynamo/estimation/tsc/utils_kinetic.py index 046afbcff..05277ba3f 100755 --- a/dynamo/estimation/tsc/utils_kinetic.py +++ b/dynamo/estimation/tsc/utils_kinetic.py @@ -9,12 +9,15 @@ class LinearODE: """A general class for linear odes.""" - def __init__(self, n_species: int, x0: Optional[np.ndarray] = None) -> None: + def __init__(self, n_species: int, x0: Optional[np.ndarray] = None): """Initialize the LinearODE object. Args: n_species: the number of species. x0: the initial condition of variable x. + + Returns: + An instance of LinearODE. """ self.n_species = n_species # solution @@ -32,7 +35,7 @@ def ode_func(self, x: np.ndarray, t: np.ndarray) -> np.ndarray: Args: x: the variable. - t: the aray of time. + t: the array of time. Returns: The derivatives dx. @@ -134,12 +137,15 @@ def integrate_matrix(self, t: np.ndarray, x0: Optional[np.ndarray] = None) -> np class MixtureModels: """The base class for mixture models.""" - def __init__(self, models: LinearODE, param_distributor: List) -> None: + def __init__(self, models: LinearODE, param_distributor: List): """Initialize the MixtureModels class. Args: models: the models to mix. param_distributor: the index to assign parameters. + + Returns: + An instance of MixtureModels. """ self.n_models = len(models) self.models = models @@ -220,12 +226,15 @@ class LambdaModels_NoSwitching(MixtureModels): """Linear ODEs for the lambda mixture model. The order of params is: parameter order: alpha, lambda, (beta), gamma distributor order: alpha_1, alpha_2, (beta), gamma""" - def __init__(self, model1: LinearODE, model2: LinearODE) -> None: + def __init__(self, model1: LinearODE, model2: LinearODE): """Initialize the LambdaModels_NoSwitching class. Args: model1: the first model to mix. model2: the second model to mix. + + Returns: + An instance of LambdaModels_NoSwitching. """ models = [model1, model2] if type(model1) in nosplicing_models and type(model2) in nosplicing_models: @@ -261,7 +270,7 @@ def __init__( beta: Optional[np.ndarray] = None, gamma: Optional[np.ndarray] = None, x0: Optional[np.ndarray] = None, - ) -> None: + ): """Initialize the Moments object. Args: @@ -272,6 +281,9 @@ def __init__( beta: splicing rate. gamma: degradation rate. x0: the initial conditions. + + Returns: + An instance of Moments. """ # species self.ua = 0 @@ -301,7 +313,7 @@ def ode_func(self, x: np.ndarray, t: np.ndarray) -> np.ndarray: Args: x: the variable. - t: the aray of time. + t: the array of time. Returns: The derivatives dx. @@ -523,7 +535,7 @@ def __init__( alpha_i: Optional[np.ndarray] = None, gamma: Optional[np.ndarray] = None, x0: Optional[np.ndarray] = None - ) -> None: + ): """Initialize the Moments_Nosplicing object. Args: @@ -533,6 +545,9 @@ def __init__( alpha_i: transcription rate for inactive promoter. gamma: degradation rate. x0: the initial conditions. + + Returns: + An instance of Moments_Nosplicing. """ # species self.ua = 0 @@ -553,7 +568,7 @@ def ode_func(self, x: np.ndarray, t: np.ndarray) -> np.ndarray: Args: x: the variable. - t: the aray of time. + t: the array of time. Returns: The derivatives dx. @@ -677,7 +692,7 @@ def __init__( beta: Optional[np.ndarray] = None, gamma: Optional[np.ndarray] = None, x0: Optional[np.ndarray] = None, - ) -> None: + ): """Initialize the Moments_NoSwitching object. Args: @@ -685,6 +700,9 @@ def __init__( beta: splicing rate. gamma: degradation rate. x0: the initial conditions. + + Returns: + An instance of Moments_NoSwitching. """ # species self.u = 0 @@ -707,7 +725,7 @@ def ode_func(self, x: np.ndarray, t: np.ndarray) -> np.ndarray: Args: x: the variable. - t: the aray of time. + t: the array of time. Returns: The derivatives dx. @@ -862,13 +880,16 @@ def __init__( alpha: Optional[np.ndarray] = None, gamma: Optional[np.ndarray] = None, x0: Optional[np.ndarray] = None, - ) -> None: + ): """Initialize the Moments_NoSwitchingNoSplicing object. Args: alpha: transcription rate. gamma: degradation rate. x0: the initial conditions. + + Returns: + An instance of Moments_NoSwitchingNoSplicing. """ # species self.u = 0 @@ -888,7 +909,7 @@ def ode_func(self, x: np.ndarray, t: np.ndarray) -> np.ndarray: Args: x: the variable. - t: the aray of time. + t: the array of time. Returns: The derivatives dx. @@ -982,7 +1003,7 @@ def __init__( beta: Optional[np.ndarray] = None, gamma: Optional[np.ndarray] = None, x0: Optional[np.ndarray] = None, - ) -> None: + ): """Initialize the Deterministic object. Args: @@ -990,6 +1011,9 @@ def __init__( beta: the splicing rate. gamma: the degradation rate. x0: the inital conditions. + + Returns: + An instance of Deterministic. """ # species self.u = 0 @@ -1112,13 +1136,16 @@ def __init__( alpha: Optional[np.ndarray] = None, gamma: Optional[np.ndarray] = None, x0: Optional[np.ndarray] = None, - ) -> None: + ): """Initialize the Deterministic_NoSplicing object. Args: alpha: the transcription rate. gamma: the degradation rate. x0: the initial conditions. + + Returns: + An instance of Deterministic_NoSplicing. """ # species self.u = 0 From 43ab31bf87befdc843cf33ba096af0ad59c27826 Mon Sep 17 00:00:00 2001 From: sichao Date: Thu, 2 Nov 2023 18:55:30 -0400 Subject: [PATCH 34/36] Capitalize arguments --- dynamo/estimation/csc/utils_velocity.py | 221 +++++++------- dynamo/estimation/csc/velocity.py | 26 +- dynamo/estimation/fit_jacobian.py | 90 +++--- dynamo/estimation/tsc/estimation_kinetic.py | 312 ++++++++++---------- dynamo/estimation/tsc/twostep.py | 44 +-- dynamo/estimation/tsc/utils_kinetic.py | 192 ++++++------ 6 files changed, 442 insertions(+), 443 deletions(-) diff --git a/dynamo/estimation/csc/utils_velocity.py b/dynamo/estimation/csc/utils_velocity.py index dceb0b3d8..3040f61e0 100755 --- a/dynamo/estimation/csc/utils_velocity.py +++ b/dynamo/estimation/csc/utils_velocity.py @@ -14,10 +14,10 @@ def sol_u(t: np.ndarray, u0: float, alpha: float, beta: float) -> np.ndarray: """The analytical solution of unspliced mRNA kinetics. Args: - t: a vector of time points. - u0: initial value of u. - alpha: transcription rate. - beta: splicing rate constant. + t: A vector of time points. + u0: Initial value of u. + alpha: Transcription rate. + beta: Splicing rate constant. Returns: Unspliced mRNA counts at given time points. @@ -29,12 +29,12 @@ def sol_u_2p(t: np.ndarray, u0: float, t1: np.ndarray, alpha0: float, alpha1: fl """The combined 2-piece analytical solution of unspliced mRNA kinetics. Args: - t: a vector of time points for both steady state and stimulation labeling. - u0: initial value of u. - t1: the time point when the cells switch from steady state to stimulation. - alpha0: transcription rate for steady state labeling. - alpha1: transcription rate for stimulation based labeling. - beta: splicing rate constant. + t: A vector of time points for both steady state and stimulation labeling. + u0: Initial value of u. + t1: The time point when the cells switch from steady state to stimulation. + alpha0: Transcription rate for steady state labeling. + alpha1: Transcription rate for stimulation based labeling. + beta: Splicing rate constant. Returns: Unspliced mRNA counts at given time points. @@ -50,12 +50,12 @@ def sol_s(t: np.ndarray, s0: float, u0: float, alpha: float, beta: float, gamma: """The analytical solution of spliced mRNA kinetics. Args: - t: a vector of time points. - s0: initial value of s. - u0: initial value of u. - alpha: transcription rate. - beta: splicing rate constant. - gamma: degradation rate constant for spliced mRNA. + t: A vector of time points. + s0: Initial value of s. + u0: Initial value of u. + alpha: Transcription rate. + beta: Splicing rate constant. + gamma: Degradation rate constant for spliced mRNA. Returns: Spliced mRNA counts at given time points. @@ -86,15 +86,15 @@ def sol_p( """The analytical solution of protein kinetics. Args: - t: a vector of time points. - p0: initial value of p. - s0: initial value of s. - u0: initial value of u. - alpha: transcription rate. - beta: splicing rate constant. - gamma: degradation rate constant for spliced mRNA. - eta: synthesis rate constant for protein. - delta: degradation rate constant for protein. + t: A vector of time points. + p0: Initial value of p. + s0: Initial value of s. + u0: Initial value of u. + alpha: Transcription rate. + beta: Splicing rate constant. + gamma: Degradation rate constant for spliced mRNA. + eta: Synthesis rate constant for protein. + delta: Degradation rate constant for protein. Returns: Protein counts at given time points (p), spliced mRNA counts at given time points (s) and unspliced mRNA counts @@ -113,9 +113,9 @@ def solve_gamma(t: float, old: np.ndarray, total: np.ndarray) -> np.ndarray: """Analytical solution to calculate gamma (degradation rate) using first-order degradation kinetics. Args: - t: metabolic labeling time period. - old: a vector of old RNA amount in each cell - total: a vector of total RNA amount in each cell + t: Metabolic labeling time period. + old: A vector of old RNA amount in each cell + total: A vector of total RNA amount in each cell Returns: The degradation rate (gamma). @@ -132,11 +132,11 @@ def solve_alpha_2p(t0: float, t1: float, alpha0: float, beta: float, u1: Union[c labeling experiment. Args: - t0: time period for steady state labeling. - t1: time period for stimulation labeling. - alpha0: steady state transcription rate calculated from one-shot experiment mode. - beta: steady state (and simulation) splicing rate calculated from one-shot experiment mode. - u1: a vector of labeled RNA amount in each cell observed at time t0 + t1. + t0: Time period for steady state labeling. + t1: Time period for stimulation labeling. + alpha0: Steady state transcription rate calculated from one-shot experiment mode. + beta: Steady state (and simulation) splicing rate calculated from one-shot experiment mode. + u1: A vector of labeled RNA amount in each cell observed at time t0 + t1. Returns: The transcription rate (alpha1) for the stimulation period in the data. @@ -165,11 +165,11 @@ def solve_alpha_2p_mat( labeling experiment in a matrix form. Args: - t0: time period for steady state labeling for each cell. - t1: time period for stimulation labeling for each cell. - alpha0: steady state transcription rate calculated from one-shot experiment mode for each gene. - beta: steady state (and simulation) splicing rate calculated from one-shot experiment mode for each gene. - u1: a vector of labeled RNA amount in each cell observed at time t0 + t1. + t0: Time period for steady state labeling for each cell. + t1: Time period for stimulation labeling for each cell. + alpha0: Steady state transcription rate calculated from one-shot experiment mode for each gene. + beta: Steady state (and simulation) splicing rate calculated from one-shot experiment mode for each gene. + u1: A vector of labeled RNA amount in each cell observed at time t0 + t1. Returns: The transcription rate (alpha1) for the stimulation period in the data. @@ -202,13 +202,13 @@ def fit_linreg( """Simple linear regression: y = kx + b. Args: - x: a vector of independent variables. - y: a vector of dependent variables. - mask: the index to filter the data before regression. - intercept: if using steady state assumption for fitting, then: + x: A vector of independent variables. + y: A vector of dependent variables. + mask: The index to filter the data before regression. + intercept: If using steady state assumption for fitting, then: True -- the linear regression is performed with an unfixed intercept; False -- the linear regression is performed with a fixed zero intercept. - r2: the R squared value calculated from the regression. + r2: The R squared value calculated from the regression. Returns: The estimated slope (k) , the estimated intercept (b), coefficient of determination or r square calculated with @@ -266,12 +266,12 @@ def fit_linreg_robust( """Apply robust linear regression of y w.r.t x. Args: - x: a vector of independent variables. - y: a vector of dependent variables. - intercept: if using steady state assumption for fitting, then: + x: A vector of independent variables. + y: A vector of dependent variables. + intercept: If using steady state assumption for fitting, then: True -- the linear regression is performed with an unfixed intercept; False -- the linear regression is performed with a fixed zero intercept. - est_method: the linear regression estimation method that will be used. + est_method: The linear regression estimation method that will be used. Returns: The estimated slope (k), the estimated intercept (b), coefficient of determination or r square calculated with @@ -345,10 +345,10 @@ def fit_stochastic_linreg( """Generalized method of moments: [u, 2*us + u] = gamma * [s, 2*ss - s]. Args: - u: a vector of first moments (mean) of unspliced (or new) RNA expression. - s: a vector of first moments (mean) of spliced (or total) RNA expression. - us: a vector of second moments (uncentered co-variance) of unspliced/spliced (or new/total) RNA expression. - ss: a vector of second moments (uncentered variance) of spliced (or total) RNA expression. + u: A vector of first moments (mean) of unspliced (or new) RNA expression. + s: A vector of first moments (mean) of spliced (or total) RNA expression. + us: A vector of second moments (uncentered co-variance) of unspliced/spliced (or new/total) RNA expression. + ss: A vector of second moments (uncentered variance) of spliced (or total) RNA expression. Returns: The estimated gamma. @@ -395,12 +395,12 @@ def fit_first_order_deg_lsq( """Estimate beta with degradation data using least squares method. Args: - t: a vector of time points. - l: a vector of unspliced, labeled mRNA counts for each time point. - bounds: the bound for beta. The default is beta > 0. - fix_l0: whether to used fixed l0. If ture, l0 will be calculated by averaging the first column of l; + t: A vector of time points. + l: A vector of unspliced, labeled mRNA counts for each time point. + bounds: The bound for beta. The default is beta > 0. + fix_l0: Whether to used fixed l0. If ture, l0 will be calculated by averaging the first column of l; If false, l0 is a parameter that will be estimated all together with beta using lsq. - beta_0: initial guess for beta. + beta_0: Initial guess for beta. Returns: The estimated value for beta (beta) and the estimated value for the initial spliced, labeled mRNA count (l0). @@ -427,8 +427,8 @@ def solve_first_order_deg(t: np.ndarray, l: Union[csr_matrix, np.ndarray]) -> Tu under first-order degradation kinetics model. Args: - t: a vector of time points. - l: a vector of labeled mRNA counts for each time point. + t: A vector of time points. + l: A vector of labeled mRNA counts for each time point. Returns: The initial counts of the species (for example, labeled mRNA), degradation rate constant and half-life the species. @@ -457,12 +457,12 @@ def fit_gamma_lsq( """Estimate gamma with degradation data using least squares method. Args: - t: a vector of time points. - s: a vector of spliced, labeled mRNA counts for each time point. - beta: the value of beta. - u0: initial number of unspliced mRNA. - bounds: the bound for gamma. The default is gamma > 0. - fix_s0: whether to use fixed s0. If true, s0 will be calculated by averaging the first column of s; + t: A vector of time points. + s: A vector of spliced, labeled mRNA counts for each time point. + beta: The value of beta. + u0: Initial number of unspliced mRNA. + bounds: The bound for gamma. The default is gamma > 0. + fix_s0: Whether to use fixed s0. If true, s0 will be calculated by averaging the first column of s; If false, s0 is a parameter that will be estimated all together with gamma using lsq. Returns: @@ -494,9 +494,9 @@ def fit_alpha_synthesis(t: np.ndarray, u: Union[csr_matrix, np.ndarray], beta: f It is assumed that u(0) = 0. Args: - t: a vector of time points. - u: a matrix of unspliced mRNA counts. Dimension: cells x time points. - beta: the value of beta. + t: A vector of time points. + u: A matrix of unspliced mRNA counts. Dimension: cells x time points. + beta: The value of beta. Returns: The estimated value for alpha. @@ -522,10 +522,10 @@ def fit_alpha_degradation( constrains u0 to be larger than 0. Args: - t: a vector of time points. - u: a matrix of unspliced mRNA counts. Dimension: cells x time points. - beta: the value of beta. - intercept: whether to steady state assumption for fitting, then: + t: A vector of time points. + u: A matrix of unspliced mRNA counts. Dimension: cells x time points. + beta: The value of beta. + intercept: Whether to steady state assumption for fitting, then: True -- the linear regression is performed with an unfixed intercept; False -- the linear regression is performed with a fixed zero intercept. @@ -558,10 +558,10 @@ def solve_alpha_degradation( """Estimate alpha with degradation data using linear regression. Args: - t: a vector of time points. - u: a matrix of unspliced mRNA counts. Dimension: cells x time points. - beta: the value of beta. - intercept: if using steady state assumption for fitting, then: + t: A vector of time points. + u: A matrix of unspliced mRNA counts. Dimension: cells x time points. + beta: The value of beta. + intercept: If using steady state assumption for fitting, then: True -- the linear regression is performed with an unfixed intercept; False -- the linear regression is performed with a fixed zero intercept. @@ -608,11 +608,11 @@ def fit_alpha_beta_synthesis( It is assumed that u(0) = 0. Args: - t: a vector of time points. - l: a matrix of labeled mRNA counts. Dimension: cells x time points. - bounds: the bound for alpha and beta. The default is alpha / beta > 0. - alpha_0: initial guess for alpha. - beta_0: initial guess for beta. + t: A vector of time points. + l: A matrix of labeled mRNA counts. Dimension: cells x time points. + bounds: The bound for alpha and beta. The default is alpha / beta > 0. + alpha_0: Initial guess for alpha. + beta_0: Initial guess for beta. Returns: The estimated value for alpha and the estimated value for beta. @@ -639,12 +639,12 @@ def fit_all_synthesis( It is assumed that u(0) = 0 and s(0) = 0. Args: - t: a vector of time points. - l: a matrix of labeled mRNA counts. Dimension: cells x time points. - bounds: the bound for alpha and beta. The default is alpha / beta > 0. - alpha_0: initial guess for alpha. - beta_0: initial guess for beta. - gamma_0: initial guess for gamma. + t: A vector of time points. + l: A matrix of labeled mRNA counts. Dimension: cells x time points. + bounds: The bound for alpha and beta. The default is alpha / beta > 0. + alpha_0: Initial guess for alpha. + beta_0: Initial guess for beta. + gamma_0: Initial guess for gamma. Returns: The estimated value for alpha, beta and gamma. @@ -666,8 +666,8 @@ def concat_time_series_matrices( """Concatenate a list of gene x cell matrices into a single matrix. Args: - mats: a list of gene x cell matrices. The length of the list equals the number of time points - t: a vector or list of time points. + mats: A list of gene x cell matrices. The length of the list equals the number of time points + t: A vector or list of time points. Returns: Concatenated gene x cell matrix (ret_mat) and a vector of time point for each cell (ret_t). @@ -687,8 +687,8 @@ def compute_dispersion(mX: Union[np.ndarray, csr_matrix], varX: Union[np.ndarray variance = mean + phi * mean^2 Args: - mX: the mean of variable X. - varX: the variance of variable X. + mX: The mean of variable X. + varX: The variance of variable X. Returns: The dispersion parameter. @@ -710,12 +710,12 @@ def fit_k_negative_binomial( k^2 * variance(r) = k * n - phi * n * 2 Args: - n: unspliced mRNA or new/labeled mRNA of one gene. - r: spliced mRNA or total mRNA of one gene. - var: the variance of r. - phi: the dispersion parameter of the negative binomial distribution of one gene. - k0: the initial guess for k. - return_k0: whether to return k0. + n: Unspliced mRNA or new/labeled mRNA of one gene. + r: Spliced mRNA or total mRNA of one gene. + var: The variance of r. + phi: The dispersion parameter of the negative binomial distribution of one gene. + k0: The initial guess for k. + return_k0: Whether to return k0. Returns: The resolution from the least squares optimizer. @@ -744,12 +744,12 @@ def fit_K_negbin( """Fit parameter K of negative binomial distribution to each gene. Args: - N: unspliced mRNA or new/labeled mRNA. - R: spliced mRNA or total mRNA. - varR: the variance of R. - perc_left: left percentile for selecting extreme data points. - perc_right: right percentile for selecting extreme data points. - return_phi: whether to return dispersion parameter Phi. + N: Unspliced mRNA or new/labeled mRNA. + R: Spliced mRNA or total mRNA. + varR: The variance of R. + perc_left: Left percentile for selecting extreme data points. + perc_right: Right percentile for selecting extreme data points. + return_phi: Whether to return dispersion parameter Phi. Returns: Fitted parameter K for each gene. @@ -782,10 +782,10 @@ def compute_velocity_labeling( """Compute velocity for labeling data by: velocity = gamma / k * new - gamma * total. Args: - N: new or labeled mRNA. - R: total mRNA. - K: fitted slope k. - tau: time information. + N: New or labeled mRNA. + R: Total mRNA. + K: Fitted slope k. + tau: Time information. Returns: Computed velocity. @@ -807,14 +807,13 @@ def compute_bursting_properties( Anton J.M. Larsson et al. Genomic encoding of transcriptional burst kinetics, Nature volume 565, pages251–254(2019) Args: - M_t: the first moment of the number of total mRNA. If an array is passed, a cell-wise bursting size will be + M_t: The first moment of the number of total mRNA. If an array is passed, a cell-wise bursting size will be calculated. - phi: the reciprocal dispersion parameter. - gamma: degradation rate constant. + phi: The reciprocal dispersion parameter. + gamma: Degradation rate constant. Returns: - bs: bursting size - bf: bursting frequency + Bursting size and bursting frequency """ bs = M_t * phi bf = gamma / phi diff --git a/dynamo/estimation/csc/velocity.py b/dynamo/estimation/csc/velocity.py index 59cb7c12b..7cf27afac 100755 --- a/dynamo/estimation/csc/velocity.py +++ b/dynamo/estimation/csc/velocity.py @@ -36,13 +36,13 @@ def __init__( """Initialize the velocity class. Args: - alpha: a matrix of transcription rate. - beta: a vector of splicing rate constant for each gene. - gamma: a vector of spliced mRNA degradation rate constant for each gene. - eta: a vector of protein synthesis rate constant for each gene. - delta: a vector of protein degradation rate constant for each gene. - t: a vector of the measured time points for cells - estimation: a ss_estimation instance of the estimation class. If this not None, the parameters will be + alpha: A matrix of transcription rate. + beta: A vector of splicing rate constant for each gene. + gamma: A vector of spliced mRNA degradation rate constant for each gene. + eta: A vector of protein synthesis rate constant for each gene. + delta: A vector of protein degradation rate constant for each gene. + t: A vector of the measured time points for cells + estimation: A ss_estimation instance of the estimation class. If this not None, the parameters will be taken from this class instead of the input arguments. Returns: @@ -75,8 +75,8 @@ def vel_u( """Calculate the unspliced mRNA velocity. Args: - U: a matrix of unspliced mRNA count. Dimension: genes x cells. - repeat: whether to use average alpha or cell-wise alpha with the formula: + U: A matrix of unspliced mRNA count. Dimension: genes x cells. + repeat: Whether to use average alpha or cell-wise alpha with the formula: $a = \frac{n \gamma}{1 - e^{-\gamma t}}$. Returns: @@ -167,8 +167,8 @@ def vel_s(self, U: Union[csr_matrix, np.ndarray], S: Union[csr_matrix, np.ndarra """Calculate the unspliced mRNA velocity. Args: - U: a matrix of unspliced mRNA counts. Dimension: genes x cells. - S: a matrix of spliced mRNA counts. Dimension: genes x cells. + U: A matrix of unspliced mRNA counts. Dimension: genes x cells. + S: A matrix of spliced mRNA counts. Dimension: genes x cells. Returns: Each column of V is a velocity vector for the corresponding cell. Dimension: genes x cells. @@ -231,8 +231,8 @@ def vel_p(self, S: Union[csr_matrix, np.ndarray], P: Union[csr_matrix, np.ndarra """Calculate the protein velocity. Args: - S: a matrix of spliced mRNA counts. Dimension: genes x cells. - P: a matrix of protein counts. Dimension: genes x cells. + S: A matrix of spliced mRNA counts. Dimension: genes x cells. + P: A matrix of protein counts. Dimension: genes x cells. Returns: Each column of V is a velocity vector for the corresponding cell. Dimension: genes x cells. diff --git a/dynamo/estimation/fit_jacobian.py b/dynamo/estimation/fit_jacobian.py index 700685036..5cdf5c494 100644 --- a/dynamo/estimation/fit_jacobian.py +++ b/dynamo/estimation/fit_jacobian.py @@ -16,11 +16,11 @@ def hill_inh_func( """The inhibitory Hill equation function. Args: - x: the input transcription factor. - A: the maximum transcription rate. - K: is the ligand concentration producing half occupation. - n: the Hill coefficient. - g: the degradation constant as the offset. + x: The input transcription factor. + A: The maximum transcription rate. + K: Is the ligand concentration producing half occupation. + n: The Hill coefficient. + g: The degradation constant as the offset. Returns: The output of the Hill equation representing the inhibitory effect of the transcription factor. @@ -39,11 +39,11 @@ def hill_inh_grad( """The gradient function of the inhibitory Hill equation. Args: - x: the input transcription factor. - A: the maximum transcription rate. - K: is the ligand concentration producing half occupation. - n: the Hill coefficient. - g: the degradation constant as the offset. + x: The input transcription factor. + A: The maximum transcription rate. + K: Is the ligand concentration producing half occupation. + n: The Hill coefficient. + g: The degradation constant as the offset. Returns: The gradient of the Hill equation. @@ -62,11 +62,11 @@ def hill_act_func( """The activation Hill equation function. Args: - x: the input transcription factor. - A: the maximum transcription rate. - K: is the ligand concentration producing half occupation. - n: the Hill coefficient. - g: the degradation constant as the offset. + x: The input transcription factor. + A: The maximum transcription rate. + K: Is the ligand concentration producing half occupation. + n: The Hill coefficient. + g: The degradation constant as the offset. Returns: The output of the Hill equation representing the activation effect of the transcription factor. @@ -85,11 +85,11 @@ def hill_act_grad( """The gradient function of the activation Hill equation. Args: - x: the input transcription factor. - A: the maximum transcription rate. - K: is the ligand concentration producing half occupation. - n: the Hill coefficient. - g: the degradation constant as the offset. + x: The input transcription factor. + A: The maximum transcription rate. + K: Is the ligand concentration producing half occupation. + n: The Hill coefficient. + g: The degradation constant as the offset. Returns: The gradient of the Hill equation. @@ -108,11 +108,11 @@ def calc_mean_squared_deviation( """Calculate the mean squared deviation of the fit. Args: - func: the function to evaluate. - x_data: an array of the x data. - y_mean: an array of the mean of y data. - y_sigm: an array of the standard deviation of the y data. - weighted: whether to use weighted mean squared deviation. + func: The function to evaluate. + x_data: An array of the x data. + y_mean: An array of the mean of y data. + y_sigm: An array of the standard deviation of the y data. + weighted: Whether to use weighted mean squared deviation. Returns: The mean squared deviation. @@ -143,14 +143,14 @@ def fit_hill_grad( fitting the derivatives of inhibitory or activation Hill equations to the corresponding Jacobian elements. Args: - x_data: an array of the x data. - y_mean: an array of the mean of y data. - type: the type of the Hill equation, either `act` or `inh`. - y_sigm: an array of the standard deviation of the y data. - fix_g: the degradation constant as the offset. - n_num: the number of Hill coefficients to try. - x_tol: the tolerance of the x data. - x_shift: the shift of the x data. + x_data: An array of the x data. + y_mean: An array of the mean of y data. + type: The type of the Hill equation, either `act` or `inh`. + y_sigm: An array of the standard deviation of the y data. + fix_g: The degradation constant as the offset. + n_num: The number of Hill coefficients to try. + x_tol: The tolerance of the x data. + x_shift: The shift of the x data. Returns: The fitted parameters and the mean squared deviation. @@ -222,12 +222,12 @@ def fit_hill_inh_grad( """Fit the inhibitory Hill equation to the data. Args: - x_data: an array of the x data. - y_mean: an array of the mean of y data. - y_sigm: an array of the standard deviation of the y data. - n_num: the number of Hill coefficients to try. - x_tol: the tolerance of the x data. - x_shift: the shift of the x data. + x_data: An array of the x data. + y_mean: An array of the mean of y data. + y_sigm: An array of the standard deviation of the y data. + n_num: The number of Hill coefficients to try. + x_tol: The tolerance of the x data. + x_shift: The shift of the x data. Returns: The fitted parameters and the mean squared deviation. @@ -282,12 +282,12 @@ def fit_hill_act_grad( """Fit the activation Hill equation to the data. Args: - x_data: an array of the x data. - y_mean: an array of the mean of y data. - y_sigm: an array of the standard deviation of the y data. - n_num: the number of Hill coefficients to try. - x_tol: the tolerance of the x data. - x_shift: the shift of the x data. + x_data: An array of the x data. + y_mean: An array of the mean of y data. + y_sigm: An array of the standard deviation of the y data. + n_num: The number of Hill coefficients to try. + x_tol: The tolerance of the x data. + x_shift: The shift of the x data. Returns: The fitted parameters and the mean squared deviation. diff --git a/dynamo/estimation/tsc/estimation_kinetic.py b/dynamo/estimation/tsc/estimation_kinetic.py index 25d9caf24..d89fea0fe 100755 --- a/dynamo/estimation/tsc/estimation_kinetic.py +++ b/dynamo/estimation/tsc/estimation_kinetic.py @@ -15,8 +15,8 @@ def guestimate_alpha(x_data: np.ndarray, time: np.ndarray) -> np.ndarray: """Roughly estimate alpha for kinetics data, which is simply the averaged ratio of new RNA and labeling time. Args: - x_data: a matrix representing RNA data. - time: a matrix of labeling time. + x_data: A matrix representing RNA data. + time: A matrix of labeling time. Returns: The estimated alpha. @@ -30,8 +30,8 @@ def guestimate_gamma(x_data: np.ndarray, time: np.ndarray) -> np.ndarray: """Roughly estimate gamma0 with the assumption that time starts at 0. Args: - x_data: a matrix representing RNA data. - time: a matrix of labeling time. + x_data: A matrix representing RNA data. + time: A matrix of labeling time. Returns: The estimated gamma. @@ -44,7 +44,7 @@ def guestimate_init_cond(x_data: np.ndarray) -> np.ndarray: """Roughly estimate x0 for degradation data. Args: - x_data: a matrix representing RNA data. + x_data: A matrix representing RNA data. Returns: The estimated x0. @@ -58,8 +58,8 @@ def guestimate_p0_kinetic_chase(x_data: np.ndarray, time: np.ndarray) -> Tuple[n the average abundance of labeled RNAs across all cells belonging to the initial labeling time point. Args: - x_data: a matrix representing RNA data. - time: a matrix of labeling time. + x_data: A matrix representing RNA data. + time: A matrix of labeling time. Returns: The estimated alpha0, gamma0 and x0. @@ -118,8 +118,8 @@ def sample_p0(self, samples: int = 1, method: str = "lhs") -> np.ndarray: """Sample the initial parameters with either Latin Hypercube Sampling or random method. Args: - samples: the number of samples. - method: the sampling method. Only support "lhs" and random sampling. + samples: The number of samples. + method: The sampling method. Only support "lhs" and random sampling. Returns: The sampled array. @@ -140,7 +140,7 @@ def _lhsclassic(self, samples: int) -> np.ndarray: """Run the Latin Hypercube Sampling function. Args: - samples: the number of samples. + samples: The number of samples. Returns: The sampled data array. @@ -156,7 +156,7 @@ def get_bound(self, axis: int) -> np.ndarray: """Get the bounds of the specified axis for all parameters. Args: - axis: the index of axis. + axis: The index of axis. Returns: An array containing the bounds of the specified axis for all parameters. @@ -170,7 +170,7 @@ def normalize_data(self, X: np.ndarray) -> np.ndarray: """Perform log1p normalization on the data. Args: - X: target data to normalize. + X: Target data to normalize. Returns: The normalized data. @@ -181,8 +181,8 @@ def extract_data_from_simulator(self, t: Optional[np.ndarray] = None, **kwargs) """Extract data from the ODE simulator. Args: - t: the time information. If provided, the data will be integrated with time information. - kwargs: additional keyword arguments. + t: The time information. If provided, the data will be integrated with time information. + kwargs: Additional keyword arguments. Returns: The variable from ODE simulator. @@ -197,7 +197,7 @@ def assemble_kin_params(self, unfixed_params: np.ndarray) -> np.ndarray: """Assemble the kinetic parameters array. Args: - unfixed_params: array of unfixed parameters. + unfixed_params: Array of unfixed parameters. Returns: The assembled kinetic parameters. @@ -210,7 +210,7 @@ def assemble_x0(self, unfixed_params: np.ndarray) -> np.ndarray: """Assemble the initial conditions array. Args: - unfixed_params: array of unfixed parameters. + unfixed_params: Array of unfixed parameters. Returns: The assembled initial conditions. @@ -223,7 +223,7 @@ def set_params(self, params: np.ndarray) -> None: """Set the parameters of the simulator using assembled kinetic parameters. Args: - params: array of assembled kinetic parameters. + params: Array of assembled kinetic parameters. """ self.simulator.set_params(*self.assemble_kin_params(params)) @@ -260,11 +260,11 @@ def f_lsq( """Calculate the difference between simulated and observed data for least squares fitting. Args: - params: array of parameters for the simulation. - t: array of time values. - x_data: the input array. - method: method for integration. - normalize: whether to normalize data. + params: Array of parameters for the simulation. + t: Array of time values. + x_data: The input array. + method: Method for integration. + normalize: Whether to normalize data. Returns: Residuals representing the differences between simulated and observed data (flattened). @@ -294,16 +294,16 @@ def fit_lsq( the optimized parameters and associated cost. Args: - t: a numpy array of n time points. - x_data: an m-by-n numpy array of m species, each having n values for the n time points. - p0: initial guesses of parameters. If None, a random number is generated within the bounds. - n_p0: number of initial guesses. - bounds: lower and upper bounds for parameters. - sample_method: method used for sampling initial guesses of parameters: - `lhs`: latin hypercube sampling; - `uniform`: uniform random sampling. - method: method used for solving ODEs. See options in simulator classes. - normalize: whether to normalize values in x_data across species, so that large values do + t: A numpy array of n time points. + x_data: An m-by-n numpy array of m species, each having n values for the n time points. + p0: Initial guesses of parameters. If None, a random number is generated within the bounds. + n_p0: Number of initial guesses. + bounds: Lower and upper bounds for parameters. + sample_method: Method used for sampling initial guesses of parameters: + `lhs`: Latin hypercube sampling; + `uniform`: Uniform random sampling. + method: Method used for solving ODEs. See options in simulator classes. + normalize: Whether to normalize values in x_data across species, so that large values do not dominate the optimizer. Returns: @@ -348,7 +348,7 @@ def export_model(self, reinstantiate: bool = True) -> LinearODE: """Export the simulator model. Args: - reinstantiate: whether to reinstantiate the model class (default: True). + reinstantiate: Whether to reinstantiate the model class (default: True). Returns: Exported simulator model. @@ -416,9 +416,9 @@ def __init__(self, ranges: np.ndarray, x0: np.ndarray, simulator: LinearODE): """Initialize the Estimation_Degradation object. Args: - ranges: the lower and upper ranges of parameters. - x0: initial conditions. - simulator: instance of the Python class to solve ODEs. + ranges: The lower and upper ranges of parameters. + x0: Initial conditions. + simulator: Instance of the Python class to solve ODEs. Returns: An instance of the Estimation_Degradation class. @@ -430,7 +430,7 @@ def guestimate_init_cond(self, x_data: np.ndarray) -> np.ndarray: """Roughly estimate initial conditions for parameter estimation. Args: - x_data: a matrix representing RNA data. + x_data: A matrix representing RNA data. Returns: Estimated initial conditions. @@ -441,8 +441,8 @@ def guestimate_gamma(self, x_data: np.ndarray, time: np.ndarray) -> np.ndarray: """Roughly estimate initial conditions for parameter estimation. Args: - x_data: a matrix representing RNA data. - time: a matrix of time information. + x_data: A matrix representing RNA data. + time: A matrix of time information. Returns: Estimated gamma. @@ -453,7 +453,7 @@ def get_param(self, key: str) -> np.ndarray: """Get the estimated parameter value by key. Args: - key: the key of parameter. + key: The key of parameter. Returns: The estimated parameter value. @@ -464,7 +464,7 @@ def calc_half_life(self, key: str) -> np.ndarray: """Calculate half-life of a parameter. Args: - key: the key of parameter. + key: The key of parameter. Returns: The half-life value. @@ -503,9 +503,9 @@ def __init__( """Initialize the Estimation_DeterministicDeg object. Args: - beta: the splicing rate. - gamma: the degradation rate. - x0: the initial conditions. + beta: The splicing rate. + gamma: The degradation rate. + x0: The initial conditions. Returns: An instance of the Estimation_DeterministicDeg class. @@ -518,9 +518,9 @@ def _initialize(self, beta: np.ndarray, gamma: np.ndarray, x0: np.ndarray) -> No """Initialize the parameters to the default value. Args: - beta: the splicing rate. - gamma: the degradation rate. - x0: the initial conditions. + beta: The splicing rate. + gamma: The degradation rate. + x0: The initial conditions. """ ranges = np.zeros((2, 2)) ranges[0] = beta * np.ones(2) if np.isscalar(beta) else beta @@ -531,9 +531,9 @@ def auto_fit(self, time: np.ndarray, x_data: np.ndarray, **kwargs) -> Tuple[np.n """Estimate the parameters. Args: - time: the time information. - x_data: a matrix representing RNA data. - kwargs: the additional keyword arguments. + time: The time information. + x_data: A matrix representing RNA data. + kwargs: The additional keyword arguments. Returns: The optimized parameters and the cost. @@ -557,8 +557,8 @@ def __init__(self, gamma: Optional[np.ndarray] = None, x0: Optional[np.ndarray] """Initialize the Estimation_DeterministicDegNosp object. Args: - gamma: the degradation rate. - x0: the initial conditions. + gamma: The degradation rate. + x0: The initial conditions. Returns: An instance of the Estimation_DeterministicDegNosp class. @@ -570,8 +570,8 @@ def _initialize(self, gamma: np.ndarray, x0: np.ndarray) -> None: """Initialize the parameters to the default value. Args: - gamma: the degradation rate. - x0: the initial conditions. + gamma: The degradation rate. + x0: The initial conditions. """ ranges = gamma * np.ones(2) if np.isscalar(gamma) else gamma if np.isscalar(x0) or x0.ndim > 1: @@ -591,13 +591,13 @@ def auto_fit( """Estimate the parameters. Args: - time: the time information. - x_data: a matrix representing RNA data. - sample_method: method used for sampling initial guesses of parameters: - `lhs`: latin hypercube sampling; - `uniform`: uniform random sampling. - method: method used for solving ODEs. See options in simulator classes. - normalize: whether to normalize the data. + time: The time information. + x_data: A matrix representing RNA data. + sample_method: Method used for sampling initial guesses of parameters: + `lhs`: Latin hypercube sampling; + `uniform`: Uniform random sampling. + method: Method used for solving ODEs. See options in simulator classes. + normalize: Whether to normalize the data. Returns: The optimized parameters and the cost. @@ -635,10 +635,10 @@ def __init__( """Initialize the Estimation_MomentDeg object. Args: - beta: the splicing rate. - gamma: the degradation rate. - x0: the initial conditions. - include_cov: whether to consider covariance when estimating. + beta: The splicing rate. + gamma: The degradation rate. + x0: The initial conditions. + include_cov: Whether to consider covariance when estimating. Returns: An instance of the Estimation_MomentDeg class. @@ -652,9 +652,9 @@ def _initialize(self, beta: np.ndarray, gamma: np.ndarray, x0: np.ndarray) -> No """Initialize the parameters to the default value. Args: - beta: the splicing rate. - gamma: the degradation rate. - x0: the initial conditions. + beta: The splicing rate. + gamma: The degradation rate. + x0: The initial conditions. """ ranges = np.zeros((2, 2)) ranges[0] = beta * np.ones(2) if np.isscalar(beta) else beta @@ -686,8 +686,8 @@ def __init__(self, gamma: Optional[np.ndarray] = None, x0: Optional[np.ndarray] """Initialize the Estimation_MomentDeg object. Args: - gamma: the degradation rate. - x0: the initial conditions. + gamma: The degradation rate. + x0: The initial conditions. Returns: An instance of the Estimation_MomentDeg class. @@ -699,8 +699,8 @@ def _initialize(self, gamma: np.ndarray, x0: np.ndarray) -> None: """Initialize the parameters to the default value. Args: - gamma: the degradation rate. - x0: the initial conditions. + gamma: The degradation rate. + x0: The initial conditions. """ ranges = gamma * np.ones(2) if np.isscalar(gamma) else gamma super().__init__(ranges, x0, Moments_NoSwitchingNoSplicing()) @@ -716,13 +716,13 @@ def auto_fit( """Estimate the parameters. Args: - time: the time information. - x_data: a matrix representing RNA data. - sample_method: method used for sampling initial guesses of parameters: - `lhs`: latin hypercube sampling; - `uniform`: uniform random sampling. - method: method used for solving ODEs. See options in simulator classes. - normalize: whether to normalize the data. + time: The time information. + x_data: A matrix representing RNA data. + sample_method: Method used for sampling initial guesses of parameters: + `lhs`: Latin hypercube sampling; + `uniform`: Uniform random sampling. + method: Method used for solving ODEs. See options in simulator classes. + normalize: Whether to normalize the data. Returns: The optimized parameters and the cost. @@ -753,13 +753,13 @@ def __init__(self, a, b, alpha_a, alpha_i, beta, gamma, include_cov=True): """Initialize the Estimation_MomentKin object. Args: - a: switching rate from active promoter state to inactive promoter state. - b: switching rate from inactive promoter state to active promoter state. - alpha_a: transcription rate for active promoter. - alpha_i: transcription rate for inactive promoter. - beta: splicing rate. - gamma: degradation rate. - include_cov: whether to include the covariance when estimating. + a: Switching rate from active promoter state to inactive promoter state. + b: Switching rate from inactive promoter state to active promoter state. + alpha_a: Transcription rate for active promoter. + alpha_i: Transcription rate for inactive promoter. + beta: Splicing rate. + gamma: Degradation rate. + include_cov: Whether to include the covariance when estimating. Returns: An instance of the Estimation_MomentKin class. @@ -848,11 +848,11 @@ def __init__(self, a, b, alpha_a, alpha_i, gamma): """Initialize the Estimation_MomentKinNosp object. Args: - a: switching rate from active promoter state to inactive promoter state. - b: switching rate from inactive promoter state to active promoter state. - alpha_a: transcription rate for active promoter. - alpha_i: transcription rate for inactive promoter. - gamma: degradation rate. + a: Switching rate from active promoter state to inactive promoter state. + b: Switching rate from inactive promoter state to active promoter state. + alpha_a: Transcription rate for active promoter. + alpha_i: Transcription rate for inactive promoter. + gamma: Degradation rate. Returns: An instance of the Estimation_MomentKinNosp class. @@ -921,9 +921,9 @@ def __init__(self, alpha, gamma, x0=0): """Initialize the Estimation_DeterministicKinNosp object. Args: - alpha: transcription rate. - gamma: degradation rate. - x0: the initial condition. + alpha: Transcription rate. + gamma: Degradation rate. + x0: The initial condition. Returns: An instance of the Estimation_DeterministicKinNosp class. @@ -975,10 +975,10 @@ def __init__(self, alpha, beta, gamma, x0=np.zeros(2)): """Initialize the Estimation_DeterministicKin object. Args: - alpha: transcription rate. - beta: splicing rate. - gamma: degradation rate. - x0: the initial condition. + alpha: Transcription rate. + beta: Splicing rate. + gamma: Degradation rate. + x0: The initial condition. Returns: An instance of the Estimation_DeterministicKin class. @@ -1041,12 +1041,12 @@ def __init__(self, model1, model2, alpha=None, gamma=None, x0=None, beta=None): """Initialize the Mixture_KinDeg_NoSwitching object. Args: - model1: the first model to mix. - model2: the second model to mix. - alpha: transcription rate. - gamma: degradation rate. - x0: the initial condition. - beta: splicing rate. + model1: The first model to mix. + model2: The second model to mix. + alpha: Transcription rate. + gamma: Degradation rate. + x0: The initial condition. + beta: Splicing rate. Returns: An instance of the Mixture_KinDeg_NoSwitching class. @@ -1061,10 +1061,10 @@ def _initialize(self, alpha, gamma, x0, beta=None): """Initialize the parameters to the default value. Args: - alpha: transcription rate. - gamma: degradation rate. - x0: the initial condition. - beta: splicing rate. + alpha: Transcription rate. + gamma: Degradation rate. + x0: The initial condition. + beta: Splicing rate. """ if type(self.model1) in nosplicing_models: self.param_distributor = [[0, 2], [1, 2]] @@ -1090,8 +1090,8 @@ def normalize_deg_data(self, x_data, weight): scaling factors to ensure the data's range remains within a certain limit. Args: - x_data: a matrix representing RNA data. - weight: weight for scaling. + x_data: A matrix representing RNA data. + weight: Weight for scaling. Returns: A tuple containing the normalized degradation data and the scaling factor. @@ -1111,14 +1111,14 @@ def auto_fit(self, time, x_data, alpha_min=0.1, beta_min=50, gamma_min=10, kin_w """Estimate the parameters. Args: - time: the time information. - x_data: a matrix representing RNA data. - alpha_min: the minimum limitation on transcription rate. - beta_min: the minimum limitation on splicing rate. - gamma_min: the minimum limitation on degradation rate. - kin_weight: weight for scaling during normalization. - use_p0: whether to use initial parameters when estimating. - kwargs: the additional keyword arguments. + time: The time information. + x_data: A matrix representing RNA data. + alpha_min: The minimum limitation on transcription rate. + beta_min: The minimum limitation on splicing rate. + gamma_min: The minimum limitation on degradation rate. + kin_weight: Weight for scaling during normalization. + use_p0: Whether to use initial parameters when estimating. + kwargs: The additional keyword arguments. Returns: The optimized parameters and the cost. @@ -1163,7 +1163,7 @@ def export_model(self, reinstantiate=True): """Export the mixture model. Args: - reinstantiate: whether to reinstantiate the model. + reinstantiate: Whether to reinstantiate the model. Returns: MixtureModels or LinearODE. @@ -1221,13 +1221,13 @@ def __init__( """Initialize the Lambda_NoSwitching object. Args: - model1: the first model to mix. - model2: the second model to mix. - alpha: transcription rate. - lambd: the lambd value. - gamma: degradation rate. - x0: the initial condition. - beta: splicing rate. + model1: The first model to mix. + model2: The second model to mix. + alpha: Transcription rate. + lambd: The lambd value. + gamma: Degradation rate. + x0: The initial condition. + beta: Splicing rate. Returns: An instance of the Lambda_NoSwitching class. @@ -1242,10 +1242,10 @@ def _initialize(self, alpha, gamma, x0, beta=None): """Initialize the parameters to the default value. Args: - alpha: transcription rate. - gamma: degradation rate. - x0: the initial condition. - beta: splicing rate. + alpha: Transcription rate. + gamma: Degradation rate. + x0: The initial condition. + beta: Splicing rate. """ if type(self.model1) in nosplicing_models and type(self.model2) in nosplicing_models: self.param_keys = ["alpha", "lambda", "gamma"] @@ -1268,9 +1268,9 @@ def auto_fit(self, time, x_data, **kwargs): """Estimate the parameters. Args: - time: the time information. - x_data: a matrix representing RNA data. - kwargs: the additional keyword arguments. + time: The time information. + x_data: A matrix representing RNA data. + kwargs: The additional keyword arguments. Returns: The optimized parameters and the cost. @@ -1281,7 +1281,7 @@ def export_model(self, reinstantiate=True): """Export the mixture model. Args: - reinstantiate: whether to reinstantiate the model. + reinstantiate: Whether to reinstantiate the model. Returns: MixtureModels or LinearODE. @@ -1298,9 +1298,9 @@ def __init__(self, alpha=None, gamma=None, x0=None): """Initialize the Estimation_KineticChase object. Args: - alpha: transcription rate. - gamma: degradation rate. - x0: the initial condition. + alpha: Transcription rate. + gamma: Degradation rate. + x0: The initial condition. Returns: An instance of the Estimation_KineticChase class. @@ -1313,9 +1313,9 @@ def _initialize(self, alpha, gamma, x0): """Initialize the parameters to the default value. Args: - alpha: transcription rate. - gamma: degradation rate. - x0: the initial condition. + alpha: Transcription rate. + gamma: Degradation rate. + x0: The initial condition. """ ranges = np.zeros((2, 2)) ranges[0] = alpha * np.ones(2) if np.isscalar(alpha) else alpha @@ -1326,9 +1326,9 @@ def auto_fit(self, time, x_data, **kwargs): """Estimate the parameters. Args: - time: the time information. - x_data: a matrix representing RNA data. - kwargs: the additional keyword arguments. + time: The time information. + x_data: A matrix representing RNA data. + kwargs: The additional keyword arguments. Returns: The optimized parameters and the cost. @@ -1390,9 +1390,9 @@ def __init__(self, simulator, params=None, x0=None): """Initialize the GoodnessOfFit object. Args: - simulator: the linearODE class. - params: the parameters. - x0: the initial conditions. + simulator: The linearODE class. + params: The parameters. + x0: The initial conditions. Returns: An instance of the GoodnessOfFit class. @@ -1411,7 +1411,7 @@ def extract_data_from_simulator(self, species=None): """Extract data from the simulator's results. Args: - species: index of the species to extract. + species: Index of the species to extract. Returns: Extracted data from the simulator's results. @@ -1433,12 +1433,12 @@ def prepare_data( """Prepare data for evaluation. Args: - t: the time information. - x_data: the RNA data. - species: index of the species to consider. - method: integration method. - normalize: whether to normalize data. - reintegrate: whether to reintegrate the model. + t: The time information. + x_data: The RNA data. + species: Index of the species to consider. + method: Integration method. + normalize: Whether to normalize data. + reintegrate: Whether to reintegrate the model. """ if reintegrate: self.simulator.integrate(t, method=method) @@ -1461,8 +1461,8 @@ def normalize(self, x_data, x_model, scale=None): """Normalize data and model predictions. Args: - x_data: the RNA data. - x_model: predictions from model. + x_data: The RNA data. + x_model: Predictions from model. scale: Scaling factors for normalization. Returns: @@ -1505,7 +1505,7 @@ def calc_mean_squared_deviation(self, weighted=True): """Calculate the mean squared deviation between model predictions and observations. Args: - weighted: whether to weight the output. + weighted: Whether to weight the output. Returns: Mean squared deviation. diff --git a/dynamo/estimation/tsc/twostep.py b/dynamo/estimation/tsc/twostep.py index 426ee0b29..7f76a31cc 100644 --- a/dynamo/estimation/tsc/twostep.py +++ b/dynamo/estimation/tsc/twostep.py @@ -20,12 +20,12 @@ def fit_slope_stochastic( answer is denoted as gamma_k most of the time, which equals gamma/beta under steady state. Args: - S: a matrix of the first moments of the spliced RNA. - U: a matrix of the first moments of the unspliced RNA. - US: a matrix of the cross moments of unspliced/spliced RNA. - S2: a matrix of the second moments of spliced RNA. - perc_left: the left percentile limitation to find extreme data points. - perc_right: the right percentile limitation to find extreme data points. + S: A matrix of the first moments of the spliced RNA. + U: A matrix of the first moments of the unspliced RNA. + US: A matrix of the cross moments of unspliced/spliced RNA. + S2: A matrix of the second moments of spliced RNA. + perc_left: The left percentile limitation to find extreme data points. + perc_right: The right percentile limitation to find extreme data points. Returns: The slope, intercept, R squared and log likelihood. @@ -62,12 +62,12 @@ def fit_labeling_synthesis( """Calculate the slope of total and new RNA under steady-state assumption. Args: - new: a matrix representing new RNA. Can be expression or the first moments. - total: a matrix representing total RNA. Can be expression or the first moments. - t: a matrix of time information. - intercept: whether to perform the linear regression with intercept. - perc_left: the left percentile limitation to find extreme data points. - perc_right: the right percentile limitation to find extreme data points. + new: A matrix representing new RNA. Can be expression or the first moments. + total: A matrix representing total RNA. Can be expression or the first moments. + t: A matrix of time information. + intercept: Whether to perform the linear regression with intercept. + perc_left: The left percentile limitation to find extreme data points. + perc_right: The right percentile limitation to find extreme data points. Returns: The slope K and R squared of linear regression. @@ -90,8 +90,8 @@ def compute_gamma_synthesis( """Calculate gamma as the linear regression results of given time and log(1 - slope k). Args: - K: a matrix of the slope k. - T: a matrix of time information. + K: A matrix of the slope k. + T: A matrix of time information. Returns: The gamma and R squared of linear regression. @@ -109,10 +109,10 @@ def compute_velocity_synthesis( """Calculate the velocity of total RNA with a physical time unit: velocity = (gamma / k) N - gamma * R. Args: - N: a matrix representing new RNA. - R: a matrix representing total RNA. - gamma: a matrix of degradation rate. - t: a matrix of time information. + N: A matrix representing new RNA. + R: A matrix representing total RNA. + gamma: A matrix of degradation rate. + t: A matrix of time information. Returns: The velocity. @@ -133,10 +133,10 @@ def lin_reg_gamma_synthesis( l(t) = (1 - exp(- gamma * t)) alpha / gamma Args: - R: a matrix representing total RNA. Can be expression or the first moments. - N: a matrix representing new RNA. Can be expression or the first moments. - time: a matrix with time information. - perc_right: the percentile limitation to find extreme data points. + R: A matrix representing total RNA. Can be expression or the first moments. + N: A matrix representing new RNA. Can be expression or the first moments. + time: A matrix with time information. + perc_right: The percentile limitation to find extreme data points. Returns: Gamma, R squared, the slope k, the mean of R squared and the fitted k by time and gamma. diff --git a/dynamo/estimation/tsc/utils_kinetic.py b/dynamo/estimation/tsc/utils_kinetic.py index 05277ba3f..3b806e993 100755 --- a/dynamo/estimation/tsc/utils_kinetic.py +++ b/dynamo/estimation/tsc/utils_kinetic.py @@ -13,8 +13,8 @@ def __init__(self, n_species: int, x0: Optional[np.ndarray] = None): """Initialize the LinearODE object. Args: - n_species: the number of species. - x0: the initial condition of variable x. + n_species: The number of species. + x0: The initial condition of variable x. Returns: An instance of LinearODE. @@ -34,8 +34,8 @@ def ode_func(self, x: np.ndarray, t: np.ndarray) -> np.ndarray: """ODE functions to be implemented in the derived class such that dx=f(x, t). Args: - x: the variable. - t: the array of time. + x: The variable. + t: The array of time. Returns: The derivatives dx. @@ -47,9 +47,9 @@ def integrate(self, t: np.ndarray, x0: Optional[np.ndarray] = None, method: Opti """Integrate the ODE using the given time values. Args: - t: array of time values. - x0: array of initial conditions. - method: the method to integrate, including "matrix" and "numerical". + t: Array of time values. + x0: Array of initial conditions. + method: The method to integrate, including "matrix" and "numerical". Returns: Array containing the integrated solution over the specified time values. @@ -69,8 +69,8 @@ def integrate_numerical(self, t: np.ndarray, x0: Optional[np.ndarray] = None) -> """Numerically integrate the ODE using the given time values. Args: - t: array of time values. - x0: array of initial conditions. + t: Array of time values. + x0: Array of initial conditions. Returns: Array containing the integrated solution over the specified time values. @@ -100,8 +100,8 @@ def integrate_matrix(self, t: np.ndarray, x0: Optional[np.ndarray] = None) -> np """Integrate the system of ordinary differential equations (ODEs) using matrix exponential. Args: - t: array of time values. - x0: array of initial conditions. + t: Array of time values. + x0: Array of initial conditions. Returns: Array containing the integrated solution over the specified time values. @@ -141,8 +141,8 @@ def __init__(self, models: LinearODE, param_distributor: List): """Initialize the MixtureModels class. Args: - models: the models to mix. - param_distributor: the index to assign parameters. + models: The models to mix. + param_distributor: The index to assign parameters. Returns: An instance of MixtureModels. @@ -162,9 +162,9 @@ def integrate(self, t: np.ndarray, x0: Optional[np.ndarray] = None, method: Opti """Integrate with time values for all models. Args: - t: array of time values. - x0: array of initial conditions. - method: the method or methods to integrate, including "matrix" and "numerical". + t: Array of time values. + x0: Array of initial conditions. + method: The method or methods to integrate, including "matrix" and "numerical". """ self.x = np.zeros((len(t), np.sum(self.n_species))) for i, mdl in enumerate(self.models): @@ -178,7 +178,7 @@ def get_model_species(self, model_index: int) -> int: """Get the indices of species associated with the specified model. Args: - model_index: index of the model. + model_index: Index of the model. Returns: Array containing the indices of species associated with the specified model. @@ -199,7 +199,7 @@ def param_mixer(self, *params: Tuple) -> Tuple: """Unpack the given parameters. Args: - params: tuple of parameters. + params: Tuple of parameters. Returns: The unpacked tuple. @@ -210,7 +210,7 @@ def set_params(self, *params: Tuple) -> None: """Set parameters for all models. Args: - params: tuple of parameters. + params: Tuple of parameters. """ params = self.param_mixer(*params) for i, mdl in enumerate(self.models): @@ -230,8 +230,8 @@ def __init__(self, model1: LinearODE, model2: LinearODE): """Initialize the LambdaModels_NoSwitching class. Args: - model1: the first model to mix. - model2: the second model to mix. + model1: The first model to mix. + model2: The second model to mix. Returns: An instance of LambdaModels_NoSwitching. @@ -249,7 +249,7 @@ def param_mixer(self, *params) -> np.ndarray: """Set parameters for all models. Args: - params: tuple of parameters. + params: Tuple of parameters. """ lam = params[1] alp_1 = params[0] * lam @@ -274,13 +274,13 @@ def __init__( """Initialize the Moments object. Args: - a: switching rate from active promoter state to inactive promoter state. - b: switching rate from inactive promoter state to active promoter state. - alpha_a: transcription rate for active promoter. - alpha_i: transcription rate for inactive promoter. - beta: splicing rate. - gamma: degradation rate. - x0: the initial conditions. + a: Switching rate from active promoter state to inactive promoter state. + b: Switching rate from inactive promoter state to active promoter state. + alpha_a: Transcription rate for active promoter. + alpha_i: Transcription rate for inactive promoter. + beta: Splicing rate. + gamma: Degradation rate. + x0: The initial conditions. Returns: An instance of Moments. @@ -312,8 +312,8 @@ def ode_func(self, x: np.ndarray, t: np.ndarray) -> np.ndarray: The second moments is calculated from the variance and covariance of variable u and x. Args: - x: the variable. - t: the array of time. + x: The variable. + t: The array of time. Returns: The derivatives dx. @@ -344,8 +344,8 @@ def fbar(self, x_a: np.ndarray, x_i: np.ndarray) -> np.ndarray: """Calculate the count of a variable by averaging active and inactive states. Args: - x_a: the variable x under the active state. - x_i: the variable x under the inactive state. + x_a: The variable x under the active state. + x_i: The variable x under the inactive state. Returns: The count of variable x. @@ -364,12 +364,12 @@ def set_params( """Set the parameters. Args: - a: switching rate from active promoter state to inactive promoter state. - b: switching rate from inactive promoter state to active promoter state. - alpha_a: transcription rate for active promoter. - alpha_i: transcription rate for inactive promoter. - beta: splicing rate. - gamma: degradation rate. + a: Switching rate from active promoter state to inactive promoter state. + b: Switching rate from inactive promoter state to active promoter state. + alpha_a: Transcription rate for active promoter. + alpha_i: Transcription rate for inactive promoter. + beta: Splicing rate. + gamma: Degradation rate. """ self.a = a self.b = b @@ -539,12 +539,12 @@ def __init__( """Initialize the Moments_Nosplicing object. Args: - a: switching rate from active promoter state to inactive promoter state. - b: switching rate from inactive promoter state to active promoter state. - alpha_a: transcription rate for active promoter. - alpha_i: transcription rate for inactive promoter. - gamma: degradation rate. - x0: the initial conditions. + a: Switching rate from active promoter state to inactive promoter state. + b: Switching rate from inactive promoter state to active promoter state. + alpha_a: Transcription rate for active promoter. + alpha_i: Transcription rate for inactive promoter. + gamma: Degradation rate. + x0: The initial conditions. Returns: An instance of Moments_Nosplicing. @@ -567,8 +567,8 @@ def ode_func(self, x: np.ndarray, t: np.ndarray) -> np.ndarray: """ODE functions to solve. Ignore the splicing part in the base class. Args: - x: the variable. - t: the array of time. + x: The variable. + t: The array of time. Returns: The derivatives dx. @@ -594,8 +594,8 @@ def fbar(self, x_a: np.ndarray, x_i: np.ndarray) -> np.ndarray: """Calculate the count of a variable by averaging active and inactive states. Args: - x_a: the variable x under the active state. - x_i: the variable x under the inactive state. + x_a: The variable x under the active state. + x_i: The variable x under the inactive state. Returns: The count of variable x. @@ -606,11 +606,11 @@ def set_params(self, a: np.ndarray, b: np.ndarray, alpha_a: np.ndarray, alpha_i: """Set the parameters. Args: - a: switching rate from active promoter state to inactive promoter state. - b: switching rate from inactive promoter state to active promoter state. - alpha_a: transcription rate for active promoter. - alpha_i: transcription rate for inactive promoter. - gamma: degradation rate. + a: Switching rate from active promoter state to inactive promoter state. + b: Switching rate from inactive promoter state to active promoter state. + alpha_a: Transcription rate for active promoter. + alpha_i: Transcription rate for inactive promoter. + gamma: Degradation rate. """ self.a = a self.b = b @@ -696,10 +696,10 @@ def __init__( """Initialize the Moments_NoSwitching object. Args: - alpha: transcription rate. - beta: splicing rate. - gamma: degradation rate. - x0: the initial conditions. + alpha: Transcription rate. + beta: Splicing rate. + gamma: Degradation rate. + x0: The initial conditions. Returns: An instance of Moments_NoSwitching. @@ -724,8 +724,8 @@ def ode_func(self, x: np.ndarray, t: np.ndarray) -> np.ndarray: """ODE functions to solve. Ignore the switching part in the base class. Args: - x: the variable. - t: the array of time. + x: The variable. + t: The array of time. Returns: The derivatives dx. @@ -751,9 +751,9 @@ def set_params(self, alpha: np.ndarray, beta: np.ndarray, gamma: np.ndarray) -> """Set the parameters. Args: - alpha: transcription rate. - beta: splicing rate. - gamma: degradation rate. + alpha: Transcription rate. + beta: Splicing rate. + gamma: Degradation rate. """ self.al = alpha self.be = beta @@ -884,9 +884,9 @@ def __init__( """Initialize the Moments_NoSwitchingNoSplicing object. Args: - alpha: transcription rate. - gamma: degradation rate. - x0: the initial conditions. + alpha: Transcription rate. + gamma: Degradation rate. + x0: The initial conditions. Returns: An instance of Moments_NoSwitchingNoSplicing. @@ -908,8 +908,8 @@ def ode_func(self, x: np.ndarray, t: np.ndarray) -> np.ndarray: """ODE functions to solve. Both splicing and switching part in the base class are ignored. Args: - x: the variable. - t: the array of time. + x: The variable. + t: The array of time. Returns: The derivatives dx. @@ -931,8 +931,8 @@ def set_params(self, alpha: np.ndarray, gamma: np.ndarray) -> None: """Set the parameters. Args: - alpha: transcription rate. - gamma: degradation rate. + alpha: Transcription rate. + gamma: Degradation rate. """ self.al = alpha self.ga = gamma @@ -1007,10 +1007,10 @@ def __init__( """Initialize the Deterministic object. Args: - alpha: the transcription rate. - beta: the splicing rate. - gamma: the degradation rate. - x0: the inital conditions. + alpha: The transcription rate. + beta: The splicing rate. + gamma: The degradation rate. + x0: The initial conditions. Returns: An instance of Deterministic. @@ -1037,8 +1037,8 @@ def ode_func(self, x: np.ndarray, t: np.ndarray) -> np.ndarray: dx[v] = beta * u - gamma * s Args: - x: the x variable. - t: the time information. + x: The x variable. + t: The time information. Returns: An array containing the ODEs' output. @@ -1059,9 +1059,9 @@ def set_params(self, alpha: np.ndarray, beta: np.ndarray, gamma: np.ndarray) -> """Set the parameters. Args: - alpha: the transcription rate. - beta: the splicing rate. - gamma: the degradation rate. + alpha: The transcription rate. + beta: The splicing rate. + gamma: The degradation rate. """ self.al = alpha self.be = beta @@ -1074,9 +1074,9 @@ def integrate(self, t: np.ndarray, x0: Optional[np.ndarray] = None, method: str """Integrate the ODE using the given time values. Args: - t: array of time values. - x0: array of initial conditions. - method: the method to integrate, including "matrix" and "numerical". + t: Array of time values. + x0: Array of initial conditions. + method: The method to integrate, including "matrix" and "numerical". """ method = self.default_method if method is None else method if method == "analytical": @@ -1114,8 +1114,8 @@ def integrate_analytical(self, t: np.ndarray, x0: Optional[np.ndarray] = None) - """Integrate the odes with the analytical solution. Args: - t: the time information. - x0: the initial conditions. + t: The time information. + x0: The initial conditions. Returns: The solution of unspliced and splcied mRNA wrapped in an array. @@ -1140,9 +1140,9 @@ def __init__( """Initialize the Deterministic_NoSplicing object. Args: - alpha: the transcription rate. - gamma: the degradation rate. - x0: the initial conditions. + alpha: The transcription rate. + gamma: The degradation rate. + x0: The initial conditions. Returns: An instance of Deterministic_NoSplicing. @@ -1167,8 +1167,8 @@ def ode_func(self, x: np.ndarray, t: np.ndarray) -> np.ndarray: dx[u] = alpha - gamma * u Args: - x: the x variable. - t: the time information. + x: The x variable. + t: The time information. Returns: An array containing the ODEs' output. @@ -1187,8 +1187,8 @@ def set_params(self, alpha: np.ndarray, gamma: np.ndarray) -> np.ndarray: """Set the parameters. Args: - alpha: the transcription rate. - gamma: the degradation rate. + alpha: The transcription rate. + gamma: The degradation rate. """ self.al = alpha self.ga = gamma @@ -1200,9 +1200,9 @@ def integrate(self, t: np.ndarray, x0: Optional[np.ndarray] = None, method: str """Integrate the ODE using the given time values. Args: - t: array of time values. - x0: array of initial conditions. - method: the method to integrate, including "matrix" and "numerical". + t: Array of time values. + x0: Array of initial conditions. + method: The method to integrate, including "matrix" and "numerical". """ method = self.default_method if method is None else method if method == "analytical": @@ -1235,8 +1235,8 @@ def integrate_analytical(self, t: np.ndarray, x0: Optional[np.ndarray] = None) - """Integrate the odes with the analytical solution. Args: - t: the time information. - x0: the initial conditions. + t: The time information. + x0: The initial conditions. Returns: The solution of unspliced mRNA as an array. From 9ab65e0434759b02402448938a3c1b23c61bcfc5 Mon Sep 17 00:00:00 2001 From: sichao Date: Fri, 3 Nov 2023 14:48:37 -0400 Subject: [PATCH 35/36] add missing typing in estimation_linetic --- dynamo/estimation/tsc/estimation_kinetic.py | 375 +++++++++++++++----- 1 file changed, 282 insertions(+), 93 deletions(-) diff --git a/dynamo/estimation/tsc/estimation_kinetic.py b/dynamo/estimation/tsc/estimation_kinetic.py index d89fea0fe..441b0f548 100755 --- a/dynamo/estimation/tsc/estimation_kinetic.py +++ b/dynamo/estimation/tsc/estimation_kinetic.py @@ -3,6 +3,7 @@ import numpy as np from scipy.optimize import least_squares +from scipy.sparse import csr_matrix from scipy.stats import chi2 from ...dynamo_logger import main_warning @@ -749,7 +750,16 @@ class Estimation_MomentKin(kinetic_estimation): Order of species: , , , , """ - def __init__(self, a, b, alpha_a, alpha_i, beta, gamma, include_cov=True): + def __init__( + self, + a: np.ndarray, + b: np.ndarray, + alpha_a: np.ndarray, + alpha_i: np.ndarray, + beta: np.ndarray, + gamma: np.ndarray, + include_cov: bool = True, + ): """Initialize the Estimation_MomentKin object. Args: @@ -775,8 +785,12 @@ def __init__(self, a, b, alpha_a, alpha_i, beta, gamma, include_cov=True): super().__init__(ranges, np.zeros((7, 2)), Moments()) self.include_cov = include_cov - def extract_data_from_simulator(self): - """Get corresponding data from the LinearODE class.""" + def extract_data_from_simulator(self) -> np.ndarray: + """Get corresponding data from the LinearODE class. + + Returns: + The variable from ODE simulator as an array. + """ if self.include_cov: ret = np.zeros((5, len(self.simulator.t))) ret[0] = self.simulator.get_nu() @@ -792,36 +806,64 @@ def extract_data_from_simulator(self): ret[3] = self.simulator.x[:, self.simulator.xx] return ret - def get_alpha_a(self): - """Get the transcription rate for active promoter.""" + def get_alpha_a(self) -> np.ndarray: + """Get the transcription rate for active promoter. + + Returns: + The transcription rate for active promoter. + """ return self.popt[2] - def get_alpha_i(self): - """Get the transcription rate for inactive promoter.""" + def get_alpha_i(self) -> np.ndarray: + """Get the transcription rate for inactive promoter. + + Returns: + The transcription rate for inactive promoter. + """ return self.popt[3] - def get_alpha(self): - """Get all transcription rates.""" + def get_alpha(self) -> np.ndarray: + """Get all transcription rates. + + Returns: + All transcription rates. + """ alpha = self.simulator.fbar(self.get_alpha_a(), self.get_alpha_i()) return alpha - def get_beta(self): - """Get the splicing rate.""" + def get_beta(self) -> np.ndarray: + """Get the splicing rate. + + Returns: + The splicing rate. + """ return self.popt[4] - def get_gamma(self): - """Get the degradation rate.""" + def get_gamma(self) -> np.ndarray: + """Get the degradation rate. + + Returns: + The degradation rate. + """ return self.popt[5] - def calc_spl_half_life(self): - """Calculate the half life of splicing.""" + def calc_spl_half_life(self) -> np.ndarray: + """Calculate the half life of splicing. + + Returns: + The half life of splicing. + """ return np.log(2) / self.get_beta() - def calc_deg_half_life(self): - """Calculate the half life of degradation.""" + def calc_deg_half_life(self) -> np.ndarray: + """Calculate the half life of degradation. + + Returns: + The half life of degradation. + """ return np.log(2) / self.get_gamma() - def export_dictionary(self): + def export_dictionary(self) -> Dict: """Export parameter estimation results as a dictionary. Returns: @@ -844,7 +886,14 @@ class Estimation_MomentKinNosp(kinetic_estimation): Order of species: , """ - def __init__(self, a, b, alpha_a, alpha_i, gamma): + def __init__( + self, + a: np.ndarray, + b: np.ndarray, + alpha_a: np.ndarray, + alpha_i: np.ndarray, + gamma: np.ndarray, + ): """Initialize the Estimation_MomentKinNosp object. Args: @@ -866,35 +915,59 @@ def __init__(self, a, b, alpha_a, alpha_i, gamma): ranges[4] = gamma * np.ones(2) if np.isscalar(gamma) else gamma super().__init__(ranges, np.zeros((3, 2)), Moments_Nosplicing()) - def extract_data_from_simulator(self): - """Get corresponding data from the LinearODE class.""" + def extract_data_from_simulator(self) -> np.ndarray: + """Get corresponding data from the LinearODE class. + + Returns: + The variable from ODE simulator as an array. + """ ret = np.zeros((2, len(self.simulator.t))) ret[0] = self.simulator.get_nu() ret[1] = self.simulator.x[:, self.simulator.uu] return ret - def get_alpha_a(self): - """Get the transcription rate for active promoter.""" + def get_alpha_a(self) -> np.ndarray: + """Get the transcription rate for active promoter. + + Returns: + The transcription rate for active promoter. + """ return self.popt[2] - def get_alpha_i(self): - """Get the transcription rate for inactive promoter.""" + def get_alpha_i(self) -> np.ndarray: + """Get the transcription rate for inactive promoter. + + Returns: + The transcription rate for inactive promoter. + """ return self.popt[3] - def get_alpha(self): - """Get all transcription rates.""" + def get_alpha(self) -> np.ndarray: + """Get all transcription rates. + + Returns: + All transcription rates. + """ alpha = self.simulator.fbar(self.get_alpha_a().self.get_alpha_i()) return alpha - def get_gamma(self): - """Get the degradation rate.""" + def get_gamma(self) -> np.ndarray: + """Get the degradation rate. + + Returns: + The degradation rate. + """ return self.popt[4] - def calc_deg_half_life(self): - """Calculate the half life of degradation.""" + def calc_deg_half_life(self) -> np.ndarray: + """Calculate the half life of degradation. + + Returns: + The half life of degradation. + """ return np.log(2) / self.get_gamma() - def export_dictionary(self): + def export_dictionary(self) -> Dict: """Export parameter estimation results as a dictionary. Returns: @@ -917,7 +990,12 @@ class Estimation_DeterministicKinNosp(kinetic_estimation): Order of species: , """ - def __init__(self, alpha, gamma, x0=0): + def __init__( + self, + alpha: np.ndarray, + gamma: np.ndarray, + x0: Union[int, np.ndarray] = 0, + ): """Initialize the Estimation_DeterministicKinNosp object. Args: @@ -936,19 +1014,27 @@ def __init__(self, alpha, gamma, x0=0): x0 = np.ones((1, 2)) * x0 super().__init__(ranges, x0, Deterministic_NoSplicing()) - def get_alpha(self): - """Get the transcription rate.""" + def get_alpha(self) -> np.ndarray: + """Get the transcription rate. + + Returns: + The transcription rate. + """ return self.popt[0] - def get_gamma(self): - """Get the degradation rate.""" + def get_gamma(self) -> np.ndarray: + """Get the degradation rate. + + Returns: + The degradation rate. + """ return self.popt[1] - def calc_half_life(self, key): + def calc_half_life(self, key: str) -> np.ndarray: """Calculate the half life.""" return np.log(2) / self.get_param(key) - def export_dictionary(self): + def export_dictionary(self) -> Dict: """Export parameter estimation results as a dictionary. Returns: @@ -971,7 +1057,13 @@ class Estimation_DeterministicKin(kinetic_estimation): Order of species: , """ - def __init__(self, alpha, beta, gamma, x0=np.zeros(2)): + def __init__( + self, + alpha: np.ndarray, + beta: np.ndarray, + gamma: np.ndarray, + x0: Union[int, np.ndarray] = np.zeros(2), + ): """Initialize the Estimation_DeterministicKin object. Args: @@ -994,27 +1086,47 @@ def __init__(self, alpha, beta, gamma, x0=np.zeros(2)): super().__init__(ranges, x0, Deterministic()) - def get_alpha(self): - """Get the transcription rate.""" + def get_alpha(self) -> np.ndarray: + """Get the transcription rate. + + Returns: + The transcription rate. + """ return self.popt[0] - def get_beta(self): - """Get the splicing rate.""" + def get_beta(self) -> np.ndarray: + """Get the splicing rate. + + Returns: + The splicing rate. + """ return self.popt[1] - def get_gamma(self): - """Get the degradation rate.""" + def get_gamma(self) -> np.ndarray: + """Get the degradation rate. + + Returns: + The degradation rate. + """ return self.popt[2] - def calc_spl_half_life(self): - """Calculate the half life of splicing.""" + def calc_spl_half_life(self) -> np.ndarray: + """Calculate the half life of splicing. + + Returns: + The half life of splicing. + """ return np.log(2) / self.get_beta() - def calc_deg_half_life(self): - """Calculate the half life of degradation.""" + def calc_deg_half_life(self) -> np.ndarray: + """Calculate the half life of degradation. + + Returns: + The half life of degradation. + """ return np.log(2) / self.get_gamma() - def export_dictionary(self): + def export_dictionary(self) -> Dict: """Export parameter estimation results as a dictionary. Returns: @@ -1037,7 +1149,15 @@ class Mixture_KinDeg_NoSwitching(kinetic_estimation): If beta is None, it is assumed that the data does not have the splicing process. """ - def __init__(self, model1, model2, alpha=None, gamma=None, x0=None, beta=None): + def __init__( + self, + model1: LinearODE, + model2: LinearODE, + alpha: Optional[np.ndarray] = None, + gamma: Optional[np.ndarray] = None, + x0: Optional[Union[int, np.ndarray]] = None, + beta: Optional[np.ndarray] = None, + ): """Initialize the Mixture_KinDeg_NoSwitching object. Args: @@ -1057,7 +1177,13 @@ def __init__(self, model1, model2, alpha=None, gamma=None, x0=None, beta=None): if alpha is not None and gamma is not None: self._initialize(alpha, gamma, x0, beta) - def _initialize(self, alpha, gamma, x0, beta=None): + def _initialize( + self, + alpha: np.ndarray, + gamma: np.ndarray, + x0: Union[int, np.ndarray], + beta: Optional[np.ndarray] = None, + ): """Initialize the parameters to the default value. Args: @@ -1065,6 +1191,9 @@ def _initialize(self, alpha, gamma, x0, beta=None): gamma: Degradation rate. x0: The initial condition. beta: Splicing rate. + + Returns: + An instance of the Mixture_KinDeg_NoSwitching class. """ if type(self.model1) in nosplicing_models: self.param_distributor = [[0, 2], [1, 2]] @@ -1085,7 +1214,9 @@ def _initialize(self, alpha, gamma, x0, beta=None): x0_ = np.vstack((np.zeros((self.model1.n_species, 2)), x0)) super().__init__(ranges, x0_, model) - def normalize_deg_data(self, x_data, weight): + def normalize_deg_data( + self, x_data: Union[csr_matrix, np.ndarray], weight: Union[float, int] + ) -> Tuple[np.ndarray, np.ndarray]: """Normalize the degradation data while preserving the relative proportions between species. It calculates scaling factors to ensure the data's range remains within a certain limit. @@ -1107,7 +1238,17 @@ def normalize_deg_data(self, x_data, weight): return x_data_norm, scale - def auto_fit(self, time, x_data, alpha_min=0.1, beta_min=50, gamma_min=10, kin_weight=2, use_p0=True, **kwargs): + def auto_fit( + self, + time: np.ndarray, + x_data: Union[csr_matrix, np.ndarray], + alpha_min: Union[float, int] = 0.1, + beta_min: Union[float, int] = 50, + gamma_min: Union[float, int] = 10, + kin_weight: Union[float, int] = 2, + use_p0: bool = True, + **kwargs, + ) -> Tuple[np.ndarray, float]: """Estimate the parameters. Args: @@ -1159,7 +1300,7 @@ def auto_fit(self, time, x_data, alpha_min=0.1, beta_min=50, gamma_min=10, kin_w popt, cost = self.fit_lsq(time, x_data_norm, **kwargs) return popt, cost - def export_model(self, reinstantiate=True): + def export_model(self, reinstantiate: bool = True) -> Union[MixtureModels, LinearODE]: """Export the mixture model. Args: @@ -1173,7 +1314,7 @@ def export_model(self, reinstantiate=True): else: return self.simulator - def export_x0(self): + def export_x0(self) -> Optional[np.ndarray]: """Export optimized initial conditions for the mixture of models analysis. Returns: @@ -1183,7 +1324,7 @@ def export_x0(self): x[self.model1.n_species :] *= self.scale return x - def export_dictionary(self): + def export_dictionary(self) -> Dict: """Export parameter estimation results as a dictionary. Returns: @@ -1210,13 +1351,13 @@ class Lambda_NoSwitching(Mixture_KinDeg_NoSwitching): def __init__( self, - model1, - model2, - alpha=None, - lambd=None, - gamma=None, - x0=None, - beta=None, + model1: LinearODE, + model2: LinearODE, + alpha: Optional[np.ndarray] = None, + lambd: Optional[np.ndarray] = None, + gamma: Optional[np.ndarray] = None, + x0: Optional[Union[int, np.ndarray]] = None, + beta: Optional[np.ndarray] = None, ): """Initialize the Lambda_NoSwitching object. @@ -1238,7 +1379,13 @@ def __init__( if alpha is not None and gamma is not None: self._initialize(alpha, gamma, x0, beta) - def _initialize(self, alpha, gamma, x0, beta=None): + def _initialize( + self, + alpha: np.ndarray, + gamma: np.ndarray, + x0: Union[int, np.ndarray], + beta: Optional[np.ndarray] = None, + ): """Initialize the parameters to the default value. Args: @@ -1246,6 +1393,9 @@ def _initialize(self, alpha, gamma, x0, beta=None): gamma: Degradation rate. x0: The initial condition. beta: Splicing rate. + + Returns: + An instance of the Lambda_NoSwitching class. """ if type(self.model1) in nosplicing_models and type(self.model2) in nosplicing_models: self.param_keys = ["alpha", "lambda", "gamma"] @@ -1264,7 +1414,7 @@ def _initialize(self, alpha, gamma, x0, beta=None): x0_ = np.vstack((np.zeros((self.model1.n_species, 2)), x0)) super(Mixture_KinDeg_NoSwitching, self).__init__(ranges, x0_, model) - def auto_fit(self, time, x_data, **kwargs): + def auto_fit(self, time: np.ndarray, x_data: Union[csr_matrix, np.ndarray], **kwargs) -> Tuple[np.ndarray, float]: """Estimate the parameters. Args: @@ -1277,7 +1427,7 @@ def auto_fit(self, time, x_data, **kwargs): """ return super().auto_fit(time, x_data, kin_weight=None, use_p0=False, **kwargs) - def export_model(self, reinstantiate=True): + def export_model(self, reinstantiate: bool = True) -> Union[LambdaModels_NoSwitching, LinearODE]: """Export the mixture model. Args: @@ -1294,7 +1444,12 @@ def export_model(self, reinstantiate=True): class Estimation_KineticChase(kinetic_estimation): """An estimation class for kinetic chase experiment.""" - def __init__(self, alpha=None, gamma=None, x0=None): + def __init__( + self, + alpha: Optional[np.ndarray] = None, + gamma: Optional[np.ndarray] = None, + x0: Optional[Union[int, np.ndarray]] = None, + ): """Initialize the Estimation_KineticChase object. Args: @@ -1309,20 +1464,28 @@ def __init__(self, alpha=None, gamma=None, x0=None): if alpha is not None and gamma is not None and x0 is not None: self._initialize(alpha, gamma, x0) - def _initialize(self, alpha, gamma, x0): + def _initialize( + self, + alpha: np.ndarray, + gamma: np.ndarray, + x0: Union[int, np.ndarray], + ): """Initialize the parameters to the default value. Args: alpha: Transcription rate. gamma: Degradation rate. x0: The initial condition. + + Returns: + An instance of the Estimation_KineticChase class. """ ranges = np.zeros((2, 2)) ranges[0] = alpha * np.ones(2) if np.isscalar(alpha) else alpha ranges[1] = gamma * np.ones(2) if np.isscalar(gamma) else gamma super().__init__(ranges, np.atleast_2d(x0), KineticChase()) - def auto_fit(self, time, x_data, **kwargs): + def auto_fit(self, time: np.ndarray, x_data: Union[csr_matrix, np.ndarray], **kwargs) -> Tuple[np.ndarray, float]: """Estimate the parameters. Args: @@ -1346,23 +1509,39 @@ def auto_fit(self, time, x_data, **kwargs): popt, cost = self.fit_lsq(time, x_data, p0=np.hstack((al0, ga0, x0)), normalize=False, **kwargs) return popt, cost - def get_param(self, key): - """Get corresponding parameter according to the key.""" + def get_param(self, key: str) -> np.ndarray: + """Get corresponding parameter according to the key. + + Returns: + The corresponding parameter. + """ return self.popt[np.where(self.kin_param_keys == key)[0][0]] - def get_alpha(self): - """Get the transcription rate.""" + def get_alpha(self) -> np.ndarray: + """Get the transcription rate. + + Returns: + The transcription rate. + """ return self.popt[0] - def get_gamma(self): - """Get the degradation rate.""" + def get_gamma(self) -> np.ndarray: + """Get the degradation rate. + + Returns: + The degradation rate. + """ return self.popt[1] - def calc_half_life(self, key): - """Calculate the half life.""" + def calc_half_life(self, key: str) -> np.ndarray: + """Calculate the half life. + + Returns: + The half life. + """ return np.log(2) / self.get_param(key) - def export_dictionary(self): + def export_dictionary(self) -> Dict: """Export parameter estimation results as a dictionary. Returns: @@ -1386,7 +1565,12 @@ class GoodnessOfFit: This class provides methods for assessing the quality of predictions, using various metrics including Gaussian likelihood, Gaussian log-likelihood, and mean squared deviation. """ - def __init__(self, simulator, params=None, x0=None): + def __init__( + self, + simulator: LinearODE, + params: Optional[Tuple] = None, + x0: Optional[Union[int, np.ndarray]] = None, + ): """Initialize the GoodnessOfFit object. Args: @@ -1407,7 +1591,7 @@ def __init__(self, simulator, params=None, x0=None): self.sigm = None self.pred = None - def extract_data_from_simulator(self, species=None): + def extract_data_from_simulator(self, species: Optional[int] = None) -> np.ndarray: """Extract data from the simulator's results. Args: @@ -1423,13 +1607,13 @@ def extract_data_from_simulator(self, species=None): def prepare_data( self, - t, - x_data, - species=None, - method=None, - normalize=True, - reintegrate=True, - ): + t: np.ndarray, + x_data: Union[csr_matrix, np.ndarray], + species: Optional[int] = None, + method: Optional[str] = None, + normalize: bool = True, + reintegrate: bool = True, + ) -> None: """Prepare data for evaluation. Args: @@ -1457,7 +1641,12 @@ def prepare_data( self.sigm = strat_mom(x_data_norm.T, t, np.std) self.pred = strat_mom(x_model_norm.T, t, np.mean) - def normalize(self, x_data, x_model, scale=None): + def normalize( + self, + x_data: Union[csr_matrix, np.ndarray], + x_model: np.ndarray, + scale: Optional[Union[float, int, np.ndarray]] = None, + ) -> Tuple[np.ndarray, np.ndarray]: """Normalize data and model predictions. Args: @@ -1473,7 +1662,7 @@ def normalize(self, x_data, x_model, scale=None): x_model_norm = (x_model.T / scale).T return x_data_norm, x_model_norm - def calc_gaussian_likelihood(self): + def calc_gaussian_likelihood(self) -> float: """ Calculate the Gaussian likelihood between model predictions and observations. Returns: @@ -1487,7 +1676,7 @@ def calc_gaussian_likelihood(self): ret = 1 / (np.sqrt((2 * np.pi) ** len(err)) * np.prod(sig)) * np.exp(-0.5 * (err).dot(err)) return ret - def calc_gaussian_loglikelihood(self): + def calc_gaussian_loglikelihood(self) -> float: """Calculate the Gaussian log-likelihood between model predictions and observations. Returns: @@ -1501,7 +1690,7 @@ def calc_gaussian_loglikelihood(self): ret = -len(err) / 2 * np.log(2 * np.pi) - np.sum(np.log(sig)) - 0.5 * err.dot(err) return ret - def calc_mean_squared_deviation(self, weighted=True): + def calc_mean_squared_deviation(self, weighted: bool = True) -> float: """Calculate the mean squared deviation between model predictions and observations. Args: From 8d09107fc10ac7596e095d19350e73556365ea03 Mon Sep 17 00:00:00 2001 From: sichao Date: Tue, 14 Nov 2023 11:46:47 -0500 Subject: [PATCH 36/36] fix typo --- dynamo/estimation/tsc/utils_kinetic.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/dynamo/estimation/tsc/utils_kinetic.py b/dynamo/estimation/tsc/utils_kinetic.py index 3b806e993..9dff4d051 100755 --- a/dynamo/estimation/tsc/utils_kinetic.py +++ b/dynamo/estimation/tsc/utils_kinetic.py @@ -1118,7 +1118,7 @@ def integrate_analytical(self, t: np.ndarray, x0: Optional[np.ndarray] = None) - x0: The initial conditions. Returns: - The solution of unspliced and splcied mRNA wrapped in an array. + The solution of unspliced and spliced mRNA wrapped in an array. """ if x0 is None: x0 = self.x0