Remove mix-max dead code and --method CLI flag

crema/confidence.py

@@ -661,184 +661,6 @@ def _assign_confidence(self):
             self.decoy_confidence_estimates[level] = df.loc[~targets, :]
 
 
-class MixmaxConfidence(Confidence):
-    """Assign confidence estimates using mix-max competition.
-
-    Estimates qvalues using the mix-max competition method. To use this
-    method a separate target and decoy database search using a calibrated
-    score function must be used.
-
-    # TODO Describe how mixmax works here
-
-    Additional details can be found in this manuscript.
-    U. Keich, A. Kertesz-Farkas, and W. S. Noble. Improved false
-    discovery rate estimation procedure for shotgun proteomics.
-    Journal of Proteome Research, 14(8):3148-3161, 2015.
-
-    Parameters
-    ----------
-    psms : a PsmDataset object
-        A collection of PSMs
-    score_column : str, optional
-        The score by which to rank the PSMs for confidence estimation. If
-        :code:`None`, the score that yields the most PSMs at the specified
-        false discovery rate threshold (`eval_fdr`), will be used.
-    desc : bool, optional
-        True if higher scores better, False if lower scores are better.
-        If None, crema will try both and use the
-        choice that yields the most PSMs at the specified false discovery
-        rate threshold (`eval_fdr`). If `score_column` is :code:`None`,
-        this parameter is ignored.
-    eval_fdr : float, optional
-        The false discovery rate threshold used to evaluate the best
-        `score_column` and `desc` to choose. This should range from 0 to 1.
-    pep_fdr_type : {"psm-only","peptide-only",psm-peptide"}, optional
-        The method for crema to use when calculating peptide level confidence
-        estimates. Default is "psm-peptide".
-    prot_fdr_type : {"best", "combine"}, optional
-        The method for crema to use when calculating protein level confidence
-        estimates. Default is "best".
-        estimates.
-    threshold : float or "q-value", optional
-        The FDR threshold for accepting discoveries. Default is 0.01. If
-        "q-value" is chosen, then "accept" column is replaced with
-        "crema q-value".
-
-    Attributes
-    ----------
-    data : pandas.DataFrame
-    dataset : crema.PsmDataset
-    levels : list of str
-    confidence_estimates : Dict
-        A dictionary containing the confidence estimates at each level, each
-        as a :py:class:`pandas.DataFrame`.
-    decoy_confidence_estimates : Dict
-        A dictionary containing the confidence estimates for the decoy hits at
-        each level, each as a :py:class:`pandas.DataFrame`
-    """
-
-    def __init__(
-        self,
-        psms,
-        score_column=None,
-        desc=None,
-        eval_fdr=0.01,
-        pep_fdr_type="psm-peptide",
-        prot_fdr_type="best",
-        threshold=0.01,
-    ):
-        """Initialize a TdcConfidence object."""
-        LOGGER.info(
-            "Assigning confidence estimates using mix-max competition..."
-        )
-
-        super().__init__(
-            psms=psms,
-            score_column=score_column,
-            desc=desc,
-            eval_fdr=eval_fdr,
-            pep_fdr_type=pep_fdr_type,
-            prot_fdr_type=prot_fdr_type,
-        )
-
-    def _assign_confidence(self):
-        """Assign confidence estimates using target-decoy competition"""
-        # TODO maybe better way to do this
-        # can not infer desc value as the wrong value will
-        # result in a divide by zero error
-        if self._desc is None:
-            raise ValueError("'desc' has to be set for mix-max.")
-
-        # TODO check if separate target-decoy search is done
-        for level, group_cols in zip(self.levels, self._level_columns):
-            if level != "psms":
-                continue
-
-            df = self.data
-            group_cols = utils.listify(group_cols)
-
-            targets = df[df[self.dataset._target_column]]
-            decoys = df[~df[self.dataset._target_column]]
-
-            if targets.shape[1] != decoys.shape[1]:
-                LOGGER.warning(
-                    "The mix-max procedure is not well behaved when "
-                    "# targets (%i) != # decoys (%i).",
-                    targets.shape[0],
-                    decoys.shape[0],
-                )
-
-            if self._desc:
-                keep = "last"
-            else:
-                keep = "first"
-
-            # sort targets by score column and keep top rank
-            targets_sorted = (
-                targets.sample(frac=1)
-                .sort_values([self._score_column, *group_cols])
-                .drop_duplicates(group_cols, keep=keep, ignore_index=True)
-            )
-
-            # sort decoys by score column and keep top rank
-            decoys_sorted = (
-                decoys.sample(frac=1)
-                .sort_values([self._score_column, *group_cols])
-                .drop_duplicates(group_cols, keep=keep, ignore_index=True)
-            )
-
-            # combine top ranked target and decoy into one dataframe
-            combined = pd.concat([targets_sorted, decoys_sorted])
-            combined_sorted = combined.sample(frac=1).sort_values(
-                [self._score_column], ascending=~self._desc, ignore_index=True
-            )
-
-            # qvalues.py::calculate_mixmax_qval expects target scores
-            # and decoy scores to be sorted from worst to best
-            # qvalues.py::mixmax expected combined_sorted scores
-            # to be sorted from best to worst
-            if self._desc:  # larger score is better
-                combined_sorted = combined_sorted[::-1]
-            else:  # smaller score is better
-                targets_sorted[self._score_column] = (
-                    targets_sorted[self._score_column] * -1.0
-                )
-                decoys_sorted[self._score_column] = (
-                    decoys_sorted[self._score_column] * -1.0
-                )
-                targets_sorted = targets_sorted[::-1]
-                decoys_sorted = decoys_sorted[::-1]
-
-            # Now calculate q-values:
-            pi0, targets_sorted["crema q-value"] = qvalues.mixmax(
-                target_scores=targets_sorted[self._score_column],
-                decoy_scores=decoys_sorted[self._score_column],
-                combined_score=combined_sorted[self._score_column],
-                combined_score_target=combined_sorted[
-                    self.dataset._target_column
-                ],
-            )
-
-            LOGGER.info("  - Estimated pi_zero = %f.", pi0)
-
-            LOGGER.info(
-                "  - Found %i %s at q<=%g.",
-                (targets_sorted["crema q-value"] <= self._eval_fdr).sum(),
-                self._level_labs[level],
-                self._eval_fdr,
-            )
-
-            # reverse rows so that best score is at top
-            targets_sorted = targets_sorted[::-1]
-
-            # undo previous multipliation by -1.0
-            if not self._desc:
-                targets_sorted[self._score_column] = (
-                    targets_sorted[self._score_column] * -1.0
-                )
-            self.confidence_estimates[level] = targets_sorted
-
-
 class DuckdbTdcConfidence(Confidence):
     """TDC confidence estimation using DuckDB as the compute engine.
 

crema/crema.py

@@ -128,7 +128,6 @@ def _run_assign_confidence(args, start_time):
         prot_fdr_type=args.prot_fdr_type,
         desc=desc,
         eval_fdr=args.eval_fdr,
-        method=args.method,
     )
 
     logging.info("Writing results...")

crema/dataset.py

@@ -3,10 +3,8 @@
 """
 
 import logging
-from typing import ClassVar
 
 from .confidence import TdcConfidence
-from .confidence import MixmaxConfidence
 from .confidence import DuckdbTdcConfidence
 from .qvalues import tdc
 from .utils import listify
@@ -59,15 +57,9 @@ class PsmDataset:
     peptide_column : str
     protein_column : str
     protein_delim : str
-    methods : dict
     peptide_pairing : dict
     """
 
-    methods: ClassVar[dict] = {
-        "tdc": TdcConfidence,
-        "mixmax": MixmaxConfidence,
-    }
-
     def __init__(
         self,
         psms,
@@ -171,7 +163,6 @@ def assign_confidence(
         prot_fdr_type="best",
         desc=None,
         eval_fdr=0.01,
-        method="tdc",
         backend="pandas",
     ):
         """Assign confidence estimates to this collection of peptide-spectrum matches.
@@ -201,11 +192,9 @@ def assign_confidence(
         eval_fdr : float, optional
             The false discovery rate threshold used to evaluate the best
             `score_column` and `desc` to choose. This should range from 0 to 1.
-        method : {"tdc"}, optional
-            The method for crema to use when calculating the confidence estimates.
         backend : {"pandas", "duckdb"}, optional
             The compute backend. ``"duckdb"`` requires the optional ``duckdb``
-            package and only supports ``method="tdc"``.
+            package.
 
         Returns
         -------
@@ -216,10 +205,6 @@ def assign_confidence(
             score_column, _, desc = self.find_best_score(eval_fdr)
 
         if backend == "duckdb":
-            if method != "tdc":
-                raise ValueError(
-                    "backend='duckdb' only supports method='tdc'."
-                )
             conf = DuckdbTdcConfidence(
                 psms=self,
                 score_column=score_column,
@@ -230,7 +215,7 @@ def assign_confidence(
                 threshold=threshold,
             )
         elif backend == "pandas":
-            conf = self.methods[method](
+            conf = TdcConfidence(
                 psms=self,
                 score_column=score_column,
                 desc=desc,

crema/params.py

@@ -172,15 +172,6 @@ def _configure_parser():
         ),
     )
 
-    ac.add_argument(
-        "-m",
-        "--method",
-        type=str,
-        default="tdc",
-        choices=["tdc"],
-        help="The confidence estimation method to use.",
-    )
-
     ac.add_argument(
         "--parquet",
         action="store_true",

crema/parsers/msamanda.py

@@ -54,7 +54,7 @@ def read_msamanda(
     scores_all = scores
 
     # Keep only MSAmanda scores that exist in all of the files.
-    skip_first_row = False
+    skip_first_row = 0
     if isinstance(txt_files, pd.DataFrame):
         scores = scores.intersection(set(txt_files.columns))
     else:
@@ -65,12 +65,15 @@ def read_msamanda(
                 raise ValueError(f"{txt_file} must be in .csv format.")
 
             with open(txt_file) as txt_ref:
-                # First line of MSAmanda output consists only of version line
-                # If statement below in case first line is removed
+                # MSAmanda output may be preceded by one or more comment
+                # lines (e.g. a version line) starting with '#'. Skip all
+                # of them to find the header line.
+                n_comment_lines = 0
                 line = txt_ref.readline().rstrip()
-                if line.startswith("#version"):
+                while line.startswith("#"):
+                    n_comment_lines += 1
                     line = txt_ref.readline().rstrip()
-                    skip_first_row = True
+                skip_first_row = n_comment_lines
                 cols = line.split("\t")
                 scores = scores.intersection(set(cols))