<jupyter_start><jupyter_text>ExtensionsThe purpose of this notebook is to test everything related to generating extensions among hits<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)

from testing_framework import testing_utils
import database
from preprocessing import merge_search
from identification import create_hits
import utils

import operator

ppm_tolerance = 20
max_peptide_length = 20
precursor_tolerance = 10

import matplotlib.pyplot as plt
datasets = testing_utils.define_data()

dataset = datasets[0]

input_spectra_path = dataset[0]
input_spectra, boundaries = testing_utils.preprocess_input_spectra(input_spectra_path, ppm_tolerance)

correct_sequences = testing_utils.generate_truth_set(datasets[0])

path = dataset[2]
db = database.build(path)
write_path = os.path.abspath(os.path.join(module_path, 'intermediate_files'))
matched_masses_b, matched_masses_y, kmer_set = merge_search.modified_match_masses(boundaries, db, max_peptide_length, True, write_path)
print('Finished matching masses')<jupyter_output>On protein 279/279 [100%]
Sorting the set of protein masses...
Sorting the set of protein masses done
Performing Merge
Done
Finished matching masses<jupyter_text>Getting initial hits<jupyter_code>from gen_spectra import gen_spectrum
spectrum_num = 0

correct_sequence = correct_sequences[spectrum_num]
print(correct_sequence)

input_spectrum = input_spectra[spectrum_num]
unique_b,unique_y = testing_utils.get_unique_matched_masses(boundaries, matched_masses_b, matched_masses_y)
location = os.path.join(os.path.abspath(os.path.join('../..')), 'intermediate_files/')
b_hits, y_hits = create_hits(spectrum_num, input_spectrum, matched_masses_b, matched_masses_y, False, location)
correct_hits = testing_utils.append_correct_hits(correct_sequence, input_spectrum, ppm_tolerance)
ion = 'b'
clusters = testing_utils.create_clusters(ion, b_hits, y_hits)
b_sorted_clusters = testing_utils.Bayes_clusters(ion, clusters, location, kmer_set, unique_b)
ion = 'y'
clusters = testing_utils.create_clusters(ion, b_hits, y_hits)
y_sorted_clusters = testing_utils.Bayes_clusters(ion, clusters, location, kmer_set, unique_y)<jupyter_output><empty_output><jupyter_text>Printing hits<jupyter_code>b_sorted_clusters = sorted(b_sorted_clusters, key=operator.attrgetter('prob', 'score', 'pid'))
# b_sorted_clusters = sorted(b_sorted_clusters, key=operator.attrgetter('score', 'prob', 'pid'), reverse = True)
testing_utils.write_b_sorted_cluster(b_sorted_clusters)
for i in range(0, 50):
    x = b_sorted_clusters[i]
    post_prob = x.prob
    score = x.score
    seq = x.seq
    indices = x.indices
    print(post_prob, score, seq)
y_sorted_clusters = sorted(y_sorted_clusters, key=operator.attrgetter('prob', 'score', 'pid'))
# y_sorted_clusters = sorted(y_sorted_clusters, key=operator.attrgetter('score', 'prob', 'pid'), reverse = True)
testing_utils.write_y_sorted_cluster(y_sorted_clusters)
for i in range(0, 50):
    x = y_sorted_clusters[i]
    post_prob = x.prob
    score = x.score
    seq = x.seq
    indices = x.indices
    print(post_prob, score, seq)<jupyter_output>4.16818804282123e-18 5 LDSFSEI
1.1616080807536527e-16 4 SSNWVGKGFFAVYEAIC
2.512036051166044e-16 4 FVDLTMPYSV
1.05226703890205e-15 4 LESYGLE
1.12385997955871e-15 4 EIPHSELD
1.12385997955871e-15 4 EISSIDEF
3.347163330282923e-15 4 IEVLETDPH
1.2869227751022117e-14 4 HPDSEL
1.2869227751022117e-14 4 CGLYEL
5.663227040807141e-14 3 LSNPTGLQESISDVTTCL
5.7487143274934e-14 3 TAEIASLDSENIDEI
6.37719274049424e-14 3 AAGCKVEAFAVQGEEL
6.601894630036116e-14 3 LLMAASIYFHDQNP
6.601894630036116e-14 3 LEDSDLKKSDIDEI
7.172157438730102e-14 3 LSIHQLVENTDETYCI
7.258536525544357e-14 3 DLNINMTSPMGTKSI
7.258536525544357e-14 3 WLKGQGVYLGMPGCL
7.258536525544357e-14 3 DTNAPAHQLIQTESP
7.258536525544357e-14 3 AKDLDTVASDMMVLL
7.660124169797852e-14 3 YFEEYGKIDTIEI
8.163489183127015e-14 3 YNYVWANCFEITLEL
8.163489183127015e-14 3 QQQPPKQQQQQQQQQ
8.163489183127015e-14 3 QQPPKQQQQQQQQQQ
8.336059504994219e-14 3 EIEGEIKRDFMAAL
8.336059504994219e-14 3 LQCYSEAIKLDPQN
9.91871119161064e-14 3 TSQDARFYAL
9.91871119161064e-14 3 DFTFV[...]<jupyter_text>To Filter data by parent protI want to be able to only view the b or y hits from a certain protein<jupyter_code># b_hits
target_pid = 274
b_target_clusters = []
for cluster in b_sorted_clusters:
    if cluster.pid == target_pid:
        b_target_clusters.append(cluster)

for cluster in b_target_clusters:
    assessment, _ = testing_utils.is_good_hit(cluster.seq, 'b', correct_sequence)
    non_indices = str(cluster.score) + '\t' + str(cluster.prob) + '\t' + str(cluster.pid) + '\t' + cluster.seq + '\t' + str(assessment)
    print(non_indices + '\t'+ '\t'.join([str(o) for o in cluster.indices]))

# y_hits
y_target_clusters = []
for cluster in y_sorted_clusters:
    if cluster.pid == target_pid:
        y_target_clusters.append(cluster)

for cluster in y_target_clusters:
    assessment, _ = testing_utils.is_good_hit(cluster.seq, 'y', correct_sequence)
    non_indices = str(cluster.score) + '\t' + str(cluster.prob) + '\t' + str(cluster.pid) + '\t' + cluster.seq + '\t' + str(assessment)
    print(non_indices + '\t'+ '\t'.join([str(o) for o in cluster.indices]))<jupyter_output>3	7.553986085145072e-13	274	DPQVEQLEL	True	(60, 68, 'DPQVEQLEL', 535.7725269681209)	(67, 68, 'EL', 261.1429443359375)	(68, 68, 'L', 132.10121154785156)
1	1.6930194586815518e-08	274	EALYLVC	False	(37, 43, 'EALYLVC', 810.3919677734375)
1	2.029622598916721e-08	274	LE	False	(66, 67, 'LE', 261.1429443359375)
1	2.029622598916721e-08	274	LE	False	(80, 81, 'LE', 261.1429443359375)
1	2.029622598916721e-08	274	LE	False	(103, 104, 'LE', 261.1429443359375)
1	3.15816271525945e-08	274	LC	False	(30, 31, 'LC', 235.1075439453125)
1	3.15816271525945e-08	274	IC	False	(97, 98, 'IC', 235.1075439453125)
1	6.105773986227816e-07	274	L	True	(3, 3, 'L', 132.10121154785156)
1	6.105773986227816e-07	274	L	True	(4, 4, 'L', 132.10121154785156)
1	6.105773986227816e-07	274	L	True	(8, 8, 'L', 132.10121154785156)
1	6.105773986227816e-07	274	L	True	(10, 10, 'L', 132.10121154785156)
1	6.105773986227816e-07	274	L	True	(11, 11, 'L', 132.10121154785156)
1	6.105773986227816e-07	274	L	True	(13, 13, 'L', 132.10121154785156)
1	6[...]<jupyter_text>Finding optimal "hybrid" combos* Hybrid is in quotation marks because all outputs will be a hybrid and then we can check if it is a non-hybrid<jupyter_code>def get_top_X(b_clusters, y_clusters, top_num):
    filtered_b = []
    filtered_y = []
    b_len = top_num if len(b_clusters) >= top_num else len(b_clusters)
    y_len = top_num if len(y_clusters) >= top_num else len(y_clusters)
    for x in range(0,b_len):
        filtered_b.append(b_clusters[x])
    for x in range(0,y_len):
        filtered_y.append(y_clusters[x])
    return filtered_b, filtered_y

def combine(b_cluster, y_cluster):
    b_start, b_end, y_start, y_end = b_cluster.start, b_cluster.end, y_cluster.start, y_cluster.end
    if b_cluster.pid == y_cluster.pid:
        score_add = 2
        hybrid = False
        if  (b_end <= y_end) and (b_start <= y_start) and (b_end >= y_start): #overlap
            overlap = True
            score_add = 2
            seq = b_cluster.seq
            rem_chars = y_start - b_end
            while (rem_chars >= 0):
                seq = seq + y_cluster.seq[len(y_cluster.seq)-1 - rem_chars]
                rem_chars = rem_chars - 1
        else:                                                                #no overlap
            hybrid = False
            overlap = False
            score_add = 2
            seq = b_cluster.seq + '-' + y_cluster.seq
    else:                                                                    #hybrid
        hybrid = True
        overlap = False
        score_add = 0
        seq = b_cluster.seq + '-' + y_cluster.seq
    return seq, score_add, hybrid, overlap

def filter_by_validity(b_cluster, y_cluster):
    valid = True
    for b in b_cluster.indices:
        for y in y_cluster.indices:
            if b[3] == y[3]:
                valid = False
    return valid

def filter_by_precursor(seq, pc, overlap, obs_prec, precursor_tolerance):
    new_seq = seq.replace("-", "") if overlap == False else seq
    tol = ppm_to_da(obs_prec, precursor_tolerance)
    if get_precursor(new_seq, charge=pc) > obs_prec + tol:
        return False
    else:
        return True

def filter_by_dist(b, y, x):
    if y.start - b.end > x:
        return False
    else:
        return True

def merge_clusters(b_clusters, y_clusters, target_precursor, precursor_tolerance):
    # filtered_b, filtered_y = get_top_X(b_sorted_clusters, y_sorted_clusters, 50)
    #Start with printing overlapping. Then will incorportate boundary overlaps between last of b and first of y
    interesting_combos = []
    for b_cluster in b_sorted_clusters:
        for y_cluster in y_sorted_clusters:
            if b_cluster.start <= y_cluster.end:
                # b_indices = parse_indices(b_cluster.indices)
                # y_indices = parse_indices(y_cluster.indices)
                if filter_by_validity(b_cluster, y_cluster):
                    comb_seq, score_add, hybrid, overlap = combine(b_cluster, y_cluster)
                    if filter_by_precursor(comb_seq, 2, overlap, target_precursor, precursor_tolerance):
                        if filter_by_dist(b_cluster, y_cluster, 10):
                            tup = (comb_seq, b_cluster.score + y_cluster.score + score_add, overlap)
                            interesting_combos.append(tup)

    interesting_combos.sort(key=lambda a: a[1], reverse=True)
    return interesting_combos
print(len(b_sorted_clusters), len(y_sorted_clusters))
def get_unique_sorted_clusters(b_sorted_clusters, y_sorted_clusters):
    b_sorted_set = set()
    for b in b_sorted_clusters:
        seq = b.seq
        b_sorted_set.add(seq)
    y_sorted_set = set()
    for y in y_sorted_clusters:
        seq = y.seq
        y_sorted_set.add(seq)
    return b_sorted_set, y_sorted_set

b_sorted_set, y_sorted_set = get_unique_sorted_clusters(b_sorted_clusters, y_sorted_clusters)
merged_seqs = merge_clusters(b_sorted_clusters, y_sorted_clusters, input_spectrum.precursor_mass, precursor_tolerance)
[print(merged_seqs[x]) for x in range(0,50)]
print(len(b_sorted_set), len(y_sorted_set))
import sys
import collections
import operator

def min_info(cluster):
    return (cluster.pid, cluster.start, cluster.end, cluster.score, cluster.seq)

def bsearch(key, Y):
        lo = -1
        hi = len(Y)
        mid = -1
        while (hi - lo > 1):
            mid = int((hi+lo) / 2)
            if Y[mid].start < key:
                lo = mid
            else:
                hi = mid
        return hi
def Ryan_merge(b_sorted_clusters, y_sorted_clusters):
    Cluster = collections.namedtuple('Cluster', 'score pid seq mass start end ion hits')
    merge_seqs = []

    B = {}
    for c in b_sorted_clusters:
        if c.pid not in B:
            B[c.pid] = []
        B[c.pid].append(c)

    Y = {}
    for c in y_sorted_clusters:
        if c.pid not in Y:
            Y[c.pid] = []
        Y[c.pid].append(c)

    for pid in B:
        if pid not in Y:
            continue

        sorted_B = sorted(B[pid], key=operator.attrgetter('pid', 'start', 'end'))
        sorted_Y = sorted(Y[pid], key=operator.attrgetter('pid', 'start', 'end'))

        for b in sorted_B:
            y_i = bsearch(b.start, sorted_Y)

            if y_i >= len(sorted_Y): break

            y = sorted_Y[y_i]

            while y_i < len(sorted_Y) and y.start - b.end < 10:
                y = sorted_Y[y_i]
                merge_seqs.append((b.score + y.score, b.end - y.start, y.end-b.start,min_info(b), min_info(y)))
                y_i += 1
    return merge_seqs
m = Ryan_merge(b_sorted_clusters, y_sorted_clusters)
m.sort(key = lambda x: x[0], reverse = True) 
[print(m[x]) for x in range(0,50)]
target_combos = merge_clusters(b_target_clusters, y_target_clusters, input_spectrum.precursor_mass, precursor_tolerance)
print(target_combos)
print(len(target_combos))
def modified_filter_by_precursor(b, y, pc, obs_prec, precursor_tolerance):
    new_seq = b + y
    tol = ppm_to_da(obs_prec, precursor_tolerance)
    if get_precursor(new_seq, charge=pc) > obs_prec + tol:
        return False
    else:
        return True

interesting_combos = []
for b in b_sorted_set:
    for y in y_sorted_set:
        if modified_filter_by_precursor(b, y, 2, input_spectrum.precursor_mass, precursor_tolerance):
            interesting_combos.append(b + '-' + y)
print(len(interesting_combos))
print(merge_clusters(b_sorted_clusters, y_sorted_clusters, input_spectrum.precursor_mass, precursor_tolerance))<jupyter_output><empty_output><jupyter_text>How big does mass have to be to get matches from both b and y ions<jupyter_code>from utils import ppm_to_da
ppm_tolerance = 20

m = 1000

print(ppm_to_da(m, ppm_tolerance))<jupyter_output>0.02<jupyter_text>OccuranceQ: How does occurance scale with the length of the peptide<jupyter_code>mz_array = []
occurance_array = []
for key in matched_masses_b:
    for tup in matched_masses_b[key]:
        mz = tup[]<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.10026445032111825, 'fraction_numerical': 0.0881752927842841, 'mean_word_length': 3.6425815503332166, 'pattern_counts': {'":': 0, '<': 31, '<?xml version=': 0, '>': 34, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 2, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '87473', 'license': 'MIT', 'n_tokens_mistral': 5740, 'n_tokens_neox': 4940, 'n_words': 1070, 'path': 'src/testing_framework/jupyter notebooks/Extensions.ipynb', 'repo_name': 'ryanlayerlab/hypedsearch', 'size': '13235'}