ValueError when using simulated data

[hannanebelung]

Hi,

I was trying to run pyrovelocity with simulated data . In particular, I was trying to test it with different kinetics parameters as was was done in another paper. The following function was used to generate the data with the scVelo simulation function

def get_adata(n_obs, mu, R, C, n_vars=1000):
    cov = C.dot(C.T) * R

    alpha, beta, gamma = np.exp(np.random.multivariate_normal(mu, cov, size=n_vars).T)  # multivariate log-normal 
    beta /= 3
    gamma /= 3

    # remove outliers
    idx = (alpha < np.percentile(alpha, 99)) & (beta < np.percentile(beta, 99)) & (gamma < np.percentile(gamma, 99))
    alpha = alpha[idx]
    beta = beta[idx]
    gamma = gamma[idx]
    n_vars = np.sum(idx)

    alpha = alpha +1

    switches = np.random.uniform(.1, .5, size=n_vars)

    adata = scv.datasets.simulation(
        n_obs=n_obs,
        t_max=20,
        n_vars=n_vars,
        switches=switches,
        alpha=alpha,
        beta=beta,
        gamma=gamma,
        noise_level=.8,
        random_seed = 0
    )
    
    return adata

with the following parameters and model

mu = np.array([1, .2, .05])

R = np.array([[1., .2, .2],
              [.2, 1., .8],
              [.2, .8, 1.]])

C = np.array([.4, .4, .4])[:, None]
adata = get_adata(n_obs=1000, mu=mu, R=R, C=C)

adata.layers['raw_spliced'] = adata.layers['spliced'].copy()
adata.layers['raw_unspliced'] = adata.layers['unspliced'].copy()
adata.obs['u_lib_size_raw'] = adata.layers['raw_unspliced'].sum(axis=1)
adata.obs['s_lib_size_raw'] = adata.layers['raw_spliced'].sum(axis=1)

scv.pp.pca(adata)
scv.tl.velocity(adata)
scv.pp.neighbors(adata)
scv.tl.velocity(adata, mode='steady_state')

num_epochs = 1000 # large data
# num_epochs = 4000 # small data
adata_model_pos = train_model(adata,
                            max_epochs=num_epochs, svi_train=True, log_every=100,
                            patient_init=45,
                            batch_size=500, 
                            use_gpu=use_gpu,
                            include_prior=True,
                            offset=False,
                            library_size=True,
                            patient_improve=1e-3,
                            model_type='auto',
                            guide_type='auto_t0_constraint',
                            train_size=1.0)

When I try to train the model with the data, the code runs into a value error when calculating a log probability at site u.

Traceback (most recent call last):
  File "/rna-velocity/src/simulations.py", line 309, in <module>
    correlation_df = correlation_df.append(
  File "/rna-velocity/src/simulations.py", line 176, in fit_pyrovelocity
    batch_size=500, 
  File ".conda/envs/pyrovelocity_bioconda/lib/python3.8/site-packages/pyrovelocity/api.py", line 64, in train_model
    losses = model.train_faster_with_batch(max_epochs=max_epochs,
  File "/.conda/envs/pyrovelocity_bioconda/lib/python3.8/site-packages/pyrovelocity/_trainer.py", line 265, in train_faster_with_batch
    loss = svi.step(*args, **kwargs)
  File "/.conda/envs/pyrovelocity_bioconda/lib/python3.8/site-packages/pyro/infer/svi.py", line 128, in step
    loss = self.loss_and_grads(self.model, self.guide, *args, **kwargs)
  File "/.conda/envs/pyrovelocity_bioconda/lib/python3.8/site-packages/pyro/infer/traceenum_elbo.py", line 404, in loss_and_grads
    for model_trace, guide_trace in self._get_traces(model, guide, args, kwargs):
  File "/.conda/envs/pyrovelocity_bioconda/lib/python3.8/site-packages/pyro/infer/traceenum_elbo.py", line 327, in _get_traces
    self._guess_max_plate_nesting(model, guide, args, kwargs)
  File "/.conda/envs/pyrovelocity_bioconda/lib/python3.8/site-packages/pyro/infer/elbo.py", line 109, in _guess_max_plate_nesting
    model_trace.compute_log_prob()
  File "/.conda/envs/pyrovelocity_bioconda/lib/python3.8/site-packages/pyro/poutine/trace_struct.py", line 220, in compute_log_prob
    raise ValueError("Error while computing log_prob at site '{}':\n{}\n{}"
  File "/.conda/envs/pyrovelocity_bioconda/lib/python3.8/site-packages/pyro/poutine/trace_struct.py", line 216, in compute_log_prob
    log_p = site["fn"].log_prob(site["value"], *site["args"], **site["kwargs"])
  File "/.conda/envs/pyrovelocity_bioconda/lib/python3.8/site-packages/torch/distributions/poisson.py", line 61, in log_prob
    self._validate_sample(value)
  File "/.conda/envs/pyrovelocity_bioconda/lib/python3.8/site-packages/torch/distributions/distribution.py", line 277, in _validate_sample
    raise ValueError('The value argument must be within the support')
ValueError: Error while computing log_prob at site 'u':
The value argument must be within the support
       Trace Shapes:          
        Param Sites:          
       Sample Sites:          
          alpha dist     972 |
               value     972 |
            log_prob     972 |
          gamma dist     972 |
               value     972 |
            log_prob     972 |
           beta dist     972 |
               value     972 |
            log_prob     972 |
             t0 dist     972 |
               value     972 |
            log_prob     972 |
   dt_switching dist     972 |
               value     972 |
            log_prob     972 |
          u_inf dist     972 |
               value     972 |
            log_prob     972 |
          s_inf dist     972 |
               value     972 |
            log_prob     972 |
      switching dist     972 |
               value     972 |
            log_prob     972 |
      cell_time dist 500   1 |
               value 500   1 |
            log_prob 500   1 |
   u_read_depth dist 500   1 |
               value 500   1 |
            log_prob 500   1 |
   s_read_depth dist 500   1 |
               value 500   1 |
            log_prob 500   1 |
time_constraint dist 500 972 |
               value         |
            log_prob 500 972 |
cell_gene_state dist 500 972 |
               value 500 972 |
            log_prob 500 972 |
             ut dist 500 972 |
               value 500 972 |
            log_prob 500 972 |
             st dist 500 972 |
               value 500 972 |
            log_prob 500 972 |
        ut_norm dist 500 972 |
               value 500 972 |
            log_prob 500 972 |
        st_norm dist 500 972 |
               value 500 972 |
            log_prob 500 972 |
              u dist 500 972 |
               value 500 972 |

I am not quite sure why this error occurs. The training does seem to work if I use a simpler simulation like this one:

scv.datasets.simulation(
            random_seed=0,
            n_obs=1000,
            n_vars=12,
            alpha=5,
            beta=0.5,
            gamma=0.3,
            noise_model="gillespie"
 )

However if I change the noise_model to 'normal', the same ValueError when calculating the log_prob as before occurs. Have you by any chance tried to run pyrovelocity with more complex simulations? Or do you know why this type of data is not supported by the distribution?

I am not sure if this is a bug or if I am just using the wrong data for the model, so I would appreciate any help.

[cameronraysmith]

Hello @hannanebelung ,

In the first simulation you are likely generating data which are not integer counts. When you set noise_model="gillespie" you are likely generating integer count data. Since the latter may be the only mode of that simulator that respects the domain of real data, we would only endeavor to support this mode. In this case, the simulation still does not capture even the so-called random telegraph approximation to single gene expression that has been validated in single-molecule fluorescence experiments. Of course, this also does not capture the impacts of regulation. This is all to say that, to date, there is not an ideal simulator representative of the phenomenon of interest if that is what you are looking for.

There may be other issues as well. Note that pyrovelocity is under active development and the model is likely to change substantially in relatively short periods of time. Thank you for bringing this point to our attention as it is something we should document clearly. Feel free to let us know if you experience further issues.

[hannanebelung]

Hi @cameronraysmith,
thanks for the answer, that is good to know. I will definitely be on a lookout for any changes or updates of the model.