Fix bug with the layer parameter (#168)

* fix layer parameter

* fix layer parameter

* [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci

* make layer an optional parameter to prevent failing tests

* Fix tests

---------

Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Co-authored-by: Gregor Sturm <mail@gregor-sturm.de>

Author: 3 people

src/infercnvpy/tl/_infercnv.py

@@ -108,7 +108,7 @@ def infercnv(
         var_mask = var_mask | adata.var["chromosome"].isin(exclude_chromosomes)
 
     tmp_adata = adata[:, ~var_mask]
-    reference = _get_reference(adata, reference_key, reference_cat, reference)[:, ~var_mask]
+    reference = _get_reference(adata, reference_key, reference_cat, reference, layer)[:, ~var_mask]
 
     expr = tmp_adata.X if layer is None else tmp_adata.layers[layer]
 
@@ -358,6 +358,7 @@ def _get_reference(
     reference_key: str | None,
     reference_cat: None | str | Sequence[str],
     reference: np.ndarray | None,
+    layer: str | None,
 ) -> np.ndarray:
     """Parameter validation extraction of reference gene expression.
 
@@ -367,13 +368,18 @@ def _get_reference(
     Returns a 2D array with reference categories in rows, cells in columns.
     If there's just one category, it's still a 2D array.
     """
+    if layer is not None:
+        X = adata.layers[layer]
+    else:
+        X = adata.X
+
     if reference is None:
         if reference_key is None or reference_cat is None:
             logging.warning(
                 "Using mean of all cells as reference. For better results, "
                 "provide either `reference`, or both `reference_key` and `reference_cat`. "
             )  # type: ignore
-            reference = np.mean(adata.X, axis=0)
+            reference = np.mean(X, axis=0)
 
         else:
             obs_col = adata.obs[reference_key]
@@ -388,7 +394,7 @@ def _get_reference(
                     f"{reference_cat[~reference_cat_in_obs]}"
                 )
 
-            reference = np.vstack([np.mean(adata.X[obs_col.values == cat, :], axis=0) for cat in reference_cat])
+            reference = np.vstack([np.mean(X[obs_col.values == cat, :], axis=0) for cat in reference_cat])
 
     if reference.ndim == 1:
         reference = reference[np.newaxis, :]

tests/test_tools.py

@@ -10,7 +10,7 @@
 
 def test_get_reference_key_and_cat(adata_mock):
     """Test that reference is correctly calculated given a reference key and category"""
-    actual = _get_reference(adata_mock, "cat", ["foo", "baz"], None)
+    actual = _get_reference(adata_mock, "cat", ["foo", "baz"], None, layer=None)
     npt.assert_almost_equal(
         actual,
         np.array(
@@ -24,19 +24,19 @@ def test_get_reference_key_and_cat(adata_mock):
 
 def test_get_reference_no_reference(adata_mock):
     """If no reference is specified, the mean of the entire adata object is taken"""
-    actual = _get_reference(adata_mock, None, None, None)
+    actual = _get_reference(adata_mock, None, None, None, layer=None)
     npt.assert_almost_equal(actual, np.array([[4.8, 4.2, 4.4, 5]]), decimal=5)
 
 
 def test_get_reference_given_reference(adata_mock):
     """Predefined reference takes precendence over reference_key and reference_cat"""
     reference = np.array([1, 2, 3, 4])
-    actual = _get_reference(adata_mock, "foo", "bar", reference)
+    actual = _get_reference(adata_mock, "foo", "bar", reference, layer=None)
     npt.assert_equal(reference, actual[0, :])
 
     with pytest.raises(ValueError):
         reference = np.array([1, 2, 3])
-        actual = _get_reference(adata_mock, "foo", "bar", reference)
+        actual = _get_reference(adata_mock, "foo", "bar", reference, layer=None)
 
 
 @pytest.mark.parametrize(
@@ -141,7 +141,7 @@ def test_calculate_gene_averages():
 
 
 def test_infercnv_chunk_with_gene_values(adata_full_mock, gene_res_actual, x_res_actual):
-    reference = _get_reference(adata_full_mock, reference_key=None, reference_cat=None, reference=None)
+    reference = _get_reference(adata_full_mock, reference_key=None, reference_cat=None, reference=None, layer=None)
     var = adata_full_mock.var.loc[:, ["chromosome", "start", "end"]]
     tmp_x = adata_full_mock.X
 
@@ -156,7 +156,7 @@ def test_infercnv_chunk_with_gene_values(adata_full_mock, gene_res_actual, x_res
 
 
 def test_infercnv_chunk_default(adata_full_mock, x_res_actual):
-    reference = _get_reference(adata_full_mock, reference_key=None, reference_cat=None, reference=None)
+    reference = _get_reference(adata_full_mock, reference_key=None, reference_cat=None, reference=None, layer=None)
     var = adata_full_mock.var.loc[:, ["chromosome", "start", "end"]]
     tmp_x = adata_full_mock.X
 
@@ -218,6 +218,27 @@ def test_workflow(adata_oligodendroma):
     cnv.pl.chromosome_heatmap_summary(adata_oligodendroma, show=False)
 
 
+def test_layer_parameter(adata_oligodendroma):
+    adata_oligodendroglioma = cnv.datasets.oligodendroglioma()
+
+    # create lognorm layer but leave X unchanged
+    adata_oligodendroglioma.layers["LogNormalize"] = adata_oligodendroglioma.X.copy()
+    sc.pp.normalize_total(adata_oligodendroglioma, layer="LogNormalize", target_sum=1e4)
+    sc.pp.log1p(adata_oligodendroglioma, layer="LogNormalize")
+
+    # copy to test if results change with layer option
+    adata_oligodendroglioma2 = adata_oligodendroglioma.copy()
+    adata_oligodendroglioma2.X = adata_oligodendroglioma.layers["LogNormalize"]
+
+    cnv.tl.infercnv(adata_oligodendroglioma, layer="LogNormalize")
+    cnv.tl.infercnv(adata_oligodendroglioma2, layer=None)
+
+    X_cnv = adata_oligodendroglioma.obsm["X_cnv"].toarray()
+    X_cnv2 = adata_oligodendroglioma2.obsm["X_cnv"].toarray()
+
+    assert np.all(X_cnv == X_cnv2), "Different results found with infercnv layer parameter"
+
+
 def test_calculate_gene_values_speed(benchmark, adata_oligodendroma):
     # Benchmark with calculate_gene_values=True
     benchmark(