welch-lab/MultiVelo

Questions about error "cannot compute length" at line 1478.

Aitar opened this issue · 5 comments

Aitar commented

Hi,MultiVelo team,Thank you for giving us such a great tool!
A error ocurred when I ran MultiVelo on my own data due to dividing by zero. It happened in the code block below at line 726 in file "dynamical_chrom_func.py". At line "s_inf = alpha / gamma", the result is "inf" because gamma is zero. Because at line 1399 self.gamma = np.dot(ss_u, ss_s) / np.dot(ss_s, ss_s) the result of dot product of vectors "ss_u" and "ss_s" is zero.

def approx_tau(u, s, u0, s0, alpha, beta, gamma):
    if gamma == beta:
        gamma -= 1e-3
    u_inf = alpha / beta
    if beta > gamma:
        b_new = beta / (gamma - beta)
        s_inf = alpha / gamma
        s_inf_new = s_inf - b_new * u_inf
        s_new = s - b_new * u
        s0_new = s0 - b_new * u0
        tau = -1.0 / gamma * log_valid((s_new - s_inf_new) / (s0_new - s_inf_new))
    else:
        tau = -1.0 / beta * log_valid((u - u_inf) / (u0 - u_inf))
    return tau

Therefore, at line 1475 rna_interval = approx_tau(u0_r, s0_r, 0, 0, alpha, self.beta, self.gamma), the result "rna_interval" is "nan", which causes the error at line 1478 for t_sw_1 in np.arange(1, rna_interval-1, 2, dtype=np.float64): "cannot compute length".
Thanks for your answering.

Hi, thanks for the error report. It's interesting to see that either the steady-state unspliced or the steady-state spliced expression value is 0. Do you know which gene(s) was causing this issue, and can you plot this gene's RNA phase portrait with scvelo.pl.scatter? I'll add a check to the method, but for now, you can exclude this gene from the "gene_list'' argument to skip it.

Aitar commented

Thanks for your reply. The reason why result of dot product is zero is not the steady-state unspliced or the steady-state spliced expression value is zero but ss_s and ss_u are orthogonal. The RNA phase portrait is below.
image

Thanks, Indeed, there isn't a well-defined steady-state region, but still, somehow this gene has enough cells in the middle that make it pass the quality filter.

The issue is fixed in the new release.

Aitar commented

Thanks!