<jupyter_start><jupyter_text>intensity ion backbone scoreThe ion backbone score, but also taking into account the relative intensity that is being made up by identified peaks<jupyter_code>import os
import sys
module_path = os.path.abspath(os.path.join('..'))
if module_path not in sys.path:
    sys.path.append(module_path)
module_path = os.path.abspath(os.path.join('../..'))
if module_path not in sys.path:
    sys.path.append(module_path)
    
from src.objects import Spectrum
from src.spectra.gen_spectra import gen_spectrum
from src.utils import ppm_to_da
sum([1 if jc == 0 else 0 for jc in jcount])
def intensity_ion_backbone_score(observed: Spectrum, reference: str, ion: str, ppm_tolerance: int) -> float:
    '''
    Scoring algorithm that factors in how much of the backbone is identified along with the abundance each peak 
    contributes to the total score. The formula is:

        count(# backbone cleavages found) + 1 for every additional hit of an already identifed cleavage * total percentage
        of intensity covered by identfied peaks

    Example:
        reference:   ABCDE, 4 junctions to describe
        observed:    ions: b1+, y1++, y2+, b4+, with relative intensities of (100, 200, 100, 50) of a total 600
        
        ions for A, E, DE, ABCD found. Coverage = A**DE with DE described by both E and ABCD
        
        Score = %(3/4) + 1 = 75 + 1 = 76 *(100 + 200 + 100 + 50)/600 = 57 

    Inputs:
        observed:       (Spectrum) spectrum being scored against
        reference:      (str) reference amino acid sequence being scored against the spectrum
        ion:            (str) the ion type to focus on. Options are 'b' or 'y'
        ppm_tolerance:  (int) tolerance to allow in ppm for each peak
    Outputs:
        (float) score according to the formula 
    '''
    # check to see if observed is nothing
    if len(observed.spectrum) == 0:
        return 0

    # keep track of the junction (bond) sites found
    jcount = [0 for _ in range(len(reference)-1)]

    # keep track of the abundances that contribute to our score
    ided_abundances = 0
    
    for charge in [1, 2]:

        # take off the trailing or leading amino acid from the reference according to ion type
        singled_seq = reference[:-1] if ion == 'b' else reference[1:]

        # get the m/z peaks
        peaks = gen_spectrum(singled_seq, charge=charge, ion=ion)['spectrum']
        
        # go through each peak and try and match it to an observed one
        for i in range(len(peaks)):

            # take tolerance into account 
            da_tol = ppm_to_da(peaks[i], ppm_tolerance)

            # get hits
            peak_hits = list(map(lambda idx_x: idx_x[0] if peaks[i] - da_tol <= idx_x[1] <= peaks[i] + da_tol else None, enumerate(observed.spectrum)))

            # remove None from peak hits
            peak_hits = [idx for idx in peak_hits if idx is not None]

            # if the len > 1, increment
            if len(peak_hits) > 0:
                jcount[i] += 1
                ided_abundances += sum([observed.abundance[idx] for idx in peak_hits])

    # if any entry has at least 1, that bond has been identified at least once
    jcoverage = sum([1 if jc > 0 else 0 for jc in jcount])

    # if an entry has more than 1, we give it extra points
    extrapoints = sum([jc - 1 if jc > 1 else 0 for jc in jcount]) // 3
    
    return (jcoverage + extrapoints) * (ided_abundances / observed.total_intensity)
def intensity_ion_backbone_score2(observed: Spectrum, reference: str, ion: str, ppm_tolerance: int) -> float:
    '''
    Scoring algorithm that factors in how much of the backbone is identified along with the abundance each peak 
    contributes to the total score. The formula is:

        count(# backbone cleavages found) + 1 for every additional hit of an already identifed cleavage * total percentage
        of intensity covered by identfied peaks

    Example:
        reference:   ABCDE, 4 junctions to describe
        observed:    ions: b1+, y1++, y2+, b4+, with relative intensities of (100, 200, 100, 50) of a total 600
        
        ions for A, E, DE, ABCD found. Coverage = A**DE with DE described by both E and ABCD
        
        Score = %(3/4) + 1 = 75 + 1 = 76 *(100 + 200 + 100 + 50)/600 = 57 

    Inputs:
        observed:       (Spectrum) spectrum being scored against
        reference:      (str) reference amino acid sequence being scored against the spectrum
        ion:            (str) the ion type to focus on. Options are 'b' or 'y'
        ppm_tolerance:  (int) tolerance to allow in ppm for each peak
    Outputs:
        (float) score according to the formula 
    '''
    # check to see if observed is nothing
    if len(observed.spectrum) == 0:
        return 0

    # keep track of the junction (bond) sites found
    jcount = [0 for _ in range(len(reference)-1)]

    # keep track of the abundances that contribute to our score
    ided_abundances = 0
    
    for charge in [1, 2]:

        # take off the trailing or leading amino acid from the reference according to ion type
        singled_seq = reference[:-1] if ion == 'b' else reference[1:]

        # get the m/z peaks
        peaks = gen_spectrum(singled_seq, charge=charge, ion=ion)['spectrum']
        
        # go through each peak and try and match it to an observed one
        for i in range(len(peaks)):

            # take tolerance into account 
            da_tol = ppm_to_da(peaks[i], ppm_tolerance)

            # get hits
            peak_hits = list(map(lambda idx_x: idx_x[0] if peaks[i] - da_tol <= idx_x[1] <= peaks[i] + da_tol else None, enumerate(observed.spectrum)))

            # remove None from peak hits
            peak_hits = [idx for idx in peak_hits if idx is not None]

            # if the len > 1, increment
            if len(peak_hits) > 0:
                jcount[i] += 1
                ided_abundances += sum([observed.abundance[idx] for idx in peak_hits])

    # if any entry has at least 1, that bond has been identified at least once
    jcoverage = sum([1 if jc > 0 else 0 for jc in jcount])
    
    # subtract points for every missing bond 
    missing_bonds = sum([1 if jc == 0 else 0 for jc in jcount])
    print(f'Missing {missing_bonds}/{len(jcount)} bonds')
    jcoverage -= missing_bonds

    # if an entry has more than 1, we give it extra points
    extrapoints = sum([jc - 1 if jc > 1 else 0 for jc in jcount]) // 3
    
    return (jcoverage + extrapoints) * (ided_abundances / observed.total_intensity)
seq1 = 'MALWARM'
seq2 = 'MALW'
seq3 = 'WARM'
seq4 = 'QQQPS'

observed_spec = gen_spectrum(seq1)['spectrum']
abundances = [1 for _ in range(len(observed_spec))]
observed = Spectrum(observed_spec, abundances, sum(abundances))

print(intensity_ion_backbone_score(observed, seq2, 'b', 20))
print(intensity_ion_backbone_score2(observed, seq2, 'b', 20))
print(intensity_ion_backbone_score(observed, seq3, 'y', 20))
print(intensity_ion_backbone_score2(observed, seq3, 'y', 20))
print(intensity_ion_backbone_score(observed, seq2, 'y', 20))
print(intensity_ion_backbone_score2(observed, seq2, 'y', 20))
intensity_ion_backbone_score(observed, seq3, 'b', 20)
intensity_ion_backbone_score(observed, seq4, 'b', 20)
intensity_ion_backbone_score(observed, seq4, 'y', 20)
from src.file_io.mzML import read

mzmlfile = '/Users/zacharymcgrath/Desktop/nod2 data/single/singleRealSpectrum.mzml'

spec = None
for e in read(mzmlfile):
    spec = e
    break
    
print(spec)
intensity_ion_backbone_score(spec, 'NFEANTTIGRIRFH', 'y', 20) * len('NFEANTTIGRIRFH') // 2
intensity_ion_backbone_score(spec, 'LPS', 'y', 20) * len('LPS') // 2
intensity_ion_backbone_score(spec, 'TDS', 'b', 20)
intensity_ion_backbone_score(spec, 'DEA', 'b', 20)
gen_spectrum('EAPNFEANTTIGRIRFHD', charge=2)
gen_spectrum('DEAPNFEANTTIGRIRFH', charge=2)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.09266802443991853, 'fraction_numerical': 0.024185336048879836, 'mean_word_length': 2.8973214285714284, 'pattern_counts': {'":': 0, '<': 9, '<?xml version=': 0, '>': 15, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '1', 'dir': 'jupyter-structured-clean-dedup', 'id': '574331', 'license': 'MIT', 'n_tokens_mistral': 2698, 'n_tokens_neox': 2398, 'n_words': 898, 'path': 'sandbox/jupyter notebooks/intensity based ion backbone score.ipynb', 'repo_name': 'jamessdixon/hypedsearch', 'size': '7856'}