def nb_exact_test(x_a, x_b, size_factor_a, size_factor_b, mu, phi):
""" Compute p-value for a pairwise exact test using the negative binomial.
Args:
x_a (int) - Total count for a single gene in group A
x_b (int) - Total count for a single gene in group B
size_factor_a (float) - Sum of size factors for group A
size_factor_b (float) - Sum of size factors for group B
mu (float) - Common mean count for this gene
phi (float) - Common dispersion for this gene
Returns:
p-value (float); the probability that a random pair of counts under the null hypothesis is more extreme
than the observed pair of counts. """
size_factor_a = float(size_factor_a)
size_factor_b = float(size_factor_b)
mu = float(mu)
phi = float(phi)
if (x_a + x_b) == 0:
return 1.0
if phi == 0:
return 1.0
if size_factor_a == 0 or size_factor_b == 0:
return 1.0
all_x_a = np.arange(0, 1+x_a+x_b, 1)
all_x_b = np.arange(x_a+x_b, -1, -1)
def log_prob(x, size_factor):
return neg_bin_log_pmf(x, size_factor*mu, phi/size_factor)
log_p_obs = log_prob(x_a, size_factor_a) + log_prob(x_b, size_factor_b)
log_p_all = log_prob(all_x_a, size_factor_a) + log_prob(all_x_b, size_factor_b)
more_extreme = log_p_all <= log_p_obs
if np.sum(more_extreme) == 0:
return 0.0
return np.exp(logsumexp(log_p_all[log_p_all <= log_p_obs]) - logsumexp(log_p_all))
评论列表
文章目录
