Continous Integration Tests

.github/workflows/ci.yml

@@ -0,0 +1,32 @@
+name: continuous-integration
+
+on: [push]
+
+jobs:
+  test:
+    runs-on: ${{ matrix.os }}
+    strategy:
+      matrix:
+        os: [ubuntu-latest]
+          #os: [ubuntu-latest, windows-latest, macos-latest] # remove mac tests
+        # 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

pyproject.toml

@@ -0,0 +1,51 @@
+[build-system]
+requires = ["setuptools", "wheel", "numpy", "cython"]  # Dependencies needed to build the package
+build-backend = "setuptools.build_meta"
+
+[project]
+name = "pyprophet"
+version = "2.2.8"
+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 = "rocksportrocker@gmail.com" }]
+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 >= 1.9.0",
+    "scipy",
+    "pandas >= 0.17",
+    "cython",
+    "numexpr >= 2.10.1",
+    "scikit-learn >= 0.17",
+    "xgboost",
+    "hyperopt",
+    "statsmodels >= 0.8.0",
+    "matplotlib",
+    "tabulate",
+    "pyarrow",
+    "pypdf"
+]
+
+# 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,8 @@ 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()
+            rng = np.random.default_rng(42)
+            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 +130,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 +141,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/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/stats.py

@@ -118,7 +118,10 @@ def posterior_chromatogram_hypotheses_fast(experiment, prior_chrom_null):
 
 
 def mean_and_std_dev(values):
-    return np.mean(values), np.std(values, ddof=1)
+    std = np.std(values, ddof=1)
+    if std == 0:
+        raise RuntimeError("Computed standard deviation is 0, cannot perform normalization")
+    return np.mean(values), std
 
 
 def pnorm(stat, stat0):
@@ -233,7 +236,7 @@ def pi0est(p_values, lambda_ = np.arange(0.05,1.0,0.05), pi0_method = "smoother"
 
 @profile
 def qvalue(p_values, pi0, pfdr = False):
-    p = np.array(p_values)
+    p = np.array(p_values).copy()
 
     qvals_out = p
     rm_na = np.isfinite(p)
@@ -277,7 +280,7 @@ def bw_nrd0(x):
 @profile
 def lfdr(p_values, pi0, trunc = True, monotone = True, transf = "probit", adj = 1.5, eps = np.power(10.0,-8)):
     """ Estimate local FDR / posterior error probability from p-values according to bioconductor/qvalue """
-    p = np.array(p_values)
+    p = np.array(p_values).copy()
 
     # Compare to bioconductor/qvalue reference implementation
     # import rpy2

requirements.txt

@@ -0,0 +1,118 @@
+#
+# This file is autogenerated by pip-compile with Python 3.11
+# by the following command:
+#
+#    pip-compile --all-extras --output-file=requirements.txt
+#
+click==8.1.7
+    # via pyprophet (setup.py)
+cloudpickle==3.1.0
+    # via hyperopt
+contourpy==1.3.0
+    # via matplotlib
+cycler==0.12.1
+    # via matplotlib
+cython==3.0.11
+    # via pyprophet (setup.py)
+duckdb==1.1.3
+    # via
+    #   duckdb-extension-sqlite-scanner
+    #   duckdb-extensions
+    #   pyprophet (setup.py)
+duckdb-extension-sqlite-scanner==1.1.3
+    # via pyprophet (setup.py)
+duckdb-extensions==1.1.3
+    # via pyprophet (setup.py)
+fonttools==4.55.0
+    # via matplotlib
+future==1.0.0
+    # via hyperopt
+hyperopt==0.2.7
+    # via pyprophet (setup.py)
+iniconfig==2.0.0
+    # via pytest
+joblib==1.4.2
+    # via scikit-learn
+kiwisolver==1.4.7
+    # via matplotlib
+matplotlib==3.9.2
+    # via pyprophet (setup.py)
+networkx==3.2.1
+    # via hyperopt
+numexpr==2.10.1
+    # via pyprophet (setup.py)
+numpy==2.0.2
+    # via
+    #   contourpy
+    #   hyperopt
+    #   matplotlib
+    #   numexpr
+    #   pandas
+    #   patsy
+    #   pyprophet (setup.py)
+    #   scikit-learn
+    #   scipy
+    #   statsmodels
+    #   xgboost
+nvidia-nccl-cu12==2.23.4
+    # via xgboost
+packaging==24.2
+    # via
+    #   matplotlib
+    #   pytest
+    #   statsmodels
+pandas==2.2.3
+    # via
+    #   pyprophet (setup.py)
+    #   statsmodels
+patsy==1.0.1
+    # via statsmodels
+pillow==11.0.0
+    # via matplotlib
+pluggy==1.5.0
+    # via pytest
+py4j==0.10.9.7
+    # via hyperopt
+pyarrow==18.0.0
+    # via pyprophet (setup.py)
+pyparsing==3.2.0
+    # via matplotlib
+pypdf==5.1.0
+    # via pyprophet (setup.py)
+pytest==8.3.3
+    # via
+    #   pyprophet (setup.py)
+    #   pytest-regtest
+pytest-regtest==2.3.3
+    # via pyprophet (setup.py)
+python-dateutil==2.9.0.post0
+    # via
+    #   matplotlib
+    #   pandas
+pytz==2024.2
+    # via pandas
+scikit-learn==1.5.2
+    # via pyprophet (setup.py)
+scipy==1.13.1
+    # via
+    #   hyperopt
+    #   pyprophet (setup.py)
+    #   scikit-learn
+    #   statsmodels
+    #   xgboost
+six==1.16.0
+    # via
+    #   hyperopt
+    #   python-dateutil
+statsmodels==0.14.4
+    # via pyprophet (setup.py)
+tabulate==0.9.0
+    # via pyprophet (setup.py)
+threadpoolctl==3.5.0
+    # via scikit-learn
+tqdm==4.67.0
+    # via hyperopt
+tzdata==2024.2
+    # via pandas
+xgboost==2.1.2
+    # via pyprophet (setup.py)