danielee0707
commented
1 year ago Hi, sorry for late reply. Which version of multivelo did you use? It seems the issue is from the mv.recover_dynamics_chrom function. It is a bit difficult to tell the exact location of this error since it's inside the generate_exp function and to find it, you may have to print parameters explicitly. We're working on enhancing the logging ability of the package and will hopefully provide more information while running. A quick workaround is to skip the gene that is causing the issue and run on the rest of the genes.
jacobrepucci
OneHitKO
adata_result = mv.recover_dynamics_chrom(adata_rna, adata_atac, max_iter=5, init_mode="invert", verbose=False, parallel=False, save_plot=False, rna_only=False, fit=True, n_anchors=500, extra_color_key='cell_type2' )
ZeroDivisionError Traceback (most recent call last)
Cell In[18], line 1
----> 1 adata_result = mv.recover_dynamics_chrom(adata_rna,
2 adata_atac,
3 max_iter=5,
4 init_mode="invert",
5 verbose=False,
6 parallel=False,
7 save_plot=False,
8 rna_only=False,
9 fit=True,
10 n_anchors=500,
11 extra_color_key='cell_type2'
12 )
13 sc.write(dir+'/01.RDSsave/2024-06-01_5W_6W_7W_8W_17W_scdata_for_multiVelo_delbatch3000_multiVelo.h5ad', adata_rna)
File ~/.local/lib/python3.9/site-packages/multivelo/dynamical_chrom_func.py:2727, in recover_dynamics_chrom(adata_rna, adata_atac, g ene_list, max_iter, init_mode, model_to_run, verbose, plot, parallel, n_jobs, save_plot, plot_dir, rna_only, fit, fit_decoupling, ex tra_color_key, embedding, n_anchors, k_dist, thresh_multiplier, weight_c, outlier, n_pcs, n_neighbors, fig_size, point_size, partial , direction, rescale_u, alpha, beta, gamma, t_sw)
2724 if verbose >= 1:
2725 print(f'@@@@@fitting {gene}')
2726 (loss, model, direct_out, parameters, initial_exp,
-> 2727 time, state, velocity, likelihood, anchors) = func_to_call(c_mat[:,i],
2728 u_mat[:,i],
2729 s_mat[:,i],
2730 model_to_run[i] if m_per_g else model_to_run,
2731 max_iter,
2732 init_mode,
2733 global_pdist,
2734 embed_coord,
2735 rna_conn,
2736 verbose,
2737 plot,
2738 save_plot,
2739 plot_dir,
2740 fit_args,
2741 gene,
2742 partial[i] if p_per_g else partial,
2743 direction[i] if d_per_g else direction,
2744 rna_only,
2745 fit, 2746 fit_decoupling, [65/1556] 2747 extra_color,
2748 ru[i] if isinstance(ru, (list, np.ndarray)) else ru,
2749 alpha[i] if isinstance(alpha, (list, np.ndarray)) else alpha,
2750 beta[i] if isinstance(beta, (list, np.ndarray)) else beta,
2751 gamma[i] if isinstance(gamma, (list, np.ndarray)) else gamma,
2752 t_sw[i] if isinstance(t_sw, (list, np.ndarray)) else t_sw)
2753 switch, rate, scale_cc, rescale_c, rescale_u, realign_ratio = parameters
2754 likelihood, l_c, ssd_c, var_c = likelihood
File ~/.local/lib/python3.9/site-packages/multivelo/dynamical_chrom_func.py:2153, in regressfunc(c, u, s, m, mi, im, gpdist, embed, conn, v, pl, sp, pdir, fa, gene, pa, di, ro, fit, fd, extra, ru, alpha, beta, gamma, t)
2127 cdc = ChromatinDynamical(c,
2128 u,
2129 s,
(...)
2150 gamma=gamma,
2151 t=t)
2152 if fit:
-> 2153 loss = cdc.fit()
2154 if loss[-1] == np.inf and v >= 1:
2155 print(f'low quality gene {gene}, skipping..')
File ~/.local/lib/python3.9/site-packages/multivelo/dynamical_chrom_func.py:1249, in ChromatinDynamical.fit(self)
1246 self.ax = self.fig.add_subplot(111, projection='3d')
1248 if not self.known_pars:
-> 1249 self.fit_dyn()
1251 self.update(self.params, perform_update=True, fit_outlier=True, plot=True)
1253 # remove long gaps in the last observed state
File ~/.local/lib/python3.9/site-packages/multivelo/dynamical_chrom_func.py:1472, in ChromatinDynamical.fit_dyn(self)
1470 self.update(new_params, adjust_time=False)
1471 if self.model == 0 or self.model == 1:
-> 1472 res = minimize(self.mse, x0=[self.params[0], self.params[1], self.params[3]], method='Nelder-Mead', tol=1e-2, callback=s elf.update, options={'maxiter':20})
1473 elif self.model == 2:
1474 res = minimize(self.mse, x0=[self.params[0], self.params[2], self.params[3]], method='Nelder-Mead', tol=1e-2, callback=s elf.update, options={'maxiter':20})
File ~/anaconda3/envs/scvi-env/lib/python3.9/site-packages/scipy/optimize/_minimize.py:684, in minimize(fun, x0, args, method, jac, hess, hessp, bounds, constraints, tol, callback, options)
681 bounds = standardize_bounds(bounds, x0, meth)
683 if meth == 'nelder-mead':
--> 684 res = _minimize_neldermead(fun, x0, args, callback, bounds=bounds,
685 options)
686 elif meth == 'powell':
687 res = _minimize_powell(fun, x0, args, callback, bounds, options)
File ~/anaconda3/envs/scvi-env/lib/python3.9/site-packages/scipy/optimize/_optimize.py:877, in _minimize_neldermead(func, x0, args, callback, maxiter, maxfev, disp, return_all, initial_simplex, xatol, fatol, adaptive, bounds, **unknown_options)
875 if bounds is not None:
876 xe = np.clip(xe, lower_bound, upper_bound)
--> 877 fxe = func(xe)
879 if fxe < fxr:
880 sim[-1] = xe
File ~/anaconda3/envs/scvi-env/lib/python3.9/site-packages/scipy/optimize/_optimize.py:569, in _wrap_scalar_function_maxfun_validati on..function_wrapper(x, wrapper_args)
567 ncalls[0] += 1
568 # A copy of x is sent to the user function (gh13740)
--> 569 fx = function(np.copy(x), (wrapper_args + args))
570 # Ideally, we'd like to a have a true scalar returned from f(x). For
571 # backwards-compatibility, also allow np.array([1.3]),
572 # np.array([[1.3]]) etc.
573 if not np.isscalar(fx):
File ~/.local/lib/python3.9/site-packages/multivelo/dynamical_chrom_func.py:1651, in ChromatinDynamical.mse(self, x, fit_outlier, pe nalize_gap)
1648 s_array = self.s_all if fit_outlier else self.s
1650 # distances and time assignments
-> 1651 res = calculate_dist_and_time(c_array,
1652 u_array,
1653 s_array,
1654 t_sw_array,
1655 r4[0],
1656 r4[1],
1657 r4[2],
1658 r4[3],
1659 rescale_c,
1660 rescale_u,
1661 scale_cc=scale_cc,
1662 scale_factor=self.scale_factor,
1663 model=self.model,
1664 direction=self.direction,
1665 conn=self.conn if fit_outlier else self.conn_sub,
1666 k=self.k_dist,
1667 t=self.n_anchors,
1668 rna_only=self.rna_only,
1669 penalize_gap=penalize_gap,
1670 all_cells=fit_outlier)
1672 if fit_outlier:
1673 min_dist, t_pred, state_pred, reach_ss, late_phase, max_u, max_s, self.anchor_t1_list, self.anchor_t2_list = res
File ~/.local/lib/python3.9/site-packages/multivelo/dynamical_chrom_func.py:467, in calculate_dist_and_time(c, u, s, t_sw_array, alp ha_c, alpha, beta, gamma, rescale_c, rescale_u, scale_cc, scale_factor, model, conn, t, k, direction, total_h, rna_only, penalize_ga p, all_cells)
465 typed_tau_list = List()
466 [typed_tau_list.append(x) for x in tau_list]
--> 467 exp_list, exp_sw_list = generate_exp(typed_tau_list,
468 t_sw_array[:switch],
469 alpha_c,
470 alpha,
471 beta,
472 gamma,
473 model=model,
474 scale_cc=scale_cc,
475 rna_only=rna_only)
476 rescale_factor = np.array([rescale_c, rescale_u, 1.0])
477 exp_list = [x*rescale_factor for x in exp_list]
ZeroDivisionError: division by zero