Repository: theislab/cellrank (language: python)

Issue:
`nan` for q values in lineage driver test
## Description

When running

```
g = GPCCA(kernel)
# ...
g.compute_lineage_drivers(return_drivers=False)
```

all q values are `nan`. Is there a reasonable explanation for this, @michalk8?


Discussion / hints:
Just checked: This is also the case in the real time kernel tutorial.
I've also observed this before. 
> Just checked: This is also the case in the real time kernel tutorial.

This can be because some gene/lineage has a standard deviation of zero, which results in division by zero, causing the correlation to be `NaN`, which are then later using for `pval/qval` calculation. Are all p-values NaN as well?
> This can be because some gene/lineage has a standard deviation of zero, which results in division by zero, causing the correlation to be NaN, which are then later using for pval/qval calculation. Are all p-values NaN as well?

No, the p values are not `nan`.
Strange, this should only happen when the p-values are NaN; running `adata.varm['terminal_lineage_drivers']['Neural_1_pval'].isnull().sum()` in the WOT tutorial,
I get `36` NaN values. Looking at those genes, they all have 0 expression.
Since the multiple tests correction can't handle NaN p-values, maybe we just subset non-NaN p-values for the q-value computation.
Ok, and what about the case in which the p-values are not NaN? In any case, we should issue a warning for zero-expressed genes (they should have been filtered out before by the user with any standard scanpy filtering anyways). 
> Ok, and what about the case in which the p-values are not NaN? In any case, we should issue a warning for zero-expressed genes (they should have been filtered out before by the user with any standard scanpy filtering anyways).

As I've suggested, for q-value calculation, let's subset only non-NaN p-values, which occurr when the gene expression for the gene is constant. We can also add the warning that some genes have no variability.
@michalk8, the q-values are also NaN when the p-values are not NaN.
> @michalk8, the q-values are also NaN when the p-values are not NaN.

If any p-value is NaN, then all q-values are NaN (that's what `statsmodels` does).
I agree, let's go for your suggestion then @michalk8 

Fix patch:
diff --git a/cellrank/tl/_utils.py b/cellrank/tl/_utils.py
index 67e346022..114065911 100644
--- a/cellrank/tl/_utils.py
+++ b/cellrank/tl/_utils.py
@@ -417,14 +417,29 @@ def _correlation_test(
         confidence_level=confidence_level,
         **kwargs,
     )
-    invalid = np.sum((corr < -1) | (corr > 1))
-    if invalid:
-        raise ValueError(f"Found `{invalid}` correlations that are not in `[0, 1]`.")
+    invalid = (corr < -1) | (corr > 1)
+    if np.any(invalid):
+        logg.warning(
+            f"Found `{np.sum(invalid)}` correlation(s) that are not in `[0, 1]`. "
+            f"This usually happens when gene expression is constant across all cells. "
+            f"Setting to `NaN`"
+        )
+        corr[invalid] = np.nan
+        pvals[invalid] = np.nan
+        ci_low[invalid] = np.nan
+        ci_high[invalid] = np.nan
 
     res = pd.DataFrame(corr, index=gene_names, columns=[f"{c}_corr" for c in Y.names])
     for idx, c in enumerate(Y.names):
-        res[f"{c}_pval"] = pvals[:, idx]
-        res[f"{c}_qval"] = multipletests(pvals[:, idx], alpha=0.05, method="fdr_bh")[1]
+        p = pvals[:, idx]
+        valid_mask = ~np.isnan(p)
+
+        res[f"{c}_pval"] = p
+        res[f"{c}_qval"] = np.nan
+        if np.any(valid_mask):
+            res.loc[gene_names[valid_mask], f"{c}_qval"] = multipletests(
+                p[valid_mask], alpha=0.05, method="fdr_bh"
+            )[1]
         res[f"{c}_ci_low"] = ci_low[:, idx]
         res[f"{c}_ci_high"] = ci_high[:, idx]
 


Test patch:
diff --git a/tests/test_gpcca.py b/tests/test_gpcca.py
index 0fbe02644..44b5751c4 100644
--- a/tests/test_gpcca.py
+++ b/tests/test_gpcca.py
@@ -1087,6 +1087,21 @@ def test_params(self, adata_large: AnnData, state: State):
         _assert_params(g, state, fwd=False)
         _assert_params(g, state, fwd=True)
 
+    def test_drivers_constant_gene_qvalues(self, adata_large: AnnData):
+        ix = 0
+        gene = adata_large.var_names[ix]
+        adata_large.X[:, ix] = 0
+        ck = ConnectivityKernel(adata_large).compute_transition_matrix()
+        g = cr.tl.estimators.GPCCA(ck)
+        g.fit(n_states=2)
+        g.predict()
+        g.compute_absorption_probabilities()
+
+        drivers = g.compute_lineage_drivers(use_raw=False)
+
+        np.testing.assert_array_equal(drivers.loc[gene], np.nan)
+        assert np.asarray(drivers.iloc[drivers.index != gene].isnull()).sum() == 0
+
 
 class TestGPCCASerialization:
     @pytest.mark.parametrize("state", list(State))