Development

selectlfq/featureengineering.py

@@ -1,9 +1,9 @@
 import numpy as np
+import pandas as pd
 from joblib import Parallel, delayed
 import scipy.stats as stats
 from sklearn.impute import SimpleImputer
 from numba import njit
-
 from selectlfq.utils import get_logger
 
 logger = get_logger()
@@ -51,27 +51,40 @@ def _calculate_variance_distance(data):
 
 
 @njit
-def _nan_correlation_matrix(data):
+def _nan_correlation_matrix(data: np.ndarray) -> np.ndarray:
+    """
+    Calculate the correlation matrix of the data.
+
+    Parameters:
+        data: numpy array of shape (n_fragments, n_samples)
+            The data to calculate the correlation matrix of.
+
+    Returns:
+        numpy array of shape (n_fragments, n_fragments)
+        containing the correlation matrix of the data.
+    """
     n = len(data)
-    correlation_matrix = np.full((n, n), np.nan)  # Initialize with NaN
+    correlation_matrix = np.zeros((n, n))
 
     for i in range(n):
-        for j in range(n):
-            if i != j:  # Only compute for off-diagonal elements
-                mask = np.isfinite(data[i]) & np.isfinite(data[j])
-                if np.sum(mask) > 1:  # Ensure there are at least two data points
-                    xi = data[i][mask]
-                    xj = data[j][mask]
-                    std_dev_i = np.std(xi)
-                    std_dev_j = np.std(xj)
-
-                    if (std_dev_i > 0) and (std_dev_j > 0):
-                        mean_i = np.mean(xi)
-                        mean_j = np.mean(xj)
-                        sparsity = np.mean(mask)
-                        covariance = np.mean((xi - mean_i) * (xj - mean_j))
-                        corr = covariance / (std_dev_i * std_dev_j)
-                        correlation_matrix[i, j] = corr * sparsity
+        for j in range(i + 1, n):  # Only compute upper triangle
+            mask = np.isfinite(data[i]) & np.isfinite(data[j])
+            if np.sum(mask) > 1:
+                xi = data[i][mask]
+                xj = data[j][mask]
+                std_dev_i = np.std(xi)
+                std_dev_j = np.std(xj)
+
+                if (std_dev_i > 0) and (std_dev_j > 0):
+                    mean_i = np.mean(xi)
+                    mean_j = np.mean(xj)
+                    sparsity = np.mean(mask)
+                    covariance = np.mean((xi - mean_i) * (xj - mean_j))
+                    corr = covariance / (std_dev_i * std_dev_j)
+                    correlation_matrix[i, j] = corr * sparsity
+                    correlation_matrix[j, i] = (
+                        corr * sparsity
+                    )  # Mirror to lower triangle
 
     return correlation_matrix
 
@@ -98,7 +111,7 @@ def _nan_correlation_w_ref(data_ref_pairs: tuple[np.ndarray, np.ndarray]) -> np.
         ref = ref[0]
 
     n = data.shape[0]
-    correlation_matrix = np.full(data.shape, np.nan)
+    correlation_matrix = np.zeros(data.shape)
     ref_mask = np.isfinite(ref)
 
     for i in range(n):  # Compute all elements
@@ -137,25 +150,59 @@ def feature_engineering_pipeline_for_unaligned_data(self, data):
         """
         logger.info("Feature engineering pipeline for unaligned data.")
 
-        total = self._feat_eng.parallel_process(
-            data, self._feat_eng.feature_engineering, func="sum", n_jobs=10, axis=1
+        ranks = self._feat_eng.parallel_process(
+            data,
+            self._feat_eng.feature_engineering,
+            func="rank_intensity",
+            n_jobs=10,
+            axis=0,
         )
-        return [total]
+        return [ranks]
+
+    def feature_engineering_pipeline_for_height_and_ms1_intensity(self, data):
+        logger.info("Calculating additional features.")
+        mean_corr, std_corr = zip(
+            *self._feat_eng.parallel_process(
+                data, self._feat_eng.feature_engineering, func="corr", n_jobs=10
+            )
+        )
+
+        return [
+            std_corr,
+            mean_corr,
+        ]
 
     def feature_engineering_pipeline_for_intensity(self, data):
         """
         Feature engineering pipeline for intensity data.
         """
-        logger.info("Calculating mean and std correlation across fragments.")
+        logger.info("Calculating additional features.")
         mean_corr, std_corr = zip(
             *self._feat_eng.parallel_process(
                 data, self._feat_eng.feature_engineering, func="corr", n_jobs=10
             )
         )
+        no_of_datapoints_across_fragments = self._feat_eng.parallel_process(
+            data,
+            self._feat_eng.feature_engineering,
+            func="count_no_of_datapoints",
+            n_jobs=10,
+            axis=0,
+        )
+
+        no_of_datapoints_across_samples = self._feat_eng.parallel_process(
+            data,
+            self._feat_eng.feature_engineering,
+            func="count_no_of_datapoints",
+            n_jobs=10,
+            axis=1,
+        )
 
         return [
             std_corr,
             mean_corr,
+            no_of_datapoints_across_fragments,
+            no_of_datapoints_across_samples,
         ]
 
 
@@ -184,9 +231,11 @@ def feature_engineering(self, data, func, **kwargs):
             "percentile": self._assign_percentiles,
             "mean_distance": _calculate_mean_distance,
             "var_distance": _calculate_variance_distance,
-            "corr": self._nan_mean_std_corr_across_fragments,
+            "corr": self.nan_mean_std_corr_across_fragments,
             "median_std_offset": self._calculate_median_std_offset,
             "derivative": self._calculate_derivative,
+            "count_no_of_datapoints": self._count_no_of_datapoints,
+            "weighted_variance": self._calculate_weighted_variance,
         }
 
         feat_eng = func_mapping[func]
@@ -200,10 +249,11 @@ def feature_engineering(self, data, func, **kwargs):
             "sum",
             "mad",
             "cv",
-            # "rank_intensity",
             "L2",
             "sparsity",
             "median_std_offset",
+            "count_no_of_datapoints",
+            "weighted_variance",
         ]:
             axis = kwargs.get("axis")
 
@@ -221,6 +271,12 @@ def _repeater(self, data, function, instance_method, **kwargs):
         else:
             return [function(subset, **kwargs) for subset in data]
 
+    def _calculate_weighted_variance(self, data, axis=0):
+        return np.nanvar(data, axis=axis) * self._sparsity(data, axis=axis)
+
+    def _count_no_of_datapoints(self, data, axis=0):
+        return np.sum(np.isfinite(data), axis=axis)
+
     def _calculate_median_std_offset(self, data, axis=1):
         stds = np.nanstd(data, axis=axis)
         median_stds = np.nanmedian(stds)
@@ -319,7 +375,7 @@ def _growth_decay_rate(self, data, axis=0):
             coefs = np.tile(coef, (imputed_data.shape[1], 1)).T
             return coefs
 
-    def _nan_mean_std_corr_across_fragments(self, data):
+    def nan_mean_std_corr_across_fragments(self, data):
         mean_corr = self._mean_tile(_nan_correlation_matrix(data), data.shape[1])
         std_corr = self._std_tile(_nan_correlation_matrix(data), data.shape[1])
 
@@ -331,12 +387,12 @@ def _nan_corrs(self, inputs):
         )
 
     def _mean_tile(self, data, samples):
-        mean = np.nanmean(data, axis=0)
+        mean = np.mean(data, axis=0)
         mean = np.tile(mean, (samples, 1)).T
         return mean
 
     def _std_tile(self, data, samples):
-        std_data = np.nanstd(data, axis=0)
+        std_data = np.std(data, axis=0)
         std_data = np.tile(std_data, (samples, 1)).T
         return std_data
 
@@ -373,3 +429,28 @@ def parallel_process(self, inputs, method, n_jobs=10, **kwargs):
         return Parallel(n_jobs=n_jobs)(
             delayed(method)(input, **kwargs) for input in inputs
         )
+
+    def calculate_ms1_ms2_corr(
+        self,
+        ms1_data_extracted: list[np.ndarray],
+        ms2_data_extracted: list[np.ndarray],
+    ) -> pd.DataFrame:
+        """
+        Calculate the correlation between the ms1 and ms2 data.
+
+        Parameters:
+            ms1_data_extracted: list[np.ndarray]
+                The ms1 data extracted from the precursor data.
+            ms2_data_extracted: list[np.ndarray]
+                The ms2 data extracted from the precursor data.
+        Returns:
+            ms1_ms2_corr_data: pd.DataFrame
+                The ms1-ms2 correlation data.
+        """
+
+        zipped_data = zip(ms2_data_extracted, ms1_data_extracted)
+        ms1_ms2_corr = self.parallel_process(
+            inputs=zipped_data, method=_nan_correlation_w_ref, n_jobs=10
+        )
+        ms1_ms2_corr = pd.DataFrame(np.vstack(ms1_ms2_corr))
+        return ms1_ms2_corr

selectlfq/ms1_features.py

@@ -203,31 +203,6 @@ class FeatureConfig:
     DEFAULT_FEATURES = [
         "intensity",
         "delta_rt",
-        "rt_observed",
-        "intensity_correlation",
-        "score",
-        "proba",
-        "base_width_rt",
-        "rt_calibrated",
-        "rt_library",
-        "delta_rt",
-        "cycle_fwhm",
-        "mz_observed",
-        "mz_library",
-        "mz_calibrated",
-        "mean_ms2_mass_error",
-        "top_3_ms2_mass_error",
-        "mean_overlapping_mass_error",
-        # "isotope_intensity_correlation",
-        # "isotope_height_correlation",
-        # "height_correlation",
-        # "fragment_scan_correlation",
-        # "template_scan_correlation",
-        # "fragment_frame_correlation",
-        "top3_frame_correlation",
-        # "template_frame_correlation",
-        "top3_b_ion_correlation",
-        # "top3_y_ion_correlation",
     ]
 
     NORMALIZATION_FEATURES = ["intensity", "mono_ms1_intensity", "sum_ms1_intensity"]