From d1e80b06db9ed2593cdd402e85ec718fae2469ba Mon Sep 17 00:00:00 2001
From: pgarrett <pgarrett@scripps.edu>
Date: Mon, 27 Nov 2023 10:27:06 -0800
Subject: [PATCH] fixed accession issue, where sequence wouldn't reset
 properly, and throw an error.

---
 app.py       | 349 +++++++++++++++++++++++++++++++++------------------
 constants.py |  13 +-
 util.py      | 141 ++++++++++++++++++++-
 3 files changed, 379 insertions(+), 124 deletions(-)

diff --git a/app.py b/app.py
index d2246ce..60cf42d 100644
--- a/app.py
+++ b/app.py
@@ -1,34 +1,54 @@
-import pickle
+import uuid
 from collections import Counter
+import random
 
 import pandas as pd
 import streamlit as st
-import numpy as np
-import requests
 from peptacular.sequence import strip_modifications, calculate_sequence_length
 from peptacular.digest import identify_cleavage_sites
 from peptacular.mass import valid_mass_sequence
 from peptacular.constants import AMINO_ACIDS
 from peptacular.spans import calculate_span_coverage
 import matplotlib as mpl
+import regex as reg
 
 from constants import *
 from wiki import *
-from util import make_clickable, generate_peptide_df, bin_aa_counts, coverage_string, create_colorbar
-
-
-def fetch_sequence_from_uniprot(accession_number):
-    url = f"https://www.uniprot.org/uniprot/{accession_number}.fasta"
-    response = requests.get(url)
-    if response.status_code != 200:
-        st.error(f"Error fetching sequence from UniProt: {response.status_code}")
-        st.stop()
-        return None
-    return ''.join(response.text.split('\n')[1:])  # Remove the header line
-
+from util import make_clickable, generate_peptide_df, coverage_string, create_colorbar, generate_app_url, \
+    fetch_sequence_from_uniprot
 
 st.set_page_config(page_title="proteincleaver", page_icon=":knife:", layout="wide")
 
+# Parse query parameters
+params = st.experimental_get_query_params()
+query_peptide_sequence = params.get('protein_sequence', [DEFAULT_PROTEIN_SEQUENCE])[0]
+query_proteases = params.get('proteases', [';'.join(DEFAULT_PROTEASES)])[0].split(',')
+query_custom_regex = params.get('custom_regex', [''])[0]
+query_missed_cleavages = int(params.get('missed_cleavages', [DEFAULT_MISSED_CLEAVAGES])[0])
+query_mass_type = params.get('mass_type', [DEFAULT_MASS_TYPE])[0]
+query_min_peptide_len = int(params.get('min_peptide_len', [DEFAULT_MIN_PEPTIDE_LEN])[0])
+query_max_peptide_len = int(params.get('max_peptide_len', [DEFAULT_MAX_PEPTIDE_LEN])[0])
+query_min_mass = float(params.get('min_mass', [DEFAULT_MIN_PEPTIDE_MASS])[0])
+query_max_mass = float(params.get('max_mass', [DEFAULT_MAX_PEPTIDE_MASS])[0])
+query_semi_enzymatic = params.get('semi_enzymatic', ['False'])[0].lower() == 'true'
+query_infer_charge = params.get('infer_charge', ['True'])[0].lower() == 'true'
+query_min_charge = int(params.get('min_charge', [DEFAULT_MIN_CHARGE])[0])
+query_max_charge = int(params.get('max_charge', [DEFAULT_MAX_CHARGE])[0])
+query_min_mz = float(params.get('min_mz', [DEFAULT_MIN_MZ])[0])
+query_max_mz = float(params.get('max_mz', [DEFAULT_MAX_MZ])[0])
+query_remove_non_proteotypic = params.get('remove_non_proteotypic', ['False'])[0].lower() == 'true'
+query_n_term_static_mod = float(params.get('n_term_static_mod', [0.0])[0])
+query_c_term_static_mod = float(params.get('c_term_static_mod', [0.0])[0])
+query_num_static_mods = int(params.get('num_static_mods', [DEFAULT_STATIC_MODS])[0])
+query_n_term_var_mod = float(params.get('n_term_var_mod', [0.0])[0])
+query_c_term_var_mod = float(params.get('c_term_var_mod', [0.0])[0])
+query_max_var_mods = int(params.get('max_var_mods', [DEFAULT_MAX_VAR_MODS])[0])
+query_num_variable_mods = int(params.get('num_variable_mods', [DEFAULT_VAR_MODS])[0])
+query_static_mods_str = params.get('static_mods', ['C:57.02146'])[0]
+query_static_mods = [(s.split(':')[0], float(s.split(':')[1])) for s in query_static_mods_str.split(';') if s]
+query_variable_mods_str = params.get('variable_mods', [''])[0]
+query_variable_mods = [(s.split(':')[0], float(s.split(':')[1])) for s in query_variable_mods_str.split(';') if s]
+
 # CSS to inject contained in a string
 hide_table_row_index = """
             <style>
@@ -39,26 +59,49 @@ def fetch_sequence_from_uniprot(accession_number):
 # Inject CSS with Markdown
 st.markdown(hide_table_row_index, unsafe_allow_html=True)
 
+if 'protein_sequence_key' not in st.session_state:
+    st.session_state['protein_sequence_key'] = 'protein_sequence'
+
 with st.sidebar:
     st.title('Protein Cleaver 🔪')
 
-    protein_id = st.text_input(label='Protein accession number / identifier',
+    c1, c2 = st.columns([3, 1])
+    protein_id = c1.text_input(label='Protein accession number / identifier',
                                value='',
                                help='A protein accession number / identifier')
 
-    raw_sequence = None
+    reset = c2.button('Reset Sequence')
+
+    if reset:
+        st.session_state['protein_sequence_key'] = 'protein_sequence' + str(uuid.uuid4())
+
+    raw_sequence = query_peptide_sequence
+    fetched_sequence = None
     if protein_id:
-        fetched_protein_sequence = fetch_sequence_from_uniprot(protein_id)
-        if fetched_protein_sequence is not None:
-            raw_sequence = fetched_protein_sequence
+        response = fetch_sequence_from_uniprot(protein_id)
+        if response.status_code != 200:
+            st.error(f"Error fetching sequence from UniProt: {response.status_code}")
+            st.stop()
+        fetched_sequence = ''.join(response.text.split('\n')[1:])
+        raw_sequence = fetched_sequence
 
     raw_sequence = st.text_area(label="Protein sequence",
-                                value=raw_sequence if raw_sequence else DEFAULT_PROTEIN_SEQUENCE,
+                                value=raw_sequence,
                                 help='An amino acid sequence to digest',
                                 max_chars=MAX_PROTEIN_INPUT_LENGTH,
-                                height=200)
+                                height=200,
+                                key=st.session_state['protein_sequence_key'])
+
+    if fetched_sequence and raw_sequence != fetched_sequence:
+        st.warning('Input protein sequence does not match fetched protein sequence from UniProt. '
+                   'Ensure that any changes made are deliberate')
 
     protein_sequence = raw_sequence.replace(' ', '').replace('\n', '')
+
+    if not protein_sequence:
+        st.error('Please enter a protein sequence')
+        st.stop()
+
     stripped_protein_sequence = strip_modifications(protein_sequence)
     protein_length = calculate_sequence_length(protein_sequence)
     st.caption(f'Length: {protein_length}')
@@ -76,10 +119,10 @@ def fetch_sequence_from_uniprot(accession_number):
     proteases_selected = c1.multiselect(label="Proteases",
                                         options=list(VALID_PROTEASES.keys()),
                                         help='The proteases to use for digestion',
-                                        default=DEFAULT_PROTEASES)
+                                        default=query_proteases)
 
     custom_regex = c2.text_input(label='(Additional) Custom protease',
-                                 value='',
+                                 value=query_custom_regex,
                                  help='A custom regular expression to use for digestion. Will be used along with '
                                       'selected proteases')
 
@@ -87,11 +130,12 @@ def fetch_sequence_from_uniprot(accession_number):
     missed_cleavages = c1.number_input(label='Max missed cleavages',
                                        min_value=MIN_MISSED_CLEAVAGES,
                                        max_value=MAX_MISSED_CLEAVAGES,
-                                       value=DEFAULT_MISSED_CLEAVAGES,
+                                       value=query_missed_cleavages,
                                        step=MISSED_CLEAVAGES_STEP,
                                        help='Number of missed cleavages to allow during digestion')
 
-    mass_type = c2.selectbox(label='Mass type', options=MASS_TYPE_OPTIONS, index=DEFAULT_MASS_TYPE_INDEX,
+    mass_type = c2.selectbox(label='Mass type', options=MASS_TYPE_OPTIONS,
+                             index=MASS_TYPE_OPTIONS.index(query_mass_type),
                              help='Mass type to use for calculations')
     is_mono = mass_type == 'monoisotopic'
 
@@ -104,14 +148,14 @@ def fetch_sequence_from_uniprot(accession_number):
     min_peptide_len = c1.number_input(label='Min length',
                                       min_value=MIN_PEPTIDE_LEN,
                                       max_value=MAX_PEPTIDE_LEN,
-                                      value=DEFAULT_MIN_PEPTIDE_LEN,
+                                      value=query_min_peptide_len,
                                       step=PEPTIDE_LEN_STEP,
                                       help='Minimum peptide length (inclusive)')
 
     max_peptide_len = c2.number_input(label='Max length',
                                       min_value=MIN_PEPTIDE_LEN,
                                       max_value=MAX_PEPTIDE_LEN,
-                                      value=DEFAULT_MAX_PEPTIDE_LEN,
+                                      value=query_max_peptide_len,
                                       step=PEPTIDE_LEN_STEP,
                                       help='Maximum peptide length (inclusive)')
 
@@ -123,14 +167,14 @@ def fetch_sequence_from_uniprot(accession_number):
     min_mass = c3.number_input(label='Min neutral mass',
                                min_value=MIN_PEPTIDE_MASS,
                                max_value=MAX_PEPTIDE_MASS,
-                               value=DEFAULT_MIN_PEPTIDE_MASS,
+                               value=query_min_mass,
                                step=PEPTIDE_MASS_STEP,
                                help='Minimum peptide neutral mass (inclusive)')
 
     max_mass = c4.number_input(label='Max neutral mass',
                                min_value=MIN_PEPTIDE_MASS,
                                max_value=MAX_PEPTIDE_MASS,
-                               value=DEFAULT_MAX_PEPTIDE_MASS,
+                               value=query_max_mass,
                                step=PEPTIDE_MASS_STEP,
                                help='Maximum peptide neutral mass (inclusive)')
 
@@ -141,8 +185,11 @@ def fetch_sequence_from_uniprot(accession_number):
     c1, c2 = st.columns(2)
 
     semi_enzymatic = c1.checkbox(label='Semi enzymatic?',
-                                 help='Allow semi enzymatic peptides?')
-    infer_charge = c2.checkbox(label='Infer charge', value=True, help='Infer charge of peptides based on L and K count')
+                                 help='Allow semi enzymatic peptides?',
+                                 value=query_semi_enzymatic)
+    infer_charge = c2.checkbox(label='Infer charge',
+                               value=query_infer_charge,
+                               help='Infer charge of peptides based on L and K count')
 
     min_charge, max_charge, min_mz, max_mz = None, None, None, None
     if infer_charge:
@@ -150,14 +197,14 @@ def fetch_sequence_from_uniprot(accession_number):
         min_charge = c1.number_input(label='Min charge',
                                      min_value=MIN_CHARGE,
                                      max_value=MAX_CHARGE,
-                                     value=DEFAULT_MIN_CHARGE,
+                                     value=query_min_charge,
                                      step=CHARGE_STEP,
                                      help='Minimum peptide charge (inclusive)')
 
         max_charge = c2.number_input(label='Max charge',
                                      min_value=MIN_CHARGE,
                                      max_value=MAX_CHARGE,
-                                     value=DEFAULT_MAX_CHARGE,
+                                     value=query_max_charge,
                                      step=CHARGE_STEP,
                                      help='Maximum peptide charge (inclusive)')
 
@@ -165,14 +212,14 @@ def fetch_sequence_from_uniprot(accession_number):
         min_mz = c1.number_input(label='Min m/z',
                                  min_value=MIN_MZ,
                                  max_value=MAX_MZ,
-                                 value=DEFAULT_MIN_MZ,
+                                 value=query_min_mz,
                                  step=MZ_STEP,
                                  help='Minimum peptide m/z (inclusive)')
 
         max_mz = c2.number_input(label='Max m/z',
                                  min_value=MIN_MZ,
                                  max_value=MAX_MZ,
-                                 value=DEFAULT_MAX_MZ,
+                                 value=query_max_mz,
                                  step=MZ_STEP,
                                  help='Maximum peptide m/z (inclusive)')
 
@@ -180,72 +227,140 @@ def fetch_sequence_from_uniprot(accession_number):
             st.error('Min charge must be less than max charge')
             st.stop()
 
-    remove_non_proteotypic = st.checkbox(label='Remove non-proteotypic peptides', value=False,
+    remove_non_proteotypic = st.checkbox(label='Remove non-proteotypic peptides',
+                                         value=query_remove_non_proteotypic,
                                          help='Remove peptides that are not proteotypic')
 
-    st.subheader('Static Modifications')
+    with st.expander('Static Modifications'):
+
+        c1, c2 = st.columns(2)
+        n_term_static_mod = c1.number_input(label='N-term mod',
+                                            value=query_n_term_static_mod,
+                                            help='Apply a static modification to the N-terminus')
+
+        c_term_static_mod = c2.number_input(label='C-term mod',
+                                            value=query_c_term_static_mod,
+                                            help='Apply a static modification to the C-terminus')
+
+        num_static_mods = st.number_input(label='Number of unique static modifications',
+                                          min_value=MIN_STATIC_MODS,
+                                          max_value=MAX_STATIC_MODS,
+                                          value=query_num_static_mods,
+                                          step=STATIC_MODS_STEP,
+                                          help='Add another modification row')
+
+        # columns to lay out the inputs
+        grid = st.columns([3, 2])
 
-    # a selection for the user to specify the number of rows
-    c1, c2 = st.columns(2)
-    num_rows = st.number_input(label='Number of unique modifications',
-                               min_value=MIN_STATIC_MODS,
-                               max_value=MAX_STATIC_MODS,
-                               value=DEFAULT_STATIC_MODS,
-                               step=STATIC_MODS_STEP,
-                               help='Add another modification row')
-
-    # columns to lay out the inputs
-    grid = st.columns([3, 2])
-
-
-    def add_row(r):
-        with grid[0]:
-            st.multiselect(label='Amino acids',
-                           key=f'residues{r}',
-                           options=list(AMINO_ACIDS),
-                           help='Select amino acids for which to apply the static modification',
-                           default=['C'] if r == 0 else [])
-        with grid[1]:
-            st.number_input(label='Modification Mass (Da)',
-                            step=0.00001, key=f'mass{r}',
-                            help='The mass of the modification (in daltons)',
-                            value=57.02146 if r == 0 else 0.0,
-                            format='%.5f')
-
-
-    # Loop to create rows of input widgets
-    for r in range(num_rows):
-        add_row(r)
-
-    mods = {}
-    for r in range(num_rows):
-        for residue in st.session_state[f'residues{r}']:
-            mods[residue] = "{:.5f}".format(st.session_state[f'mass{r}'])
+
+        def add_static_modification(r):
+
+            aas = list(query_static_mods[r][0]) if r < len(query_static_mods) else []
+            mod = query_static_mods[r][1] if r < len(query_static_mods) else 0.0
+
+            with grid[0]:
+                st.multiselect(label='Amino acids',
+                               key=f'static_mod_residue{r}',
+                               options=list(AMINO_ACIDS),
+                               help='Select amino acids for which to apply the static modification',
+                               default=aas)
+            with grid[1]:
+                st.number_input(label='Modification Mass (Da)',
+                                step=0.00001,
+                                key=f'static_mod_mass{r}',
+                                help='The mass of the modification (in daltons)',
+                                value=mod,
+                                format='%.5f')
+
+
+        # Loop to create rows of input widgets
+        for r in range(num_static_mods):
+            add_static_modification(r)
+
+        static_mods = {}
+        for r in range(num_static_mods):
+            mod = "{:.5f}".format(st.session_state[f'static_mod_mass{r}'])
+            for residue in st.session_state[f'static_mod_residue{r}']:
+                static_mods[residue] = mod
+
+    with st.expander('Variable Modifications'):
+
+        c1, c2 = st.columns(2)
+        n_term_var_mod = c1.number_input(label='N-term var mod',
+                                         value=query_n_term_var_mod,
+                                         help='Apply a variable modification to the N-terminus')
+        c_term_var_mod = c2.number_input(label='C-term var mod',
+                                         value=query_c_term_var_mod,
+                                         help='Apply a variable modification to the C-terminus')
+
+        max_var_mods = st.number_input(label='Max var mods',
+                                       min_value=MIN_MAX_VAR_MODS,
+                                       max_value=MAX_MAX_VAR_MODS,
+                                       value=query_max_var_mods,
+                                       step=MAX_VAR_MOD_STEP,
+                                       help='Maximum number of variable modifications per peptide')
+
+        num_variable_mods = st.number_input(label='Unique var mods',
+                                            min_value=MIN_VAR_MODS,
+                                            max_value=MAX_VAR_MODS,
+                                            value=query_num_variable_mods,
+                                            step=VAR_MOD_STEP,
+                                            help='Add another modification row')
+
+        # columns to lay out the inputs
+        grid = st.columns([3, 2])
+
+        def add_variable_modification(r):
+            aas = list(query_variable_mods[r][0]) if r < len(query_variable_mods) else []
+            mod = query_variable_mods[r][1] if r < len(query_variable_mods) else 0.0
+
+            with grid[0]:
+                st.multiselect(label='Amino acids',
+                               key=f'var_mod_residue{r}',
+                               options=list(AMINO_ACIDS),
+                               help='Select amino acids for which to apply the variable modification',
+                               default=aas)
+            with grid[1]:
+                st.number_input(label='Modification Mass (Da)',
+                                step=0.00001, key=f'var_mod_mass{r}',
+                                help='The mass of the modification (in daltons)',
+                                value=mod,
+                                format='%.5f')
+
+
+        # Loop to create rows of input widgets
+        for r in range(num_variable_mods):
+            add_variable_modification(r)
+
+        var_mods = {}
+        for r in range(num_variable_mods):
+            mod = "{:.5f}".format(st.session_state[f'var_mod_mass{r}'])
+            for residue in st.session_state[f'var_mod_residue{r}']:
+                var_mods[residue] = mod
+
+    url = generate_app_url(protein_id, protein_sequence, proteases_selected, custom_regex, missed_cleavages, mass_type,
+                           min_peptide_len, max_peptide_len, min_mass, max_mass, semi_enzymatic, infer_charge,
+                           min_charge, max_charge, min_mz, max_mz, remove_non_proteotypic, n_term_static_mod,
+                           c_term_static_mod, num_static_mods, n_term_var_mod, c_term_var_mod, max_var_mods,
+                           num_variable_mods, static_mods, var_mods)
 
 sites = set()
 for enzyme_regex in enzyme_regexes:
     sites.update(identify_cleavage_sites(protein_sequence, enzyme_regex))
 sites = sorted(list(sites))
 
-df = generate_peptide_df(protein_sequence, sites, missed_cleavages, min_peptide_len, max_peptide_len,
-                         semi_enzymatic, mods, min_mass, max_mass, is_mono, infer_charge, min_charge,
-                         max_charge, min_mz, max_mz)
+with st.expander('Edit Sites'):
+    sites = st.multiselect(label="Sites",
+                            options=list(range(len(stripped_protein_sequence)+1)),
+                            help='The proteases to use for digestion',
+                            default=sites)
 
-rt_model = pickle.load(open("rt_model.pkl", "rb"))
-df['RT'] = rt_model.predict(np.array([bin_aa_counts(strip_modifications(seq)) for seq in df['Sequence']]))
-df['RT'] = df['RT'].round(3)
 
-im_model = pickle.load(open("im_model.pkl", "rb"))
-df['IM'] = im_model.predict(
-    np.array([bin_aa_counts(strip_modifications(seq), c) for seq, c in df[['Sequence', 'Charge']].values]))
-df['IM'] = df['IM'].round(3)
 
-proteotypic_model = pickle.load(open("proteotypic_model.pkl", "rb"))
-df['Proteotypic'] = proteotypic_model.predict(
-    np.array([bin_aa_counts(strip_modifications(seq)) for seq in df['Sequence']])).astype(bool)
-
-if remove_non_proteotypic:
-    df = df[df['Proteotypic']]
+df = generate_peptide_df(protein_sequence, sites, missed_cleavages, min_peptide_len, max_peptide_len,
+                         semi_enzymatic, static_mods, min_mass, max_mass, is_mono, infer_charge, min_charge,
+                         max_charge, min_mz, max_mz, var_mods, max_var_mods, n_term_static_mod, c_term_static_mod,
+                         n_term_var_mod, c_term_var_mod, remove_non_proteotypic)
 
 # Start the HTML string for the site indexes
 site_indexes_html = '<span style="font-family: Courier New, monospace; font-size: 16px;">'
@@ -303,10 +418,7 @@ def add_row(r):
     cov = calculate_span_coverage(spans, protein_length)
     protein_cov_at_mass.append(sum(cov) / len(cov) * 100)
 
-df.drop(columns=['StrippedPeptide'], inplace=True)
-df.sort_values(by=['MC'], inplace=True)
-df.drop_duplicates(subset=['Start', 'Sequence', 'Semi'], inplace=True)
-df.sort_values(by=['Start', 'Len'], inplace=True)
+st.write(f'##### [Analysis URL]({url}) (copy me and send to your friends!)')
 
 t1, t2, t3, t5 = st.tabs(['Digestion Metrics', 'Cleavage & Coverage', 'Motif Analysis', 'Help'])
 
@@ -367,34 +479,35 @@ def add_row(r):
 
 with t3:
     st.header('Motif Analysis')
-    site_regex = st.text_input('Motifs Regex', '(K)')
+    motif_regex = st.text_input('Motifs Regex', '(K)')
+
+    if motif_regex:
+        motif_sites = list(reg.finditer(motif_regex, stripped_protein_sequence, overlapped=True))
 
-    if site_regex:
-        sites = identify_cleavage_sites(stripped_protein_sequence, site_regex)
+        def count_motifs(row):
+            return sum([1 for site in motif_sites if row['Start'] <= site.start() < row['End']])
 
-        site_counts = []
-        for row in df[['Start', 'End']].values:
-            site_counts.append(sum([1 for site in sites if row[0] <= site < row[1]]))
-        df['Motifs'] = site_counts
+        df['Motifs'] = df.apply(count_motifs, axis=1)
+        motif_cov_indexes = {i for site in motif_sites for i in range(site.start(), site.end())}
 
-        cov_site_mat = [0] * len(stripped_protein_sequence)
+        motif_cov_array = [0] * len(stripped_protein_sequence)
         for row in df[['Start', 'End', 'Motifs']].values:
 
             if row[2] == 0:
                 continue
 
             for i in range(row[0], row[1] + 1):
-                if i in sites:
+                if i in motif_cov_indexes:
                     if row[2] == 1:
-                        cov_site_mat[i - 1] = 1
+                        motif_cov_array[i] = 1
                     else:
-                        if cov_site_mat[i - 1] == 0:
-                            cov_site_mat[i - 1] = row[2]
+                        if motif_cov_array[i] == 0:
+                            motif_cov_array[i] = row[2]
                         else:
-                            cov_site_mat[i - 1] = min(row[2], cov_site_mat[i - 1])
+                            motif_cov_array[i] = min(row[2], motif_cov_array[i])
 
-        clickable2 = st.checkbox('Peptide Fragmenter Links', value=False, key=1)
         st.subheader('Peptides')
+        clickable2 = st.checkbox('Peptide Fragmenter Links', value=False, key=1)
         if clickable2:
             df_clickable = df.copy(deep=True)
             df_clickable['Sequence'] = [make_clickable(peptide, mass_type) for peptide in df_clickable['Sequence']]
@@ -403,13 +516,15 @@ def add_row(r):
         else:
             st.dataframe(df, use_container_width=True)
 
-        counter = Counter(site_counts)
+        counter = Counter(df['Motifs'])
 
-        st.subheader('Coverage Analysis', help='Coverage of protein based on peptides with N number of motif matches')
-        sites = [i - 1 for i in sites if i != 0]
-        protein_coverage = coverage_string(cov_site_mat, stripped_protein_sequence, cmap, sites)
+        st.subheader('Motif Site Coverage', help='The color corresponds to the peptide with the fewest number of motif '
+                                                 'matches (excluding 0 matches). Example: Lets assume that the first '
+                                                 'site is covered by two peptides, the first with one match and the '
+                                                 'second with two matches, then this site will be displayed as 1.')
+        protein_coverage = coverage_string(motif_cov_array, stripped_protein_sequence, cmap, motif_cov_indexes)
         st.markdown(protein_coverage, unsafe_allow_html=True)
-        f = create_colorbar(max(cov_site_mat), cmap, label='Min Number of Motif Matches')
+        f = create_colorbar(max(motif_cov_array), cmap, label='Min Number of Motif Matches')
         st.pyplot(f)
         with st.expander('Show Coverage Table'):
             for i, (k, v) in enumerate(sorted(counter.items())):
@@ -433,11 +548,10 @@ def add_row(r):
                     st.warning('Warning: High motif match counts may result in long runtimes. Stopping...')
                     break
 
-#with t4:
+# with t4:
 #    st.markdown(PROTEASE_WIKI)
 
 with t5:
-
     st.header('Help')
 
     st.subheader('General')
@@ -450,7 +564,7 @@ def add_row(r):
 
     with st.expander('Protease Regexes'):
         st.subheader('Protease Regexes')
-        data = [{'Name': k, 'Regex': v} for k, v in PROTEASES.items()]
+        data = [{'Name': k, 'Regex': v} for k, v in VALID_PROTEASES.items()]
         protease_df = pd.DataFrame(data)
         st.table(protease_df)
 
@@ -469,7 +583,6 @@ def add_row(r):
         st.markdown(PROTEOTYPIC_MODEL_HELP)
 
     with st.expander('How to use ML Models?'):
-
         def get_model_file_as_byte_stream(path):
             with open(path, 'rb') as file:
                 byte_stream = file.read()
@@ -500,5 +613,3 @@ def get_model_file_as_byte_stream(path):
 
         st.subheader('Example Code')
         st.code(MODEL_CODE)
-
-
diff --git a/constants.py b/constants.py
index dd7446a..064703f 100644
--- a/constants.py
+++ b/constants.py
@@ -6,7 +6,7 @@
 MAX_PROTEIN_INPUT_LENGTH = 10_000
 
 MIN_MISSED_CLEAVAGES = 0
-MAX_MISSED_CLEAVAGES = 5
+MAX_MISSED_CLEAVAGES = 7
 DEFAULT_MISSED_CLEAVAGES = 1
 MISSED_CLEAVAGES_STEP = 1
 
@@ -33,12 +33,22 @@
 DEFAULT_STATIC_MODS = 1
 STATIC_MODS_STEP = 1
 
+MIN_VAR_MODS = 0
+MAX_VAR_MODS = 10
+DEFAULT_VAR_MODS = 0
+VAR_MOD_STEP = 1
+
 MIN_CHARGE = 1
 MAX_CHARGE = 10
 DEFAULT_MIN_CHARGE = 2
 DEFAULT_MAX_CHARGE = 5
 CHARGE_STEP = 1
 
+MIN_MAX_VAR_MODS = 0
+MAX_MAX_VAR_MODS = 3
+DEFAULT_MAX_VAR_MODS = 0
+MAX_VAR_MOD_STEP = 1
+
 MIN_MZ = 0.0
 MAX_MZ = 10_000.0
 DEFAULT_MIN_MZ = 200.0
@@ -54,5 +64,6 @@
 
 VALID_PROTEASES = copy.deepcopy(PROTEASES)
 VALID_PROTEASES.pop('non-specific', None)
+VALID_PROTEASES = {k.replace(' ', '_'): v for k, v in VALID_PROTEASES.items()}
 
 
diff --git a/util.py b/util.py
index 8f361dc..c6b1997 100644
--- a/util.py
+++ b/util.py
@@ -1,17 +1,28 @@
+import pickle
 from collections import Counter
 from typing import List
 
+import numpy as np
 import pandas as pd
+import requests
 from peptacular.mass import calculate_mass
 from peptacular.sequence import span_to_sequence, calculate_sequence_length, apply_static_modifications, \
-    strip_modifications
+    strip_modifications, apply_variable_modifications
 from peptacular.spans import build_enzymatic_spans, build_semi_spans
 import matplotlib.pyplot as plt
 import matplotlib.colors as mcolors
 import matplotlib as mpl
+from peptacular.term import add_n_term_modification, add_c_term_modification
 
 from constants import LINK
 
+
+def fetch_sequence_from_uniprot(accession_number):
+    url = f"https://www.uniprot.org/uniprot/{accession_number}.fasta"
+    response = requests.get(url)
+    return response
+
+
 def make_clickable(sequence, mass_type):
     # target _blank to open new window
     # extract clickable text to display for your link
@@ -22,7 +33,9 @@ def make_clickable(sequence, mass_type):
 def generate_peptide_df(sequence: str, cleavage_sites: List, missed_cleavages: int, min_len: int,
                         max_len: int, semi_enzymatic: bool, static_mods: dict, min_mass: float, max_mass: float,
                         is_mono: bool, infer_charge: bool, min_charge: int, max_charge: int, min_mz: float,
-                        max_mz: float):
+                        max_mz: float, var_mods: dict, max_var_mods: int, n_term_static_mod: float,
+                        c_term_static_mod: float, n_term_var_mod: float, c_term_var_mod: float,
+                        remove_non_proteotypic: bool):
     cleavage_sites = sorted(cleavage_sites)
     spans = build_enzymatic_spans(calculate_sequence_length(sequence), cleavage_sites, missed_cleavages, 1, None)
     df = pd.DataFrame(spans, columns=['Start', 'End', 'MC'])
@@ -46,12 +59,61 @@ def generate_peptide_df(sequence: str, cleavage_sites: List, missed_cleavages: i
 
     df = df[(df['Len'] >= min_len) & (df['Len'] <= max_len)]
 
-    df['Sequence'] = df['Sequence'].apply(lambda x: apply_static_modifications(x, static_mods))
+    # Apply variable modifications to each sequence in the DataFrame
+    def apply_var_mods(sequence):
+
+        var_seqs = apply_variable_modifications(sequence, var_mods, max_var_mods)
+
+        if n_term_var_mod and c_term_var_mod:
+            n_term_seq = add_n_term_modification(sequence, n_term_var_mod)
+            c_term_seq = add_c_term_modification(sequence, c_term_var_mod)
+            n_c_term_seq = add_c_term_modification(n_term_seq, c_term_var_mod)
+
+            n_term_seqs = apply_variable_modifications(n_term_seq, var_mods, max_var_mods)
+            c_term_seqs = apply_variable_modifications(c_term_seq, var_mods, max_var_mods)
+            n_c_term_seqs = apply_variable_modifications(n_c_term_seq, var_mods, max_var_mods)
+
+            return ';'.join(list(set(var_seqs + n_term_seqs + c_term_seqs + n_c_term_seqs)))
+
+        elif n_term_var_mod:
+            n_term_seq = add_n_term_modification(sequence, n_term_var_mod)
+            n_term_seqs = apply_variable_modifications(n_term_seq, var_mods, max_var_mods)
+            return ';'.join(list(set(var_seqs + n_term_seqs)))
+        elif c_term_var_mod:
+            c_term_seq = add_c_term_modification(sequence, c_term_var_mod)
+            c_term_seqs = apply_variable_modifications(c_term_seq, var_mods, max_var_mods)
+            return ';'.join(list(set(var_seqs + c_term_seqs)))
+        else:
+            return ';'.join(var_seqs)
+
+    # Apply the apply_var_mods function to the 'Sequence' column
+    df['Sequence'] = df['Sequence'].apply(apply_var_mods)
+
+    # expand the sequence column into multiple rows (sequences are separated by ';')
+    df = df.assign(Sequence=df.Sequence.str.split(';')).explode('Sequence')
+
+    def apply_static_mods(sequence):
+        # Update var_mods dictionary based on conditions
+        if n_term_static_mod:
+            var_mods.update({-1: n_term_static_mod})
+        if c_term_static_mod:
+            var_mods.update({calculate_sequence_length(sequence): c_term_static_mod})
+
+        # Apply variable modifications and join them with ';'
+        return apply_static_modifications(sequence, static_mods)
+
+    df['Sequence'] = df['Sequence'].apply(apply_static_mods)
+
+    # drop duplicates
 
     df['NeutralMass'] = [round(calculate_mass(sequence, monoisotopic=is_mono), 5) for sequence in df['Sequence']]
     df = df[(df['NeutralMass'] >= min_mass) & (df['NeutralMass'] <= max_mass)]
     df['StrippedPeptide'] = df['Sequence'].apply(strip_modifications)
 
+    df.sort_values(by=['MC'], inplace=True)
+    df.drop_duplicates(subset=['Start', 'Sequence'], inplace=True)
+    df.sort_values(by=['Start', 'Len'], inplace=True)
+
     if infer_charge is True:
         # charge should be the Lysine and Arginine count + 1
         df['Charge'] = df['StrippedPeptide'].apply(lambda x: x.count('K') + x.count('R') + 1)
@@ -59,6 +121,26 @@ def generate_peptide_df(sequence: str, cleavage_sites: List, missed_cleavages: i
         df = df[(df['Charge'] >= min_charge) & (df['Charge'] <= max_charge)]
         df = df[(df['Mz'] >= min_mz) & (df['Mz'] <= max_mz)]
 
+    if len(df) == 0:
+        return df
+
+    rt_model = pickle.load(open("rt_model.pkl", "rb"))
+    df['RT'] = rt_model.predict(np.array([bin_aa_counts(strip_modifications(seq)) for seq in df['Sequence']]))
+    df['RT'] = df['RT'].round(3)
+
+    if infer_charge:
+        im_model = pickle.load(open("im_model.pkl", "rb"))
+        df['IM'] = im_model.predict(
+            np.array([bin_aa_counts(strip_modifications(seq), c) for seq, c in df[['Sequence', 'Charge']].values]))
+        df['IM'] = df['IM'].round(3)
+
+    proteotypic_model = pickle.load(open("proteotypic_model.pkl", "rb"))
+    df['Proteotypic'] = proteotypic_model.predict(
+        np.array([bin_aa_counts(strip_modifications(seq)) for seq in df['Sequence']])).astype(bool)
+
+    if remove_non_proteotypic:
+        df = df[df['Proteotypic']]
+
     return df
 
 
@@ -99,6 +181,7 @@ def coverage_string(protein_cov_arr, stripped_protein_sequence, cmap, sites=None
 
     return protein_coverage
 
+
 def create_colorbar(max_coverage, cmap, label='Coverage'):
     # Create a figure and a subplot for the colorbar
     fig, ax = plt.subplots(figsize=(15, 1))
@@ -112,4 +195,54 @@ def create_colorbar(max_coverage, cmap, label='Coverage'):
     cb.set_label(label)
 
     # Display the colorbar in Streamlit
-    return fig
\ No newline at end of file
+    return fig
+
+
+def generate_app_url(protein_id, protein_sequence, proteases, custom_regex, missed_cleavages, mass_type,
+                     min_peptide_len, max_peptide_len, min_mass, max_mass, semi_enzymatic, infer_charge,
+                     min_charge, max_charge, min_mz, max_mz, remove_non_proteotypic, n_term_static_mod,
+                     c_term_static_mod, num_static_mods, n_term_var_mod, c_term_var_mod, max_var_mods,
+                     num_variable_mods, static_mods, variable_mods):
+    base_url = 'http://localhost:8501/'
+
+    # flip the dictionary
+    static_mods_rev = {}
+    for aa, mod in static_mods.items():
+        static_mods_rev.setdefault(mod, []).append(aa)
+    static_mod_str = ';'.join([f'{"".join(aa)}:{mod}' for mod, aa in static_mods_rev.items()])
+
+    variable_mods_rev = {}
+    for aa, mod in variable_mods.items():
+        variable_mods_rev.setdefault(mod, []).append(aa)
+    variable_mod_str = ';'.join([f'{"".join(aa)}:{mod}' for mod, aa in variable_mods_rev.items()])
+
+    params = {
+        'protein_id': protein_id,
+        'protein_sequence': protein_sequence,
+        'proteases': ';'.join(proteases),
+        'custom_regex': custom_regex,
+        'missed_cleavages': missed_cleavages,
+        'mass_type': mass_type,
+        'min_peptide_len': min_peptide_len,
+        'max_peptide_len': max_peptide_len,
+        'min_mass': min_mass,
+        'max_mass': max_mass,
+        'semi_enzymatic': semi_enzymatic,
+        'infer_charge': infer_charge,
+        'min_charge': min_charge,
+        'max_charge': max_charge,
+        'min_mz': min_mz,
+        'max_mz': max_mz,
+        'remove_non_proteotypic': remove_non_proteotypic,
+        'n_term_static_mod': n_term_static_mod,
+        'c_term_static_mod': c_term_static_mod,
+        'num_static_mods': num_static_mods,
+        'n_term_var_mod': n_term_var_mod,
+        'c_term_var_mod': c_term_var_mod,
+        'max_var_mods': max_var_mods,
+        'num_variable_mods': num_variable_mods,
+        'static_mods': static_mod_str,
+        'variable_mods': variable_mod_str
+    }
+    query_string = '&'.join([f'{key}={value}' for key, value in params.items() if value is not None])
+    return f'{base_url}?{query_string}'