welch-lab / MultiVelo

Multi-omic velocity inference
BSD 3-Clause "New" or "Revised" License
105 stars 12 forks source link

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

Closed Aitar closed 2 years ago

Aitar commented 2 years ago

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.

danielee0707 commented 2 years ago

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 2 years ago

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

danielee0707 commented 2 years ago

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.

danielee0707 commented 2 years ago

The issue is fixed in the new release.

Aitar commented 2 years ago

Thanks!