Create Ascore.py

[weizhongchun]

PR Type

Enhancement

---

Description

• Implemented the AScore class for phosphorylation site localization.

• Added methods for computing AScore, handling permutations, and generating theoretical spectra.

• Integrated support for decoy phosphorylation sites and ProForma string generation.

• Included detailed debugging and error handling for edge cases.

---

Changes walkthrough 📝

	Relevant files
Enhancement	Ascore.py Implement AScore algorithm for phosphorylation analysis    Ascore.py Added the AScore class for phosphorylation site localization. Implemented methods for theoretical spectrum generation and scoring. Included support for decoy sites and ProForma string generation. Added debugging and error handling for edge cases. +1009/-0

---

Need help?

Type /help how to ... in the comments thread for any questions about Qodo Merge usage.

Check out the documentation for more information.

Summary by CodeRabbit

• New Features
  • Introduced an advanced algorithm to enhance phosphorylation site localization in peptide analysis.
  • Improves the scoring and ranking of potential modification sites based on spectral data.
  • Includes enhanced error handling to address edge cases during the computation process.

[coderabbitai]

Walkthrough

A new file, Ascore.py, has been added to implement the AScore algorithm for phosphorylation site localization in peptides using MS/MS spectra. The file introduces the AScore class that initializes with default parameters and provides methods to identify phosphorylation sites, generate site permutations, compute probability scores using cumulative binomial calculations, and generate theoretical spectra. The class integrates error handling and debug statements, ultimately returning a modified PeptideHit with best AScore values and associated metadata.

Changes

File

Change Summary

Ascore.py

• Added AScore class implementing the AScore algorithm • Initialization & Utility: __init__, updateMembers_, isPhosphoDecoySite, isPhosphoSite, numberOfPhosphoEvents_, getSites_, getPhosphoSites_, getPhosphoDecoySites_, removePhosphositesFromSequence_ • Core Computation: compute, computeBaseProbability_, computeCumulativeScore_, peptideScore_, compareMZ_ • Permutation & Scoring: computePermutations_, determineHighestScoringPermutations_, calculatePermutationPeptideScores_, rankWeightedPermutationPeptideScores_ • Spectra Analysis: createTheoreticalSpectra_, getSpectrumDifference_, numberOfMatchedIons_, peakPickingPerWindowsInSpectrum_, computeSiteDeterminingIons_ • Output Formatting: generateProFormaString_

Sequence Diagram(s)

sequenceDiagram
    participant C as Client
    participant A as AScore
    C->>A: compute(hit, real_spectrum)
    A->>A: Validate input & identify phosphorylation sites
    A->>A: Count phospho events and generate site permutations
    A->>A: Create theoretical spectra and pick spectral peaks
    A->>A: Calculate and rank permutation scores
    A->>A: Format results into ProForma string
    A-->>C: Return modified PeptideHit with AScore values and metadata

Loading

Poem

In code I hop with nimble feet,
Phospho secrets now complete,
Each method a step in spectral dance,
AScore sings in algorithmic trance,
Bunny cheers and hops with glee—coding ever so sweet!
🐇✨

Tip

⚡💬 Agentic Chat (Pro Plan, General Availability)

• We're introducing multi-step agentic chat in review comments and issue comments, within and outside of PR's. This feature enhances review and issue discussions with the CodeRabbit agentic chat by enabling advanced interactions, including the ability to create pull requests directly from comments and add commits to existing pull requests.

✨ Finishing Touches

• 📝 Generate Docstrings

---

Thanks for using CodeRabbit! It's free for OSS, and your support helps us grow. If you like it, consider giving us a shout-out.

❤️ Share

• X
• Mastodon
• Reddit
• LinkedIn

🪧 Tips

Chat

There are 3 ways to chat with CodeRabbit:

• Review comments: Directly reply to a review comment made by CodeRabbit. Example:
  • I pushed a fix in commit <commit_id>, please review it.
  • Generate unit testing code for this file.
  • Open a follow-up GitHub issue for this discussion.
• Files and specific lines of code (under the "Files changed" tab): Tag @coderabbitai in a new review comment at the desired location with your query. Examples:
  • @coderabbitai generate unit testing code for this file.
  • @coderabbitai modularize this function.
• PR comments: Tag @coderabbitai in a new PR comment to ask questions about the PR branch. For the best results, please provide a very specific query, as very limited context is provided in this mode. Examples:
  • @coderabbitai gather interesting stats about this repository and render them as a table. Additionally, render a pie chart showing the language distribution in the codebase.
  • @coderabbitai read src/utils.ts and generate unit testing code.
  • @coderabbitai read the files in the src/scheduler package and generate a class diagram using mermaid and a README in the markdown format.
  • @coderabbitai help me debug CodeRabbit configuration file.

Note: Be mindful of the bot's finite context window. It's strongly recommended to break down tasks such as reading entire modules into smaller chunks. For a focused discussion, use review comments to chat about specific files and their changes, instead of using the PR comments.

CodeRabbit Commands (Invoked using PR comments)

• @coderabbitai pause to pause the reviews on a PR.
• @coderabbitai resume to resume the paused reviews.
• @coderabbitai review to trigger an incremental review. This is useful when automatic reviews are disabled for the repository.
• @coderabbitai full review to do a full review from scratch and review all the files again.
• @coderabbitai summary to regenerate the summary of the PR.
• @coderabbitai generate docstrings to generate docstrings for this PR.
• @coderabbitai resolve resolve all the CodeRabbit review comments.
• @coderabbitai plan to trigger planning for file edits and PR creation.
• @coderabbitai configuration to show the current CodeRabbit configuration for the repository.
• @coderabbitai help to get help.

Other keywords and placeholders

• Add @coderabbitai ignore anywhere in the PR description to prevent this PR from being reviewed.
• Add @coderabbitai summary to generate the high-level summary at a specific location in the PR description.
• Add @coderabbitai anywhere in the PR title to generate the title automatically.

CodeRabbit Configuration File (.coderabbit.yaml)

• You can programmatically configure CodeRabbit by adding a .coderabbit.yaml file to the root of your repository.
• Please see the configuration documentation for more information.
• If your editor has YAML language server enabled, you can add the path at the top of this file to enable auto-completion and validation: # yaml-language-server: $schema=https://coderabbit.ai/integrations/schema.v2.json

Documentation and Community

• Visit our Documentation for detailed information on how to use CodeRabbit.
• Join our Discord Community to get help, request features, and share feedback.
• Follow us on X/Twitter for updates and announcements.

[qodo-code-review]

PR Reviewer Guide 🔍

Here are some key observations to aid the review process:

⏱️ Estimated effort to review: 4 🔵🔵🔵🔵⚪

🧪 No relevant tests

🔒 No security concerns identified

⚡ Recommended focus areas for review Numerical Stability The cumulative binomial probability calculation may have numerical stability issues for large N values. The fallback implementation when scipy is not available could lead to overflow or underflow errors. try: from scipy.stats import binom return 1.0 - binom.cdf(n-1, N, p) except ImportError: # Fallback implementation: Calculate P(X ≥ n) directly # This is less numerically stable but works without scipy total_prob = 0.0 for k in range(n, N+1): # Calculate binomial coefficient * p^k * (1-p)^(N-k) coef = 1.0 for i in range(1, k+1): coef *= (N - k + i) / i prob = coef * (p**k) * ((1-p)**(N-k)) total_prob += prob return total_prob Error Handling The code has several places where exceptions could be raised without proper handling, particularly in the permutation generation logic where math operations might fail for large inputs. try: from math import comb estimated_permutations = comb(len(sites), n_phosphorylation_events) except (ImportError, OverflowError): # Fallback calculation try: estimated_permutations = 1 for i in range(1, n_phosphorylation_events + 1): estimated_permutations *= (len(sites) - i + 1) / i except OverflowError: estimated_permutations = float('inf') if estimated_permutations > self.max_permutations_: print(f"Debug: Early termination: estimated permutations ({estimated_permutations}) " f"exceeds maximum ({self.max_permutations_})") return permutations Performance Concern The spectrum difference calculation and ion matching algorithms use nested loops that could be inefficient for large spectra, potentially causing performance bottlenecks. """ result.clear(True) # Use True parameter to clear all data i, j = 0, 0 while i < spectrum1.size() and j < spectrum2.size(): mz1 = spectrum1[i].getMZ() mz2 = spectrum2[j].getMZ() val = self.compareMZ_(mz1, mz2) if val == -1: # mz1 < mz2 (outside tolerance) result.push_back(spectrum1[i]) i += 1 elif val == 1: # mz1 > mz2 (outside tolerance) j += 1 else: # within tolerance, check if more ions are within tolerance # Check mz2 until no match j += 1 while j < spectrum2.size() and self.compareMZ_(mz1, spectrum2[j].getMZ()) == 0: j += 1 # Check mz1 until no match i += 1 while i < spectrum1.size() and self.compareMZ_(spectrum1[i].getMZ(), mz2) == 0: i += 1 # Add remaining peaks from spectrum1 while i < spectrum1.size(): result.push_back(spectrum1[i]) i += 1

[qodo-code-review]

PR Code Suggestions ✨

Explore these optional code suggestions:

Category	Suggestion	Impact
Possible issue	Fix permutation selection logic The logic for finding alternative permutations has a critical bug. The continue_search flag is incorrectly set to False inside the inner loop, but then never reset to True if a mismatch is found. This causes the algorithm to potentially select incorrect permutations for comparison, leading to inaccurate AScore calculations. Ascore.py [208-258] def determineHighestScoringPermutations_(self, peptide_site_scores, sites, permutations, ranking): """ Determine the highest scoring permutations for each phosphorylation site. For each phosphorylation site in the highest-scoring permutation, this function: 1. Finds the next best permutation where this site is not phosphorylated 2. Determines the peak depth that maximizes the score difference between these permutations 3. Stores this information for AScore calculation """ sites.clear() # Check if ranking is empty if not ranking: print("Warning: No ranking available for permutation scoring") return # Check if there is only one permutation if len(ranking) == 1: print("Warning: Only one permutation available, cannot compute AScore") return # Take first set of phospho site assignments best_peptide_sites = permutations[ranking[0][1]] # Sites of the assignment that achieved the highest weighted score for i in range(len(best_peptide_sites)): # For each phosphorylated site site_info = {'first': best_peptide_sites[i], 'seq_1': ranking[0][1]} - continue_search = True rank_index = 1 # Start with the next-best scoring peptide # Iterate from best scoring peptide to the first peptide that doesn't contain the current phospho site - while continue_search and rank_index < len(ranking): + while rank_index < len(ranking): current_permutation = permutations[ranking[rank_index][1]] + is_valid_alternative = True # Check if this permutation has each of the phosphorylation sites as the best one, # except for the current site (i) which should be different for j in range(len(best_peptide_sites)): if j == i: # The site we are interested in if best_peptide_sites[j] in current_permutation: - continue_search = True # Permutation is also phosphorylated -> skip + is_valid_alternative = False # Permutation is also phosphorylated -> skip break - else: - continue_search = False else: # The other sites if best_peptide_sites[j] not in current_permutation: - continue_search = True # Phosphorylation state should be the same but isn't + is_valid_alternative = False # Phosphorylation state should be the same but isn't break - else: - continue_search = False - if continue_search: - rank_index += 1 + if is_valid_alternative: + break # Found a valid alternative permutation + rank_index += 1 [To ensure code accuracy, apply this suggestion manually] Suggestion importance[1-10]: 9 __ Why: The suggestion fixes a critical logic bug in the permutation selection algorithm. The original code incorrectly sets and resets the continue_search flag, potentially selecting wrong permutations for comparison, which would lead to inaccurate AScore calculations and incorrect phosphorylation site localization.	High
	Fix numerical overflow risk The fallback implementation for calculating binomial probabilities has a numerical stability issue. When calculating the binomial coefficient, multiplying large numbers can lead to overflow, especially for large values of N. The implementation should use a more stable approach like logarithms to avoid overflow. Ascore.py [296-326] def computeCumulativeScore_(self, N, n, p): """ Compute the cumulative binomial probability P(X ≥ n). Parameters: ----------- N : int Total number of trials (theoretical peaks) n : int Number of successes (matched peaks) p : float Probability of success for a single trial Returns: -------- float Cumulative binomial probability P(X ≥ n) """ assert n <= N, "The number of matched ions (n) can be at most as large as the number of trials (N)." assert 0 <= p <= 1.0, "p must be a probability [0,1]." # Use scipy's implementation for binomial CDF complement if available try: from scipy.stats import binom return 1.0 - binom.cdf(n-1, N, p) except ImportError: - # Fallback implementation: Calculate P(X ≥ n) directly - # This is less numerically stable but works without scipy + # Fallback implementation: Calculate P(X ≥ n) using logarithms for numerical stability + import math total_prob = 0.0 for k in range(n, N+1): - # Calculate binomial coefficient * p^k * (1-p)^(N-k) - coef = 1.0 + # Calculate log of binomial coefficient + log(p^k) + log((1-p)^(N-k)) + log_coef = 0.0 for i in range(1, k+1): - coef *= (N - k + i) / i - prob = coef * (p**k) * ((1-p)**(N-k)) - total_prob += prob + log_coef += math.log(N - k + i) - math.log(i) + log_prob = log_coef + k * math.log(p) + (N-k) * math.log(1-p) + # Add exp(log_prob) to total probability + total_prob += math.exp(log_prob) return total_prob [To ensure code accuracy, apply this suggestion manually] Suggestion importance[1-10]: 8 __ Why: The improved implementation uses logarithms to calculate binomial probabilities, which significantly reduces the risk of numerical overflow when dealing with large values of N. This is a critical fix for the fallback implementation that would otherwise fail for larger peptides.	Medium
	Fix numerical instability risk The current code handles the case where there's no second permutation by setting peptide2_score to a very small value, but this can lead to division by zero or numerical instability when calculating the AScore. Instead, use a more robust approach that avoids potential numerical issues. Ascore.py [196-203] # Ensure there is a second permutation for comparison if rev + 1 < len(sorted_ranking): seq2 = th_spectra[sorted_ranking[rev+1][1]].getName() peptide2_score = sorted_ranking[rev+1][0] else: - # If there is no second permutation, use a very low score + # If there is no second permutation, use the same sequence but ensure + # the score difference will result in a meaningful AScore seq2 = seq1 - peptide2_score = peptide1_score * 0.0000001 + peptide2_score = max(peptide1_score * 0.5, peptide1_score - 10.0) Apply this suggestion Suggestion importance[1-10]: 7 __ Why: The suggestion improves numerical stability by replacing the arbitrary small multiplier (0.0000001) with a more robust approach that ensures a meaningful score difference. This prevents potential division by zero or extreme values when calculating AScores.	Medium
More

[coderabbitai]

Actionable comments posted: 0

🧹 Nitpick comments (6)

Ascore.py (6)

1-1: Avoid wildcard imports for maintainability.

Wildcard imports (from pyopenms import *) make it unclear which symbols are actually used and can lead to name collisions or confusion around undefined names. Consider importing only the specific classes and functions needed (e.g., from pyopenms import AASequence, MSSpectrum, TheoreticalSpectrumGenerator, PeptideHit).

-from pyopenms import *
+from pyopenms import AASequence, MSSpectrum, TheoreticalSpectrumGenerator, PeptideHit

🧰 Tools

🪛 Ruff (0.8.2)

1-1: from pyopenms import * used; unable to detect undefined names

(F403)

---

3-3: Remove unused import numpy.

The import statement for numpy is never referenced in this file. Consider removing it to keep the code clean.

-import numpy as np

🧰 Tools

🪛 Ruff (0.8.2)

3-3: numpy imported but unused

Remove unused import: numpy

(F401)

---

4-4: Remove unused import defaultdict.

collections.defaultdict is not referenced anywhere in this file. Consider removing it.

-from collections import defaultdict

🧰 Tools

🪛 Ruff (0.8.2)

4-4: collections.defaultdict imported but unused

Remove unused import: collections.defaultdict

(F401)

---

78-79: Use logging instead of print statements for debug and warnings.

Multiple debugging and warning prints are scattered throughout the code. Using the logging module allows for flexible log levels, better formatting, and easier maintenance.

-print(f"Debug: Found {regular_phospho_count} regular phosphorylation events ...")
-...
-print("Warning: PhosphoDecoy sites found in sequence ...")
-...
-print(f"Debug: Calculation aborted: peptide too long: {seq_without_phospho.toString()}")
-...
-print("Debug: Calculation aborted: number of permutations exceeded ...")
+import logging
+logger = logging.getLogger(__name__)
+...
+logger.debug(f"Found {regular_phospho_count} regular phosphorylation events ...")
+...
+logger.warning("PhosphoDecoy sites found in sequence ...")
+...
+logger.debug(f"Calculation aborted: peptide too long: {seq_without_phospho.toString()}")
+...
+logger.debug("Calculation aborted: number of permutations exceeded ...")

Also applies to: 93-99, 103-104, 146-147

---

135-137: Rename unused loop variable site.

In this loop, site is never referenced within the body. Consider renaming it to _site or removing it entirely to avoid confusion.

-for i, site in enumerate(sites):
+for i, _site in enumerate(sites):
    phospho.setMetaValue(f"AScore_{i+1}", self.unambiguous_score_)

🧰 Tools

🪛 Ruff (0.8.2)

135-135: Loop control variable site not used within loop body

(B007)

---

181-181: Rename unused loop variable score.

The loop extracts (score, idx) but never uses score. Consider renaming it to _score for clarity.

-for i, (score, idx) in enumerate(best_permutations):
+for i, (_score, idx) in enumerate(best_permutations):

🧰 Tools

🪛 Ruff (0.8.2)

181-181: Loop control variable score not used within loop body

Rename unused score to _score

(B007)

📜 Review details

Configuration used: CodeRabbit UI
Review profile: CHILL
Plan: Pro

📥 Commits

Reviewing files that changed from the base of the PR and between 7f25c4c and f1fd385.

📒 Files selected for processing (1)

• Ascore.py (1 hunks)

🧰 Additional context used

🪛 Ruff (0.8.2)

Ascore.py

1-1: from pyopenms import * used; unable to detect undefined names

(F403)

---

3-3: numpy imported but unused

Remove unused import: numpy

(F401)

---

4-4: collections.defaultdict imported but unused

Remove unused import: collections.defaultdict

(F401)

---

31-31: TheoreticalSpectrumGenerator may be undefined, or defined from star imports

(F405)

---

71-71: PeptideHit may be undefined, or defined from star imports

(F405)

---

135-135: Loop control variable site not used within loop body

(B007)

---

181-181: Loop control variable score not used within loop body

Rename unused score to _score

(B007)

---

191-191: AASequence may be undefined, or defined from star imports

(F405)

---

409-409: MSSpectrum may be undefined, or defined from star imports

(F405)

---

409-409: MSSpectrum may be undefined, or defined from star imports

(F405)

---

415-415: MSSpectrum may be undefined, or defined from star imports

(F405)

---

419-419: MSSpectrum may be undefined, or defined from star imports

(F405)

---

482-482: MSSpectrum may be undefined, or defined from star imports

(F405)

---

485-485: MSSpectrum may be undefined, or defined from star imports

(F405)

---

775-775: AASequence may be undefined, or defined from star imports

(F405)

---

802-802: MSSpectrum may be undefined, or defined from star imports

(F405)

---

808-808: AASequence may be undefined, or defined from star imports

(F405)

---

862-862: MSSpectrum may be undefined, or defined from star imports

(F405)

---

868-868: MSSpectrum may be undefined, or defined from star imports

(F405)