As far as I can tell what you’re asking for is impossible. Negative binomial is an overdispersed Poisson and sum of Poisson variables is just a Poisson with rate equal to the sum of the individual rates. Only the sum is identifiable, no information about the individuals remains. Now, overdispersion complicates the story so maybe there is something that can be learned beyond the collective parameters but it’s not much.
You’re mushing together multiple cells and then measure the average gene expression level. Why do you think you’ll be able to disentangle them with a statistical model?
EDIT:
If we add nuisance parameters to represent the negative binomial as gamma-poisson mixture it’s possible to marginalize out the discrete parameters.
data {
int N; // # of trials
int y[N]; // # of failures in trial n
int <lower=1> C; // # random variables (cells) in each N
}
parameters {
positive_ordered[C] beta_c;
real<lower=0> alpha_c[C];
real<lower=0> lambda[N,C]; // additional latent parameters
}
model {
for (n in 1:N) {
lambda[n] ~ gamma(alpha_c, beta_c);
y[n] ~ poisson(sum(lambda[n]));
}
}