Add inverse normal normalization (INT) (cytomining#727)

* added InverseNormalTransform (INT)

* added  BaseEstimator, TransformerMixin to the InverseNormalTransform class

* added InverseNormalTransform to the normalize function

* added tests

* added test for cli

* formatting

* added docs

* [pre-commit.ci lite] apply automatic fixes

* updated docs

* applied code rabbit changes

* added tests to reflect changes done from code rabbit's recommendation

* updated docs: clarify InverseNormalTransform n_quantiles behavior

* updated docs

* added test to verify default inverse normal quantiles

* updated __init__.py

* Update pycytominer/operations/transform.py

Co-authored-by: Gregory Way <gregory.way@gmail.com>

---------

Co-authored-by: pre-commit-ci-lite[bot] <117423508+pre-commit-ci-lite[bot]@users.noreply.github.com>
Co-authored-by: Gregory Way <gregory.way@gmail.com>

Author: 3 people

pycytominer/cli.py

@@ -206,6 +206,7 @@ def normalize(
         spherize_center: bool = True,
         spherize_method: str = "ZCA-cor",
         spherize_epsilon: float = 1e-6,
+        inverse_normal_n_quantiles: int = 1000,
     ) -> str:
         """Normalize profiles from a file and write the results to disk.
 
@@ -225,6 +226,9 @@ def normalize(
             spherize_center: Whether to center data before sphering.
             spherize_method: Spherize method to use.
             spherize_epsilon: Spherize epsilon parameter.
+            inverse_normal_n_quantiles: Number of cumulative distribution function
+                landmarks used for inverse normal normalization. Values larger than
+                the number of samples are capped at the number of samples.
 
         Returns:
             The output file path.
@@ -254,6 +258,7 @@ def normalize(
             spherize_center=spherize_center,
             spherize_method=spherize_method,
             spherize_epsilon=spherize_epsilon,
+            inverse_normal_n_quantiles=inverse_normal_n_quantiles,
         )
         if isinstance(result, str):
             _announce_output_file(result)

pycytominer/normalize.py

@@ -11,7 +11,7 @@
 from pycytominer.cyto_utils.features import infer_cp_features
 from pycytominer.cyto_utils.load import load_profiles
 from pycytominer.cyto_utils.util import write_to_file_if_user_specifies_output_details
-from pycytominer.operations import RobustMAD, Spherize
+from pycytominer.operations import InverseNormalTransform, RobustMAD, Spherize
 
 
 @write_to_file_if_user_specifies_output_details
@@ -33,6 +33,7 @@ def normalize(
     spherize_center: bool = True,
     spherize_method: str = "ZCA-cor",
     spherize_epsilon: float = 1e-6,
+    inverse_normal_n_quantiles: int = 1000,
 ) -> pd.DataFrame:
     """Normalize profiling features
 
@@ -103,6 +104,10 @@ def normalize(
     spherize_epsilon : float, default 1e-6.
         The sphering (aka whitening) fudge factor parameter. The function only uses
         this variable if method = "spherize".
+    inverse_normal_n_quantiles : int, default=1000
+        Number of cumulative distribution function landmarks used for the inverse
+        normal transformation. Values larger than the number of samples are capped
+        at the number of samples. Only used when ``method="inverse_normal"``.
 
     Returns
     -------
@@ -187,7 +192,13 @@ def normalize(
     # Define which scaler to use
     method = method.lower()
 
-    avail_methods = ["standardize", "robustize", "mad_robustize", "spherize"]
+    avail_methods = [
+        "standardize",
+        "robustize",
+        "mad_robustize",
+        "spherize",
+        "inverse_normal",
+    ]
     if method not in avail_methods:
         raise ValueError(f"operation must be one {avail_methods}")
 
@@ -206,6 +217,8 @@ def normalize(
             epsilon=spherize_epsilon,
             return_numpy=True,
         )
+    elif method == "inverse_normal":
+        scaler = InverseNormalTransform(n_quantiles=inverse_normal_n_quantiles)
 
     if features == "infer":
         features = infer_cp_features(profiles, image_features=image_features)

pycytominer/operations/__init__.py

@@ -2,5 +2,5 @@
 from .frequency_threshold import calculate_frequency, frequency_threshold
 from .get_na_columns import get_na_columns
 from .noise_removal import noise_removal
-from .transform import RobustMAD, Spherize
+from .transform import InverseNormalTransform, RobustMAD, Spherize
 from .variance_threshold import variance_threshold

pycytominer/operations/transform.py

@@ -11,7 +11,7 @@
 import pandas as pd
 from scipy.stats import median_abs_deviation
 from sklearn.base import BaseEstimator, TransformerMixin
-from sklearn.preprocessing import StandardScaler
+from sklearn.preprocessing import QuantileTransformer, StandardScaler
 
 Spherize_type = TypeVar("Spherize_type", bound="Spherize")
 RobustMAD_type = TypeVar("RobustMAD_type", bound="RobustMAD")
@@ -338,3 +338,92 @@ def transform(self, X: pd.DataFrame, copy: Optional[bool] = None) -> pd.DataFram
             RobustMAD transformed dataframe
         """
         return (X - self.median) / (self.mad + self.epsilon)
+
+
+class InverseNormalTransform(BaseEstimator, TransformerMixin):
+    """Inverse normal transform.
+
+    Apply a rank-based quantile transformation to each feature independently
+    and map the resulting values to a normal distribution.
+
+    1) Rank the values of each feature independently.
+    2) Map the ranks to quantiles of a normal distribution.
+    3) Return the transformed values.
+
+    This class wraps sklearn.preprocessing.QuantileTransformer with
+    output_distribution="normal".
+
+    Parameters
+    ----------
+    n_quantiles : int, default=1000
+        Number of quantiles to be computed. It corresponds to the number of landmarks
+        used to discretize the cumulative distribution function. If ``n_quantiles`` is
+        larger than the number of samples, it is set to the number of samples because
+        a larger number of quantiles does not improve the cumulative distribution
+        function estimate. The actual number used after fitting is available as
+        ``n_quantiles_``. See sklearn.preprocessing.QuantileTransformer for more details.
+    random_state : int, RandomState instance or None, default=None
+        Determines random number generation for smoothing noise. Pass an int for
+        reproducible results across multiple calls.
+
+    Notes
+    -----
+    This transform is rank-based: values are first converted to quantile ranks,
+    then mapped to a normal distribution. The transformed values are therefore
+    normal scores, not the original raw measurements, and distances between raw
+    values are not preserved.
+    """
+
+    def __init__(
+        self,
+        n_quantiles=1000,
+        random_state=None,
+    ):
+        self.n_quantiles = n_quantiles
+        self.random_state = random_state
+
+    def fit(self, x, y=None):
+        """Fit inverse normal transform.
+
+        Parameters
+        ----------
+        x : pandas.DataFrame or numpy.ndarray
+            Data to fit.
+        y : None
+            Has no effect; only used for consistency in sklearn transform API
+
+        Returns
+        -------
+        self
+            Fitted inverse normal transform.
+        """
+        # Set number of quantiles, if n_quantiles is greater than the number of samples\
+        # set it to the number of samples.
+        self.n_quantiles_ = min(self.n_quantiles, x.shape[0])
+
+        # Initialize transformer and set output distribution to normal.
+        # We set it to normal because we want to map the ranks to a normal distribution.
+        self.transformer_ = QuantileTransformer(
+            n_quantiles=self.n_quantiles_,
+            output_distribution="normal",
+            random_state=self.random_state,
+        )
+
+        self.transformer_.fit(x)
+
+        return self
+
+    def transform(self, x):
+        """Apply inverse normal transform.
+
+        Parameters
+        ----------
+        x : pandas.DataFrame or numpy.ndarray
+            Data to transform.
+
+        Returns
+        -------
+        numpy.ndarray
+            Transformed data.
+        """
+        return self.transformer_.transform(x)

tests/test_cli.py

@@ -8,6 +8,7 @@
 import numpy as np
 import pandas as pd
 import pytest
+from sklearn.preprocessing import QuantileTransformer
 
 from pycytominer import cli as pycytominer_cli
 from pycytominer.cli import PycytominerCLI, PycytominerCLIError
@@ -103,6 +104,36 @@ def test_cli_normalize(tmp_path: pathlib.Path) -> None:
     assert np.isclose(result["Feature_2"].mean(), 0.0, atol=1e-7)
 
 
+def test_cli_normalize_inverse_normal(tmp_path: pathlib.Path) -> None:
+    """Ensure CLI normalize forwards inverse normal options."""
+    df, profiles_path = _write_profiles(tmp_path)
+    output_path = tmp_path / "normalized_inverse_normal.csv"
+
+    cli = PycytominerCLI()
+    cli.normalize(
+        profiles=str(profiles_path),
+        output_file=str(output_path),
+        features="Feature_1,Feature_2",
+        meta_features="Metadata_Plate,Metadata_Well",
+        method="inverse_normal",
+        inverse_normal_n_quantiles=3,
+    )
+
+    result = pd.read_csv(output_path)
+    expected_features = QuantileTransformer(
+        n_quantiles=3,
+        output_distribution="normal",
+    ).fit_transform(df.loc[:, ["Feature_1", "Feature_2"]])
+
+    assert result.loc[:, ["Metadata_Plate", "Metadata_Well"]].equals(
+        df.loc[:, ["Metadata_Plate", "Metadata_Well"]]
+    )
+    np.testing.assert_allclose(
+        result.loc[:, ["Feature_1", "Feature_2"]],
+        expected_features,
+    )
+
+
 def test_cli_normalize_drop_cosmicqc_rows(tmp_path: pathlib.Path) -> None:
     """Ensure CLI normalize forwards drop_cosmicqc_rows."""
     profiles = pd.DataFrame({

tests/test_normalize.py

@@ -5,6 +5,7 @@
 import numpy as np
 import pandas as pd
 import pytest
+from sklearn.preprocessing import QuantileTransformer
 
 from pycytominer.normalize import normalize
 
@@ -702,6 +703,83 @@ def test_normalize_robustize_mad_allsamples_novar():
     pd.testing.assert_frame_equal(normalize_result, expected_result)
 
 
+def test_normalize_inverse_normal_inferred_image_profile_features():
+    """
+    Testing normalize pycytominer function
+    method = "inverse_normal"
+    features = "infer"
+    samples = "all"
+    """
+    # Set features to the inferred image profile features
+    features = ["Cells_x", "Cells_y", "Cytoplasm_z", "Nuclei_zz"]
+
+    # Run normalize with the inferred features
+    normalize_result = normalize(
+        profiles=data_feature_infer_df.copy(),
+        features="infer",
+        meta_features=["Metadata_plate", "Metadata_treatment"],
+        samples="all",
+        method="inverse_normal",
+        inverse_normal_n_quantiles=5,
+    )
+
+    # Compute the expected result using QuantileTransformer with the same parameters
+    expected_features = QuantileTransformer(
+        n_quantiles=5,
+        output_distribution="normal",
+    ).fit_transform(data_feature_infer_df.loc[:, features])
+    expected_result = pd.concat(
+        [
+            data_feature_infer_df.loc[:, ["Metadata_plate", "Metadata_treatment"]],
+            pd.DataFrame(expected_features, columns=features),
+        ],
+        axis="columns",
+    )
+
+    pd.testing.assert_frame_equal(normalize_result, expected_result)
+
+
+def test_normalize_inverse_normal_control_samples():
+    """
+    Testing normalize pycytominer function
+    method = "inverse_normal"
+    samples = "Metadata_treatment == 'control'"
+    """
+
+    # Set features to the inferred image profile features
+    # and add a control query to select only control samples
+    features = ["Cells_x", "Cells_y", "Cytoplasm_z", "Nuclei_zz"]
+    control_query = "Metadata_treatment == 'control'"
+
+    # Run normalize with the inferred features and control samples
+    normalize_result = normalize(
+        profiles=data_feature_infer_df.copy(),
+        features=features,
+        meta_features=["Metadata_plate", "Metadata_treatment"],
+        samples=control_query,
+        method="inverse_normal",
+        inverse_normal_n_quantiles=3,
+    )
+
+    # Set up the expected result using QuantileTransformer with the same parameters
+    expected_scaler = QuantileTransformer(
+        n_quantiles=3,
+        output_distribution="normal",
+    ).fit(data_feature_infer_df.query(control_query).loc[:, features])
+    expected_features = expected_scaler.transform(
+        data_feature_infer_df.loc[:, features]
+    )
+    expected_result = pd.concat(
+        [
+            data_feature_infer_df.loc[:, ["Metadata_plate", "Metadata_treatment"]],
+            pd.DataFrame(expected_features, columns=features),
+        ],
+        axis="columns",
+    )
+
+    pd.testing.assert_frame_equal(normalize_result, expected_result)
+
+
 def test_normalize_standardize_allsamples_fromfile():
     """
     Testing normalize pycytominer function