mirror of
https://github.com/chanzuckerberg/cellxgene.git
synced 2026-09-30 03:08:12 +08:00
move common code into server, update tests and makefile (#2425)
* move common code into server, update tests and makefile remove backend directory, refactor update smoke tests
This commit is contained in:
Madison Dunitz
@@ -0,0 +1,163 @@
|
||||
import numpy as np
|
||||
from scipy import sparse, stats
|
||||
from server.common.constants import XApproximateDistribution
|
||||
|
||||
|
||||
def diffexp_ttest(adaptor, maskA, maskB, top_n=8, diffexp_lfc_cutoff=0.01):
|
||||
"""
|
||||
Return differential expression statistics for top N variables.
|
||||
|
||||
Algorithm:
|
||||
- compute fold change
|
||||
- compute Welch's t-test statistic and pvalue (w/ Bonferroni correction)
|
||||
- return top N abs(logfoldchange) where lfc > diffexp_lfc_cutoff
|
||||
|
||||
If there are not N which meet criteria, augment by removing the logfoldchange
|
||||
threshold requirement.
|
||||
|
||||
Notes on alogrithm:
|
||||
- Welch's ttest provides basic statistics test.
|
||||
https://en.wikipedia.org/wiki/Welch%27s_t-test
|
||||
- p-values adjusted with Bonferroni correction.
|
||||
https://en.wikipedia.org/wiki/Bonferroni_correction
|
||||
|
||||
:param adaptor: DataAdaptor instance
|
||||
:param maskA: observation selection mask for set 1
|
||||
:param maskB: observation selection mask for set 2
|
||||
:param top_n: number of variables to return stats for
|
||||
:param diffexp_lfc_cutoff: minimum
|
||||
absolute value returning [ varindex, logfoldchange, pval, pval_adj ] for top N genes
|
||||
:return: for top N genes, {"positive": for top N genes, [ varindex, foldchange, pval, pval_adj ],
|
||||
"negative": for top N genes, [ varindex, foldchange, pval, pval_adj ]}
|
||||
"""
|
||||
|
||||
X_approximate_distribution = adaptor.get_X_approximate_distribution()
|
||||
dataA = adaptor.get_X_array(maskA, None)
|
||||
dataB = adaptor.get_X_array(maskB, None)
|
||||
|
||||
# mean, variance, N - calculate for both selections
|
||||
meanA, vA, nA = mean_var_n(dataA, X_approximate_distribution)
|
||||
meanB, vB, nB = mean_var_n(dataB, X_approximate_distribution)
|
||||
res = diffexp_ttest_from_mean_var(meanA, vA, nA, meanB, vB, nB, top_n, diffexp_lfc_cutoff)
|
||||
|
||||
return res
|
||||
|
||||
|
||||
def diffexp_ttest_from_mean_var(meanA, varA, nA, meanB, varB, nB, top_n, diffexp_lfc_cutoff):
|
||||
# IMPORTANT NOTE: this code assumes the data is normally distributed and/or already logged.
|
||||
|
||||
n_var = meanA.shape[0]
|
||||
top_n = min(top_n, n_var)
|
||||
|
||||
# variance / N
|
||||
vnA = varA / min(nA, nB) # overestimate variance, would normally be nA
|
||||
vnB = varB / min(nA, nB) # overestimate variance, would normally be nB
|
||||
sum_vn = vnA + vnB
|
||||
|
||||
# degrees of freedom for Welch's t-test
|
||||
with np.errstate(divide="ignore", invalid="ignore"):
|
||||
dof = sum_vn ** 2 / (vnA ** 2 / (nA - 1) + vnB ** 2 / (nB - 1))
|
||||
dof[np.isnan(dof)] = 1
|
||||
|
||||
# Welch's t-test score calculation
|
||||
with np.errstate(divide="ignore", invalid="ignore"):
|
||||
tscores = (meanA - meanB) / np.sqrt(sum_vn)
|
||||
tscores[np.isnan(tscores)] = 0
|
||||
|
||||
# p-value
|
||||
pvals = stats.t.sf(np.abs(tscores), dof) * 2
|
||||
pvals_adj = pvals * n_var
|
||||
pvals_adj[pvals_adj > 1] = 1 # cap adjusted p-value at 1
|
||||
|
||||
# log fold change. The data is normally distributed/logged, so just subtract the means.
|
||||
logfoldchanges = meanA - meanB
|
||||
|
||||
stats_to_sort = tscores
|
||||
# find all with lfc > cutoff
|
||||
lfc_above_cutoff_idx = np.nonzero(np.abs(logfoldchanges) > diffexp_lfc_cutoff)[0]
|
||||
|
||||
# derive sort order
|
||||
if lfc_above_cutoff_idx.shape[0] > top_n * 2:
|
||||
# partition top N
|
||||
rel_t_partition = np.argpartition(stats_to_sort[lfc_above_cutoff_idx], (top_n, -top_n))
|
||||
rel_t_partition_top_n = np.concatenate((rel_t_partition[-top_n:], rel_t_partition[:top_n]))
|
||||
t_partition = lfc_above_cutoff_idx[rel_t_partition_top_n]
|
||||
# sort the top N partition
|
||||
rel_sort_order = np.argsort(stats_to_sort[t_partition])[::-1]
|
||||
sort_order = t_partition[rel_sort_order]
|
||||
else:
|
||||
# partition and sort top N, ignoring lfc cutoff
|
||||
partition = np.argpartition(stats_to_sort, (top_n, -top_n))
|
||||
partition_top_n = np.concatenate((partition[-top_n:], partition[:top_n]))
|
||||
|
||||
rel_sort_order = np.argsort(stats_to_sort[partition_top_n])[::-1]
|
||||
indices = np.indices(stats_to_sort.shape)[0]
|
||||
sort_order = indices[partition_top_n][rel_sort_order]
|
||||
|
||||
# top n slice based upon sort order
|
||||
logfoldchanges_top_n = logfoldchanges[sort_order]
|
||||
pvals_top_n = pvals[sort_order]
|
||||
pvals_adj_top_n = pvals_adj[sort_order]
|
||||
|
||||
# varIndex, logfoldchange, pval, pval_adj
|
||||
result = {
|
||||
"positive": [
|
||||
[sort_order[i], logfoldchanges_top_n[i], pvals_top_n[i], pvals_adj_top_n[i]] for i in range(top_n)
|
||||
],
|
||||
"negative": [
|
||||
[sort_order[i], logfoldchanges_top_n[i], pvals_top_n[i], pvals_adj_top_n[i]]
|
||||
for i in range(-1, -1 - top_n, -1)
|
||||
],
|
||||
}
|
||||
|
||||
return result
|
||||
|
||||
|
||||
# Convenience function which handles sparse data
|
||||
def mean_var_n(X, X_approximate_distribution=XApproximateDistribution.NORMAL):
|
||||
"""
|
||||
Two-pass variance calculation. Numerically (more) stable
|
||||
than naive methods (and same method used by numpy.var())
|
||||
https://en.wikipedia.org/wiki/Algorithms_for_calculating_variance#Two-pass
|
||||
"""
|
||||
# fp_err_occurred is a flag indicating that a floating point error
|
||||
# occured somewhere in our compute. Used to trigger non-finite
|
||||
# number handling.
|
||||
fp_err_occurred = False
|
||||
|
||||
def fp_err_set(err, flag):
|
||||
nonlocal fp_err_occurred
|
||||
fp_err_occurred = True
|
||||
|
||||
with np.errstate(divide="call", invalid="call", call=fp_err_set):
|
||||
n = X.shape[0]
|
||||
if sparse.issparse(X):
|
||||
if X_approximate_distribution == XApproximateDistribution.COUNT:
|
||||
X = X.log1p()
|
||||
mean = X.mean(axis=0).A1
|
||||
dfm = X - mean
|
||||
sumsq = np.sum(np.multiply(dfm, dfm), axis=0).A1
|
||||
v = sumsq / (n - 1)
|
||||
else:
|
||||
if X_approximate_distribution == XApproximateDistribution.COUNT:
|
||||
X = np.log1p(X)
|
||||
mean = X.mean(axis=0)
|
||||
dfm = X - mean
|
||||
sumsq = np.sum(np.multiply(dfm, dfm), axis=0)
|
||||
v = sumsq / (n - 1)
|
||||
|
||||
# AnnData does not guarantee that operations on a view of X will
|
||||
# return an ndarray, so force the cast if it wasn't done for us.
|
||||
if type(mean) is not np.ndarray:
|
||||
mean = mean.toarray()
|
||||
if type(v) is not np.ndarray:
|
||||
v = v.toarray()
|
||||
|
||||
if fp_err_occurred:
|
||||
mean[np.isfinite(mean) == False] = 0 # noqa: E712
|
||||
v[np.isfinite(v) == False] = 0 # noqa: E712
|
||||
else:
|
||||
mean[np.isnan(mean)] = 0
|
||||
v[np.isnan(v)] = 0
|
||||
|
||||
return mean, v, n
|
||||
@@ -0,0 +1,108 @@
|
||||
from typing import Tuple
|
||||
import numba
|
||||
import concurrent.futures
|
||||
import numpy as np
|
||||
from scipy import sparse
|
||||
from server.common.constants import XApproximateDistribution
|
||||
|
||||
|
||||
@numba.njit(error_model="numpy", nogil=True)
|
||||
def min_max_fast(arr: np.ndarray) -> Tuple[float, float]:
|
||||
"""Return (min, max) values for the ndarray."""
|
||||
|
||||
# initialize to first finite value in array. Normally,
|
||||
# this will exit on the first value.
|
||||
for i in range(arr.size):
|
||||
min_val = max_val = arr[i]
|
||||
if np.isfinite(min_val):
|
||||
break
|
||||
|
||||
# now find min/max, unrolled by two
|
||||
odd = arr.size % 2
|
||||
unrolled_loop_limit = arr.size - 1 if odd else arr.size
|
||||
i = 0
|
||||
while i < unrolled_loop_limit:
|
||||
x = arr[i]
|
||||
y = arr[i + 1]
|
||||
|
||||
# ignore non-finites
|
||||
x = x if np.isfinite(x) else min_val
|
||||
y = y if np.isfinite(y) else min_val
|
||||
|
||||
if x > y:
|
||||
x, y = y, x
|
||||
min_val = min(x, min_val)
|
||||
max_val = max(y, max_val)
|
||||
i += 2
|
||||
|
||||
# handle the tail if any
|
||||
if odd:
|
||||
x = arr[arr.size - 1]
|
||||
|
||||
# ignore non-finites
|
||||
x = x if np.isfinite(x) else min_val
|
||||
|
||||
min_val = min(x, min_val)
|
||||
max_val = max(x, max_val)
|
||||
|
||||
return min_val, max_val
|
||||
|
||||
|
||||
def min_max_numpy(arr: np.ndarray) -> Tuple[float, float]:
|
||||
return arr.min(), arr.max()
|
||||
|
||||
|
||||
def numba_has_support_for_scalar_type(arr: np.ndarray) -> bool:
|
||||
"""Numba does not support half-floats, 128 bit floats, ints > 64 bit or non-scalars."""
|
||||
if arr.dtype == np.float32 or arr.dtype == np.float64:
|
||||
return True
|
||||
|
||||
if np.issubdtype(arr.dtype, np.integer) and arr.dtype <= np.int64:
|
||||
return True
|
||||
|
||||
if arr.dtype == np.bool_:
|
||||
return True
|
||||
|
||||
return False
|
||||
|
||||
|
||||
def estimate_approximate_distribution(X) -> XApproximateDistribution:
|
||||
"""
|
||||
Estimate the distribution (normal, count) of the X matrix.
|
||||
|
||||
Currently this is based upon the assumption that scRNA-seq data is
|
||||
exponentially distributed in its raw (count) form, and when logged,
|
||||
any (max-min) range in excess of 24 is implies tens of millions of
|
||||
observations of a single feature and so is extremely unlikely.
|
||||
"""
|
||||
if X.dtype.kind not in ["i", "u", "f"]:
|
||||
raise TypeError(f"Unsupported matrix dtype: {X.dtype.name}")
|
||||
|
||||
if X.size == 0:
|
||||
# default for empty array
|
||||
return XApproximateDistribution.NORMAL
|
||||
|
||||
if sparse.isspmatrix_csc(X) or sparse.isspmatrix_csr(X):
|
||||
Xdata = X.data
|
||||
elif type(X) is np.ndarray:
|
||||
Xdata = X.reshape(
|
||||
X.size,
|
||||
)
|
||||
else:
|
||||
raise TypeError(f"Unsupported matrix format: {str(type(X))}")
|
||||
|
||||
min_max = min_max_fast if numba_has_support_for_scalar_type(Xdata) else min_max_numpy
|
||||
|
||||
CHUNKSIZE = 1 << 24
|
||||
if Xdata.size > CHUNKSIZE:
|
||||
min_val = max_val = Xdata[0]
|
||||
with concurrent.futures.ThreadPoolExecutor() as tp:
|
||||
for (_min, _max) in tp.map(min_max, [Xdata[i : i + CHUNKSIZE] for i in range(0, Xdata.size, CHUNKSIZE)]):
|
||||
min_val = min(_min, min_val)
|
||||
max_val = max(_max, max_val)
|
||||
|
||||
else:
|
||||
min_val, max_val = min_max(Xdata)
|
||||
|
||||
excess_range = (max_val - min_val) > 24
|
||||
return XApproximateDistribution.COUNT if excess_range else XApproximateDistribution.NORMAL
|
||||
Reference in New Issue
Block a user