diff --git a/dynamo/simulation/Gillespie.py b/dynamo/simulation/Gillespie.py index 3224158d1..995533af7 100755 --- a/dynamo/simulation/Gillespie.py +++ b/dynamo/simulation/Gillespie.py @@ -60,9 +60,10 @@ def Gillespie( """ gene_num, species_num = C0.shape[0:2] - adata_no_splicing, P = None, None + adata_no_splicing, P, layers_no_splicing = None, None, None if method == "basic": + gene_num = 2 if t_eval is None: steps = (t_span[1] - t_span[0]) // dt # // int; %% remainder t_eval = np.linspace(t_span[0], t_span[1], steps) @@ -75,10 +76,10 @@ def Gillespie( si=si, be=be, ga=ga, - C0=np.zeros((5, 1)), + C0=[np.zeros((1, 5))] * 2, t_span=[0, 50], n_traj=1, - t_eval=None, + t_eval=t_eval, report=verbose, ) # unfinished, no need to interpolate now. uu, ul, su, sl = [np.transpose(trajs_C[:, :, i + 1, :].reshape((gene_num, -1))) for i in range(4)] @@ -108,20 +109,22 @@ def Gillespie( } ) obs.set_index("Cell_name", inplace=True) + true_beta, true_gamma, delta, eta = be, ga, None, None elif method == "simulate_2bifurgenes": gene_num = 2 + param_dict = { + "a": [20, 20], + "b": [20, 20], + "S": [20, 20], + "K": [20, 20], + "m": [3, 3], + "n": [3, 3], + "beta": [1, 1], + "gamma": [1, 1], + } _, trajs_C = simulate_2bifurgenes( - a1=20, - b1=20, - a2=20, - b2=20, - K=20, - n=3, - be1=1, - ga1=1, - be2=1, - ga2=1, C0=np.zeros(4), + param_dict=param_dict, t_span=t_span, n_traj=n_traj, report=verbose, @@ -150,6 +153,7 @@ def Gillespie( } ) obs.set_index("Cell_name", inplace=True) + true_beta, true_gamma, delta, eta = param_dict["beta"], param_dict["gamma"], None, None elif method == "differentiation": gene_num = 2 @@ -413,6 +417,7 @@ def Gillespie( } ) obs.set_index("cell_name", inplace=True) + true_beta, true_gamma = [beta, beta], [gamma, gamma] elif method == "oscillation": gene_num = 2 @@ -658,6 +663,7 @@ def Gillespie( } ) obs.set_index("cell_name", inplace=True) + true_beta, true_gamma = [beta, beta], [gamma, gamma] else: raise Exception("method not implemented!") # anadata: observation x variable (cells x genes) @@ -668,8 +674,8 @@ def Gillespie( var = pd.DataFrame( { "gene_short_name": ["gene_%d" % (i) for i in range(gene_num)], - "true_beta": [beta, beta], - "true_gamma": [gamma, gamma], + "true_beta": true_beta, + "true_gamma": true_gamma, "true_eta": [eta, eta], "true_delta": [delta, delta], } @@ -682,17 +688,21 @@ def Gillespie( var.copy(), layers=layers.copy(), ) - adata_no_splicing = anndata.AnnData( - scipy.sparse.csc_matrix(E.astype(int)).copy(), - obs.copy(), - var.copy(), - layers=layers_no_splicing.copy(), - ) if P is not None: adata.obsm["protein"] = P - adata_no_splicing.obsm["protein"] = P # remove cells that has no expression adata = adata[np.array(adata.X.sum(1)).flatten() > 0, :] - adata_no_splicing = adata_no_splicing[np.array(adata_no_splicing.X.sum(1)).flatten() > 0, :] + + if layers_no_splicing is not None: + adata_no_splicing = anndata.AnnData( + scipy.sparse.csc_matrix(E.astype(int)).copy(), + obs.copy(), + var.copy(), + layers=layers_no_splicing.copy(), + ) + if P is not None: + adata_no_splicing.obsm["protein"] = P + # remove cells that has no expression + adata_no_splicing = adata_no_splicing[np.array(adata_no_splicing.X.sum(1)).flatten() > 0, :] return adata, adata_no_splicing diff --git a/dynamo/simulation/ODE.py b/dynamo/simulation/ODE.py index 281e7fc57..b896eb59c 100755 --- a/dynamo/simulation/ODE.py +++ b/dynamo/simulation/ODE.py @@ -80,7 +80,7 @@ def hill_act_grad(x: float, A: float, K: float, n: float, g: float = 0) -> float def toggle( - ab: Union[np.ndarray, Tuple[float, float]], t: Optional[float] = None, beta: float = 5, gamma: float = 1, n: int = 2 + ab: Union[np.ndarray, Tuple[float, float]], beta: float = 5, gamma: float = 1, n: int = 2 ) -> np.ndarray: """Calculates the right-hand side (RHS) of the differential equations for the toggle switch system. @@ -105,7 +105,7 @@ def toggle( return res -def Ying_model(x: np.ndarray, t: Optional[float]=None): +def Ying_model(x: np.ndarray): """network used in the potential landscape paper from Ying, et. al: https://www.nature.com/articles/s41598-017-15889-2. This is also the mixture of Gaussian model. @@ -124,7 +124,7 @@ def Ying_model(x: np.ndarray, t: Optional[float]=None): return ret -def jacobian_Ying_model(x, t=None): +def jacobian_Ying_model(x): """network used in the potential landscape paper from Ying, et. al: https://www.nature.com/articles/s41598-017-15889-2. This is also the mixture of Gaussian model. @@ -333,14 +333,12 @@ def ode_neurongenesis( def neurongenesis( x, - t=None, mature_mu=0, n=4, k=1, a=4, eta=0.25, eta_m=0.125, - eta_b=0.1, a_s=2.2, a_e=6, mx=10, @@ -386,14 +384,10 @@ def neurongenesis( return dx -def hsc(): - pass - - def state_space_sampler(ode, dim, seed_num=19491001, clip=True, min_val=0, max_val=4, N=10000): """Sample N points from the dim dimension gene expression space while restricting the values to be between min_val and max_val. Velocity vector at the sampled points will be calculated according to ode function.""" - seed(seed) + seed(seed_num) X = np.array([[uniform(min_val, max_val) for _ in range(dim)] for _ in range(N)]) Y = np.clip(X + ode(X), a_min=min_val, a_max=None) if clip else X + ode(X) diff --git a/dynamo/tools/moments.py b/dynamo/tools/moments.py index 939cf5fdc..83079fd6c 100755 --- a/dynamo/tools/moments.py +++ b/dynamo/tools/moments.py @@ -366,7 +366,7 @@ def prepare_data_deterministic( ) else: tot_sfs = adata.obs.total_Size_Factor - sfs_x, sfs_y = tot_sfs[:, None], tot_sfs[:, None] + sfs_x, sfs_y = tot_sfs.values[:, None], tot_sfs.values[:, None] m = [None] * len(layers) v = [None] * len(layers) @@ -424,7 +424,7 @@ def prepare_data_deterministic( total_layers=None, CM=pair_x_layer + group_pair_x_layer, ) - sfs_x, sfs_y = sfs_x[:, None], sfs_y[:, None] + sfs_x, sfs_y = sfs_x.values[:, None], sfs_y.values[:, None] x_layer = normalize_mat_monocle( x_layer[:, adata.var_names.isin(genes)], @@ -477,7 +477,7 @@ def prepare_data_deterministic( ) x_layer = normalize_mat_monocle( x_layer[:, adata.var_names.isin(genes)], - szfactors=tot_sfs[:, None], + szfactors=tot_sfs.values[:, None], relative_expr=True, pseudo_expr=0, norm_method=None, @@ -486,7 +486,7 @@ def prepare_data_deterministic( if return_ntr: total_layer = normalize_mat_monocle( total_layer[:, adata.var_names.isin(genes)], - szfactors=tot_sfs[:, None], + szfactors=tot_sfs.values[:, None], relative_expr=True, pseudo_expr=0, norm_method=None,