Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Debug simulation module #658

Merged
merged 2 commits into from
Feb 20, 2024
Merged

Debug simulation module #658

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
6402445
update pandas reference in moments
Sichao25 Jan 31, 2024
d2ce291
debug simulation
Sichao25 Feb 1, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
56 changes: 33 additions & 23 deletions dynamo/simulation/Gillespie.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand All @@ -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)]
Expand Down Expand Up @@ -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,
Expand Down Expand Up @@ -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

Expand Down Expand Up @@ -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

Expand Down Expand Up @@ -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)
Expand All @@ -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],
}
Expand All @@ -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
14 changes: 4 additions & 10 deletions dynamo/simulation/ODE.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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.

Expand All @@ -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.
Expand All @@ -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.
Expand Down Expand Up @@ -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,
Expand Down Expand Up @@ -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)

Expand Down
8 changes: 4 additions & 4 deletions dynamo/tools/moments.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -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)],
Expand Down Expand Up @@ -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,
Expand All @@ -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,
Expand Down
Loading