diff --git a/onsite/ascore/cli.py b/onsite/ascore/cli.py index 6a354a6..108f0ea 100644 --- a/onsite/ascore/cli.py +++ b/onsite/ascore/cli.py @@ -204,10 +204,26 @@ def ascore( else: stats["errors"] += 1 + # Count how many PSMs had their modification sites reassigned + relocated_count = 0 + orig_pf_map = {} + for _, group_df in grouped: + pep_idx = group_df.iloc[0].get("peptide_identification_index") + orig_pf = str(group_df.iloc[0].get("peptidoform", "")) + if orig_pf: + orig_pf_map[pep_idx] = orig_pf + for row in result_rows: + pep_idx = row.get("peptide_identification_index") + orig_pf = orig_pf_map.get(pep_idx, "") + out_pf = str(row.get("peptidoform", "")) + if orig_pf and orig_pf != out_pf: + relocated_count += 1 + elapsed = time.time() - start_time - click.echo(f"\nProcessing Complete:") + click.echo("\nProcessing Complete:") click.echo(f" Total identifications: {stats['total']}") click.echo(f" Successfully processed: {stats['processed']}") + click.echo(f" Modification sites reassigned: {relocated_count}") click.echo(f" Phosphorylated peptides: {stats['phospho']}") click.echo(f" Processing errors: {stats['errors']}") click.echo(f" Time elapsed: {elapsed:.2f} seconds") diff --git a/onsite/idparquet.py b/onsite/idparquet.py index 1a92dcc..7a67b4b 100644 --- a/onsite/idparquet.py +++ b/onsite/idparquet.py @@ -265,7 +265,7 @@ def _pa_col(name, values): return pa.array(values, type=tbl.schema.field(name).type) updates = {} - for col in ("peptidoform", "sequence", "score", "score_type", "higher_score_better"): + for col in ("peptidoform",): updates[col] = _pa_col(col, psms_df[col].tolist()) updates["psm_metavalues"] = _pa_col("psm_metavalues", psms_df["psm_metavalues"].tolist()) diff --git a/onsite/lucxor/cli.py b/onsite/lucxor/cli.py index a26fa5d..73563c3 100644 --- a/onsite/lucxor/cli.py +++ b/onsite/lucxor/cli.py @@ -21,8 +21,7 @@ MSExperiment, SpectrumLookup, AASequence, - PeptideHit, - PeptideIdentification + PeptideHit ) from onsite.idparquet import pyopenms_to_unimod_notation, save_dataframes, load_dataframes, unimod_to_pyopenms_notation @@ -245,6 +244,7 @@ def lucxor( debug=debug, add_decoys=add_decoys, fragment_method=fragment_method, + modeling_score_threshold=modeling_score_threshold, ) try: @@ -747,7 +747,21 @@ def run( # 6. Build output DataFrame (using second round calculation results) result_rows = [] - phospho_count = 0 + relocated_count = 0 + + # Pre-filter hit_index == 0 rows once, outside the loop + _lucxor_df = getattr(self, "_psms_df_template", None) + _input_by_ref = {} + if _lucxor_df is not None: + _hit0 = _lucxor_df[_lucxor_df["hit_index"] == 0] + for _, row in _hit0.iterrows(): + ref = str(row.get("spectrum_reference", "")) + if ref: + _input_by_ref[ref] = { + "peptidoform": str(row.get("peptidoform", "")), + "psm_metavalues": row.get("psm_metavalues"), + } + for pep_idx, psm in enumerate(all_psms): # Use metadata stashed during load_input_files mz = getattr(psm, "_mz", 0.0) @@ -775,6 +789,20 @@ def run( peptidoform = pyopenms_to_unimod_notation(seq_str).upper() + # Build psm_metavalues: preserve original + add Luciphor scores + _in_rec = _input_by_ref.get(spec_ref) + if _in_rec is not None: + _orig_mv = _in_rec["psm_metavalues"] + _orig_list = list(_orig_mv) if isinstance(_orig_mv, np.ndarray) else [] + else: + _orig_list = [] + + # Track how many peptides had their modification positions changed + if _in_rec is not None: + orig_peptidoform = _in_rec["peptidoform"] + if orig_peptidoform and orig_peptidoform != peptidoform: + relocated_count += 1 + # Use pyOpenMS to get unmodified string if possible try: seq_obj = AASequence.fromString(seq_str) if seq_str else None @@ -786,17 +814,6 @@ def run( # Fallback: strip all (...) modification annotations and uppercase base_seq = re.sub(r"\([^)]*\)", "", seq_str).upper() - # Build psm_metavalues: preserve original + add Luciphor scores - _lucxor_df = getattr(self, "_psms_df_template", None) - if _lucxor_df is not None: - _best = _lucxor_df[_lucxor_df["hit_index"] == 0] - if pep_idx < len(_best): - _orig_mv = _best.iloc[pep_idx].get("psm_metavalues") - _orig_list = list(_orig_mv) if isinstance(_orig_mv, np.ndarray) else [] - else: - _orig_list = [] - else: - _orig_list = [] _lucxor_managed = {"search_engine_sequence", "Luciphor_pep_score", "Luciphor_global_flr", "Luciphor_local_flr", "Luciphor_site_scores"} _filtered_orig = [m for m in _orig_list if isinstance(m, dict) and m.get("name") not in _lucxor_managed] @@ -844,12 +861,6 @@ def run( "run_identifier": "", }) - try: - if "(Phospho)" in seq_str: - phospho_count += 1 - except Exception: - pass - # 7. Save results (idParquet format) out_df = pd.DataFrame(result_rows) if result_rows else pd.DataFrame() save_dataframes(output, out_df, prot_ids, template_df=self._psms_df_template, source_idparquet=input_id) @@ -864,7 +875,7 @@ def run( print("\nProcessing Complete:") print(f" Total identifications: {total}") print(f" Successfully processed: {processed}") - print(f" Phosphorylated peptides: {phospho_count}") + print(f" Relocated: {relocated_count}") print(f" Processing errors: {errors}") print(f" Time elapsed: {elapsed:.2f} seconds") if elapsed > 0: @@ -875,7 +886,7 @@ def run( self.logger.info({ "total": total, "processed": processed, - "phospho": phospho_count, + "reassigned": relocated_count, "errors": errors, "elapsed_sec": round(elapsed, 2), "speed_ids_per_sec": round(processed/elapsed, 2) if elapsed > 0 else None diff --git a/onsite/onsitec.py b/onsite/onsitec.py index f467771..4f6ab25 100644 --- a/onsite/onsitec.py +++ b/onsite/onsitec.py @@ -20,7 +20,7 @@ from onsite.idparquet import load_dataframes, save_dataframes, unimod_to_pyopenms_notation @click.group() -@click.version_option(version="0.0.3") +@click.version_option(version="0.0.4") def cli(): """ OnSite: Mass spectrometry post-translational modification localization tool @@ -108,6 +108,13 @@ def cli(): default=False, help="Include A (PhosphoDecoy) as potential phosphorylation site for all algorithms", ) +@click.option( + "--modeling-score-threshold", + "modeling_score_threshold", + type=float, + default=0.95, + help="Score threshold for LucXor model training (default: 0.95, auto-adjusted based on score type)", +) def all( in_file, id_file, @@ -118,6 +125,7 @@ def all( threads, debug, add_decoys, + modeling_score_threshold, ): """Run all three algorithms (AScore, PhosphoRS, LucXor) and merge results.""" try: @@ -164,17 +172,19 @@ def all( click.echo(f"\n{'='*60}") click.echo(f"[{time.strftime('%H:%M:%S')}] Running LucXor...") click.echo(f"{'='*60}") - from onsite.lucxor.cli import lucxor as lucxor_func + from onsite.lucxor.cli import PyLuciPHOr2, setup_logging as lucxor_setup_logging if add_decoys: target_mods = ("Phospho (S)", "Phospho (T)", "Phospho (Y)", "PhosphoDecoy (A)") else: target_mods = ("Phospho (S)", "Phospho (T)", "Phospho (Y)") - ctx = click.Context(lucxor_func) - ctx.invoke( - lucxor_func, - input_spectrum=in_file, input_id=id_file, output=lucxor_out, + lucxor_setup_logging(debug, None, lucxor_out) + tool = PyLuciPHOr2() + tool.run( + input_spectrum=in_file, + input_id=id_file, + output=lucxor_out, fragment_method=fragment_method, fragment_mass_tolerance=fragment_mass_tolerance, fragment_error_units=fragment_mass_unit, @@ -184,16 +194,18 @@ def all( decoy_mass=79.966331, decoy_neutral_losses=("X -H3PO4 -97.97690",), max_charge_state=5, max_peptide_length=40, max_num_perm=16384, - modeling_score_threshold=0.95, scoring_threshold=0.0, + modeling_score_threshold=modeling_score_threshold, scoring_threshold=0.0, min_num_psms_model=50, threads=threads, rt_tolerance=0.01, - debug=debug, log_file=None, disable_split_by_charge=False, + debug=debug, disable_split_by_charge=False, ) + + # Merge results click.echo(f"\n{'='*60}") click.echo(f"[{time.strftime('%H:%M:%S')}] Merging results...") click.echo(f"{'='*60}") - merge_algorithm_results(ascore_out, phosphors_out, lucxor_out, out_file) + merge_algorithm_results(ascore_out, phosphors_out, lucxor_out, out_file, id_file) elapsed = time.time() - start_time click.echo(f"\n{'='*60}") @@ -408,8 +420,24 @@ def merge_algorithm_results(ascore_file, phosphors_file, lucxor_file, output_fil out_df = full_df.reset_index() + # Count how many peptides had their modification sites reassigned by LucXor + relocated_count = 0 + input_psms_df, _, _, _ = load_dataframes(input_idparquet) + input_pf_map = {} + for _, row in input_psms_df.iterrows(): + pep_idx = row.get("peptide_identification_index") + if pep_idx is not None: + input_pf_map[pep_idx] = str(row.get("peptidoform", "")) + for _, row in out_df.iterrows(): + pep_idx = row.get("peptide_identification_index") + out_pf = str(row.get("peptidoform", "")) + orig_pf = input_pf_map.get(pep_idx, "") + if orig_pf and orig_pf != out_pf: + relocated_count += 1 + save_dataframes(output_file, out_df, proteins_df, template_df=lucxor_df, source_idparquet=input_idparquet) click.echo(f"Successfully merged {stats['merged']} peptide identifications") + click.echo(f" Modification sites reassigned: {relocated_count}") click.echo("Each peptide contains scores from all three algorithms") @@ -418,6 +446,7 @@ def run_all_algorithms_from_single_cli( fragment_mass_tolerance, fragment_mass_unit, threads, debug, add_decoys, fragment_method="CID", + modeling_score_threshold=0.95, ): """Run all three algorithms when --compute-all-scores is specified.""" try: @@ -470,7 +499,7 @@ def run_all_algorithms_from_single_cli( decoy_mass=79.966331, decoy_neutral_losses=("X -H3PO4 -97.97690",), max_charge_state=5, max_peptide_length=40, max_num_perm=16384, - modeling_score_threshold=0.95, scoring_threshold=0.0, + modeling_score_threshold=modeling_score_threshold, scoring_threshold=0.0, min_num_psms_model=50, threads=threads, rt_tolerance=0.01, debug=debug, disable_split_by_charge=False, ) diff --git a/onsite/phosphors/cli.py b/onsite/phosphors/cli.py index c493355..4cd5f32 100644 --- a/onsite/phosphors/cli.py +++ b/onsite/phosphors/cli.py @@ -253,11 +253,27 @@ def phosphors( f"Error processing identification: {res.get('reason', 'unknown')}" ) + # Count how many PSMs had their modification sites reassigned + relocated_count = 0 + orig_pf_map = {} + for _, group_df in grouped: + pep_idx = group_df.iloc[0].get("peptide_identification_index") + orig_pf = str(group_df.iloc[0].get("peptidoform", "")) + if orig_pf: + orig_pf_map[pep_idx] = orig_pf + for row in result_rows: + pep_idx = row.get("peptide_identification_index") + orig_pf = orig_pf_map.get(pep_idx, "") + out_pf = str(row.get("peptidoform", "")) + if orig_pf and orig_pf != out_pf: + relocated_count += 1 + # Report elapsed = time.time() - start_time click.echo(f"\nProcessing Complete:") click.echo(f" Total identifications: {stats['total']}") click.echo(f" Successfully processed: {stats['processed']}") + click.echo(f" Modification sites reassigned: {relocated_count}") click.echo(f" Phosphorylated peptides: {stats['phospho']}") click.echo(f" Processing errors: {stats['errors']}") click.echo(f" Time elapsed: {elapsed:.2f} seconds") diff --git a/pyproject.toml b/pyproject.toml index 6df903a..835ffd8 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,6 +1,6 @@ [project] name = "pyonsite" -version = "0.0.3" +version = "0.0.4" description = "onsite: mass spectrometry post-translational localization tool" authors = [ {name = "BigBio Stack"} diff --git a/tests/test_output_validation.py b/tests/test_output_validation.py index faacc09..19ee221 100644 --- a/tests/test_output_validation.py +++ b/tests/test_output_validation.py @@ -57,6 +57,7 @@ def test_output_content_validation(self, idparquet_dir, mzml_file): "--input-id", idparquet_dir, "--output", out, "--target-modifications", "Phospho(S),Phospho(T),Phospho(Y)", + "--modeling-score-threshold", "0.3", ]) if r.exit_code != 0: pytest.skip(f"LucXor required but failed: {r.output[:200]}")