Continous Integration Tests

.github/workflows/ci.yml

@@ -0,0 +1,31 @@
+name: continuous-integration
+
+on: [push]
+
+jobs:
+  test:
+    runs-on: ${{ matrix.os }}
+    strategy:
+      matrix:
+        os: [ubuntu-latest]
+        # Requirements file generated with python=3.11 
+        python-version: ["3.11"]
+    steps:
+    - uses: actions/checkout@v4
+
+    - name: Set up Python ${{ matrix.python-version }}
+      uses: actions/setup-python@v5
+      with:
+        python-version: ${{ matrix.python-version }}
+    - name: Install dependencies
+      run: |
+        python -m pip install --upgrade pip
+        pip install -r requirements.txt # test with requirements file so can easily bump with dependabot
+        pip install .
+
+    - name: Compile cython module
+      run: python setup.py build_ext --inplace 
+
+    - name: Test
+      run: |
+        python -m pytest tests/ 

.github/workflows/dependabot.yml

@@ -0,0 +1,9 @@
+version: 2
+updates:
+  - package-ecosystem: "pip"
+    directory: "/"  # Location of your pyproject.toml or requirements.txt
+    schedule:
+      interval: "weekly"  # Checks for updates every week
+    commit-message:
+      prefix: "deps"  # Prefix for pull request titles
+    open-pull-requests-limit: 5  # Limit the number of open PRs at a time

.gitignore

@@ -36,3 +36,4 @@ nosetests.xml
 
 # vim
 *.sw[opqrs]
+*~

pyproject.toml

@@ -0,0 +1,55 @@
+[build-system]
+requires = ["setuptools", "wheel", "numpy", "cython"]  # Dependencies needed to build the package
+build-backend = "setuptools.build_meta"
+
+[project]
+name = "pyprophet"
+version = "2.2.9"
+description = "PyProphet: Semi-supervised learning and scoring of OpenSWATH results."
+readme = { file = "README.md", content-type = "text/markdown" }
+license = { text = "BSD" }
+authors = [{ name = "The PyProphet Developers", email = "[email protected]" }]
+classifiers = [
+    "Development Status :: 3 - Alpha",
+    "Environment :: Console",
+    "Intended Audience :: Science/Research",
+    "License :: OSI Approved :: BSD License",
+    "Operating System :: OS Independent",
+    "Topic :: Scientific/Engineering :: Bio-Informatics",
+    "Topic :: Scientific/Engineering :: Chemistry"
+]
+keywords = ["bioinformatics", "openSWATH", "mass spectrometry"]
+
+# Dependencies required for runtime
+dependencies = [
+    "Click",
+    "duckdb",
+    "duckdb-extensions",
+    "duckdb-extension-sqlite-scanner",
+    "numpy >= 2.0",
+    "scipy",
+    "pandas >= 0.17",
+    "cython",
+    "numexpr >= 2.10.1",
+    "scikit-learn >= 0.17",
+    "xgboost",
+    "hyperopt",
+    "statsmodels >= 0.8.0",
+    "matplotlib",
+    "tabulate",
+    "pyarrow",
+    "pypdf"
+]
+
+# Optional dependencies 
+[project.optional-dependencies]
+testing = ["pytest", "pytest-regtest"]
+
+# Define console entry points
+[project.scripts]
+pyprophet = "pyprophet.main:cli"
+
+[tool.setuptools]
+packages = { find = { exclude = ["ez_setup", "examples", "tests"] } }
+include-package-data = true
+zip-safe = false

pyprophet/classifiers.py

@@ -110,7 +110,7 @@ def objective(params):
 
             clf = xgb.XGBClassifier(random_state=42, verbosity=0, objective='binary:logitraw', eval_metric='auc', **params)
 
-            score = cross_val_score(clf, X, y, scoring='roc_auc', n_jobs=self.threads, cv=KFold(n_splits=3, shuffle=True, random_state=np.random.RandomState(42))).mean()
+            score = cross_val_score(clf, X, y, scoring='roc_auc', n_jobs=self.threads, cv=KFold(n_splits=3, shuffle=True, random_state=42)).mean()
             # click.echo("Info: AUC: {:.3f} hyperparameters: {}".format(score, params))
             return score
 
@@ -129,7 +129,8 @@ def objective(params):
         xgb_params_complexity = self.xgb_params_tuned
         xgb_params_complexity.update({k: self.xgb_params_space[k] for k in ('max_depth', 'min_child_weight')})
 
-        best_complexity = fmin(fn=objective, space=xgb_params_complexity, algo=tpe.suggest, max_evals=self.xgb_hyperparams['autotune_num_rounds'], rstate=np.random.RandomState(42))
+        rng = np.random.default_rng(42)
+        best_complexity = fmin(fn=objective, space=xgb_params_complexity, algo=tpe.suggest, max_evals=self.xgb_hyperparams['autotune_num_rounds'], rstate=rng)
         best_complexity['max_depth'] = int(best_complexity['max_depth'])
         best_complexity['min_child_weight'] = int(best_complexity['min_child_weight'])
 
@@ -139,31 +140,31 @@ def objective(params):
         xgb_params_gamma = self.xgb_params_tuned
         xgb_params_gamma['gamma'] = self.xgb_params_space['gamma']
 
-        best_gamma = fmin(fn=objective, space=xgb_params_gamma, algo=tpe.suggest, max_evals=self.xgb_hyperparams['autotune_num_rounds'], rstate=np.random.RandomState(42))
+        best_gamma = fmin(fn=objective, space=xgb_params_gamma, algo=tpe.suggest, max_evals=self.xgb_hyperparams['autotune_num_rounds'], rstate=rng)
 
         self.xgb_params_tuned.update(best_gamma)
 
         # Tune subsampling hyperparameters
         xgb_params_subsampling = self.xgb_params_tuned
         xgb_params_subsampling.update({k: self.xgb_params_space[k] for k in ('subsample', 'colsample_bytree', 'colsample_bylevel', 'colsample_bynode')})
 
-        best_subsampling = fmin(fn=objective, space=xgb_params_subsampling, algo=tpe.suggest, max_evals=self.xgb_hyperparams['autotune_num_rounds'], rstate=np.random.RandomState(42))
+        best_subsampling = fmin(fn=objective, space=xgb_params_subsampling, algo=tpe.suggest, max_evals=self.xgb_hyperparams['autotune_num_rounds'], rstate=rng)
 
         self.xgb_params_tuned.update(best_subsampling)
 
         # Tune regularization hyperparameters
         xgb_params_regularization = self.xgb_params_tuned
         xgb_params_regularization.update({k: self.xgb_params_space[k] for k in ('lambda', 'alpha')})
 
-        best_regularization = fmin(fn=objective, space=xgb_params_regularization, algo=tpe.suggest, max_evals=self.xgb_hyperparams['autotune_num_rounds'], rstate=np.random.RandomState(42))
+        best_regularization = fmin(fn=objective, space=xgb_params_regularization, algo=tpe.suggest, max_evals=self.xgb_hyperparams['autotune_num_rounds'], rstate=rng)
 
         self.xgb_params_tuned.update(best_regularization)
 
         # Tune learning rate
         xgb_params_learning = self.xgb_params_tuned
         xgb_params_learning['eta'] = self.xgb_params_space['eta']
 
-        best_learning = fmin(fn=objective, space=xgb_params_learning, algo=tpe.suggest, max_evals=self.xgb_hyperparams['autotune_num_rounds'], rstate=np.random.RandomState(42))
+        best_learning = fmin(fn=objective, space=xgb_params_learning, algo=tpe.suggest, max_evals=self.xgb_hyperparams['autotune_num_rounds'], rstate=rng)
 
         self.xgb_params_tuned.update(best_learning)
         click.echo("Info: Optimal hyperparameters: {}".format(self.xgb_params_tuned))

pyprophet/data_handling.py

@@ -5,6 +5,7 @@
 import sys
 import os
 import multiprocessing
+from .stats import mean_and_std_dev 
 
 from .optimized import find_top_ranked, rank
 
@@ -336,6 +337,21 @@ def get_top_target_peaks(self):
     def get_feature_matrix(self, use_main_score):
         min_col = 5 if use_main_score else 6
         return self.df.iloc[:, min_col:-1].values
+
+    def normalize_score_by_decoys(self, score_col_name):
+        '''
+        normalize the decoy scores to mean 0 and std 1, scale the targets accordingly
+        Args:
+            score_col_name: str, the name of the score column
+        '''
+        td_scores = self.get_top_decoy_peaks()[score_col_name]
+        mu, nu = mean_and_std_dev(td_scores)
+
+        if nu == 0:
+            raise Exception("Warning: Standard deviation of decoy scores is zero. Cannot normalize scores.")
+
+        self.df.loc[:, score_col_name] = (self.df[score_col_name] - mu) / nu
+
 
     def filter_(self, idx):
         return Experiment(self.df[idx])
@@ -344,7 +360,7 @@ def filter_(self, idx):
     def add_peak_group_rank(self):
         ids = self.df.tg_num_id.values
         scores = self.df.d_score.values
-        peak_group_ranks = rank(ids, scores)
+        peak_group_ranks = rank(ids, scores.astype(np.float32, copy=False))
         self.df["peak_group_rank"] = peak_group_ranks
 
     @profile

pyprophet/export_parquet.py

@@ -172,7 +172,7 @@ def export_to_parquet(infile, outfile, transitionLevel, onlyFeatures=False):
 
     # transition level
     if transitionLevel:
-        columns['FEATURE_TRANSITION'] = ['AREA_INTENSITY', 'TOTAL_AREA_INTENSITY', 'APEX_INTENSITY', 'TOTAL_MI'] + getVarColumnNames(condb, 'FEATURE_TRANSITION')
+        columns['FEATURE_TRANSITION'] = ['AREA_INTENSITY', 'TOTAL_AREA_INTENSITY', 'APEX_INTENSITY', 'TOTAL_MI'] + getVarColumnNames(con, 'FEATURE_TRANSITION')
         columns['TRANSITION'] = ['TRAML_ID', 'PRODUCT_MZ', 'CHARGE', 'TYPE', 'ORDINAL', 'DETECTING', 'IDENTIFYING', 'QUANTIFYING', 'LIBRARY_INTENSITY']
         columns['TRANSITION_PRECURSOR_MAPPING'] = ['TRANSITION_ID']
 

pyprophet/levels_contexts.py

@@ -33,7 +33,12 @@ def statistics_report(data, outfile, context, analyte, parametric, pfdr, pi0_lam
         outfile = outfile + "_" + str(data['run_id'].unique()[0])
 
     # export PDF report
-    save_report(outfile + "_" + context + "_" + analyte + ".pdf", outfile + ": " + context + " " + analyte + "-level error-rate control", data[data.decoy==1]["score"], data[data.decoy==0]["score"], stat_table["cutoff"], stat_table["svalue"], stat_table["qvalue"], data[data.decoy==0]["p_value"], pi0, color_palette)
+    save_report(outfile + "_" + context + "_" + analyte + ".pdf", 
+                outfile + ": " + context + " " + analyte + "-level error-rate control", 
+                data[data.decoy==1]["score"].values, data[data.decoy==0]["score"].values, stat_table["cutoff"].values, 
+                stat_table["svalue"].values, stat_table["qvalue"].values, data[data.decoy==0]["p_value"].values, 
+                pi0, 
+                color_palette)
 
     return(data)
 
@@ -184,7 +189,7 @@ def infer_proteins(infile, outfile, context, parametric, pfdr, pi0_lambda, pi0_m
     con.close()
 
     if context == 'run-specific':
-        data = data.groupby('run_id').apply(statistics_report, outfile, context, "protein", parametric, pfdr, pi0_lambda, pi0_method, pi0_smooth_df, pi0_smooth_log_pi0, lfdr_truncate, lfdr_monotone, lfdr_transformation, lfdr_adj, lfdr_eps, color_palette).reset_index()
+        data = data.groupby('run_id').apply(statistics_report, outfile, context, "protein", parametric, pfdr, pi0_lambda, pi0_method, pi0_smooth_df, pi0_smooth_log_pi0, lfdr_truncate, lfdr_monotone, lfdr_transformation, lfdr_adj, lfdr_eps, color_palette)
 
     elif context in ['global', 'experiment-wide']:
         data = statistics_report(data, outfile, context, "protein", parametric, pfdr, pi0_lambda, pi0_method, pi0_smooth_df, pi0_smooth_log_pi0, lfdr_truncate, lfdr_monotone, lfdr_transformation, lfdr_adj, lfdr_eps, color_palette)
@@ -257,7 +262,7 @@ def infer_peptides(infile, outfile, context, parametric, pfdr, pi0_lambda, pi0_m
     con.close()
 
     if context == 'run-specific':
-        data = data.groupby('run_id').apply(statistics_report, outfile, context, "peptide", parametric, pfdr, pi0_lambda, pi0_method, pi0_smooth_df, pi0_smooth_log_pi0, lfdr_truncate, lfdr_monotone, lfdr_transformation, lfdr_adj, lfdr_eps, color_palette).reset_index()
+        data = data.groupby('run_id').apply(statistics_report, outfile, context, "peptide", parametric, pfdr, pi0_lambda, pi0_method, pi0_smooth_df, pi0_smooth_log_pi0, lfdr_truncate, lfdr_monotone, lfdr_transformation, lfdr_adj, lfdr_eps, color_palette)
 
     elif context in ['global', 'experiment-wide']:
         data = statistics_report(data, outfile, context, "peptide", parametric, pfdr, pi0_lambda, pi0_method, pi0_smooth_df, pi0_smooth_log_pi0, lfdr_truncate, lfdr_monotone, lfdr_transformation, lfdr_adj, lfdr_eps, color_palette)

pyprophet/main.py

@@ -106,6 +106,8 @@ def score(infile, outfile, classifier, xgb_autotune, apply_weights, xeval_fracti
     xgb_hyperparams = {'autotune': xgb_autotune, 'autotune_num_rounds': 10, 'num_boost_round': 100, 'early_stopping_rounds': 10, 'test_size': 0.33}
 
     xgb_params = {'eta': 0.3, 'gamma': 0, 'max_depth': 6, 'min_child_weight': 1, 'subsample': 1, 'colsample_bytree': 1, 'colsample_bylevel': 1, 'colsample_bynode': 1, 'lambda': 1, 'alpha': 0, 'scale_pos_weight': 1, 'verbosity': 0, 'objective': 'binary:logitraw', 'nthread': 1, 'eval_metric': 'auc'}
+    if test:
+        xgb_params['tree_method'] = 'exact'
 
     xgb_params_space = {'eta': hp.uniform('eta', 0.0, 0.3), 'gamma': hp.uniform('gamma', 0.0, 0.5), 'max_depth': hp.quniform('max_depth', 2, 8, 1), 'min_child_weight': hp.quniform('min_child_weight', 1, 5, 1), 'subsample': 1, 'colsample_bytree': 1, 'colsample_bylevel': 1, 'colsample_bynode': 1, 'lambda': hp.uniform('lambda', 0.0, 1.0), 'alpha': hp.uniform('alpha', 0.0, 1.0), 'scale_pos_weight': 1.0, 'verbosity': 0, 'objective': 'binary:logitraw', 'nthread': 1, 'eval_metric': 'auc'}
 

pyprophet/semi_supervised.py

@@ -3,7 +3,7 @@
 
 from .data_handling import Experiment, update_chosen_main_score_in_table
 from .classifiers import AbstractLearner, XGBLearner
-from .stats import mean_and_std_dev, find_cutoff
+from .stats import find_cutoff
 
 try:
     profile
@@ -64,13 +64,9 @@ def learn_randomized(self, experiment, score_columns, working_thread_number):
 
         # after semi supervised iteration: classify full dataset
         clf_scores = self.score(experiment, params)
-        mu, nu = mean_and_std_dev(clf_scores)
         experiment.set_and_rerank("classifier_score", clf_scores)
 
-        td_scores = experiment.get_top_decoy_peaks()["classifier_score"]
-
-        mu, nu = mean_and_std_dev(td_scores)
-        experiment["classifier_score"] = (experiment["classifier_score"] - mu) / nu
+        experiment.normalize_score_by_decoys('classifier_score')
         experiment.rank_by("classifier_score")
 
         top_test_peaks = experiment.get_top_test_peaks()
@@ -92,13 +88,9 @@ def learn_final(self, experiment):
 
         # after semi supervised iteration: classify full dataset
         clf_scores = self.score(experiment, params)
-        mu, nu = mean_and_std_dev(clf_scores)
         experiment.set_and_rerank("classifier_score", clf_scores)
 
-        td_scores = experiment.get_top_decoy_peaks()["classifier_score"]
-
-        mu, nu = mean_and_std_dev(td_scores)
-        experiment["classifier_score"] = (experiment["classifier_score"] - mu) / nu
+        experiment.normalize_score_by_decoys('classifier_score')
         experiment.rank_by("classifier_score")
 
         return params