Helper_Functions.py

from HXMS_IO import HxmsData
import csv
from typing import Optional
import zipfile
from pigeon_feather.spectra import *
from pigeon_feather.tools import custom_pad
from pigeon_feather.hxio import *
import re

def _create_peptide_chemical_formula(peptide:str)->dict:
    total_formula = {}
    residue_formulas = {
        'A': {'C': 3, 'H': 5, 'N': 1, 'O': 1},  # Alanine
        'R': {'C': 6, 'H': 12, 'N': 4, 'O': 1},  # Arginine
        'N': {'C': 4, 'H': 6, 'N': 2, 'O': 2},  # Asparagine
        'D': {'C': 4, 'H': 5, 'N': 1, 'O': 3},  # Aspartic Acid
        'C': {'C': 3, 'H': 5, 'N': 1, 'O': 1, 'S': 1},  # Cysteine
        'Q': {'C': 5, 'H': 8, 'N': 2, 'O': 2},  # Glutamine
        'E': {'C': 5, 'H': 7, 'N': 1, 'O': 3},  # Glutamic Acid
        'G': {'C': 2, 'H': 3, 'N': 1, 'O': 1},  # Glycine
        'H': {'C': 6, 'H': 7, 'N': 3, 'O': 1},  # Histidine
        'I': {'C': 6, 'H': 11, 'N': 1, 'O': 1},  # Isoleucine
        'L': {'C': 6, 'H': 11, 'N': 1, 'O': 1},  # Leucine
        'K': {'C': 6, 'H': 12, 'N': 2, 'O': 1},  # Lysine
        'M': {'C': 5, 'H': 9, 'N': 1, 'O': 1, 'S': 1},  # Methionine
        'F': {'C': 9, 'H': 9, 'N': 1, 'O': 1},  # Phenylalanine
        'P': {'C': 5, 'H': 7, 'N': 1, 'O': 1},  # Proline
        'S': {'C': 3, 'H': 5, 'N': 1, 'O': 2},  # Serine
        'T': {'C': 4, 'H': 7, 'N': 1, 'O': 2},  # Threonine
        'W': {'C': 11, 'H': 10, 'N': 2, 'O': 1},  # Tryptophan
        'Y': {'C': 9, 'H': 9, 'N': 1, 'O': 2},  # Tyrosine
        'V': {'C': 5, 'H': 9, 'N': 1, 'O': 1}  # Valine
    }
    water_formula = {'H': 2, 'O': 1}
    for AA in peptide.upper():
        if AA not in residue_formulas:
            continue
        residue = residue_formulas[AA]
        for element, count in residue.items():
            if element not in total_formula:
                total_formula[element] = 0
            total_formula[element] = total_formula[element] + count
    for element, count in water_formula.items():
        total_formula[element] = total_formula[element] + count
    return total_formula

def _convert_dformula_to_dict(dformula:str)->dict:
    pattern = re.compile(r'([A-Z][a-z]*)(\d*)')
    initial_counts = {}
    matches = pattern.findall(dformula)
    for element, count in matches:
        count = int(count) if count else 1
        initial_counts[element] = initial_counts.get(element, 0) + count
    final_formula = initial_counts.copy()
    if 'D' in final_formula:
        deuterium_count = final_formula.pop('D')
        final_formula['H'] = final_formula.get('H', 0) + deuterium_count

    return final_formula

def _subtract_chemical_formula_dicts(form1:dict,form2:dict)->str:
    diff = ""
    used = []
    for key in form1:
        if key in form2:
            used.append(key)
            difference = math.fabs(form1[key]-form2[key])
            if difference != 0:
                diff = diff + key + str(int(math.fabs(form1[key]-form2[key])))
        else:
            diff = diff + key + str(form1[key])
    for key in form2:
        if key not in used:
            diff = diff + key + str(form2[key])
    return diff


def _conver_PFhxms_to_hxms(dataset_file_path,protein_sequence,saturation,ph,temperature,envelope_file_path,protein_name,protein_state:str=None):
    hxms = HxmsData()
    hdxms_data_list = []
    cleaned = read_hdx_tables([dataset_file_path], [dataset_file_path], exclude=False, states_subset=[protein_state])
    hdxms_data = load_dataframe_to_hdxmsdata(
        cleaned,
        n_fastamides=0,
        protein_sequence=protein_sequence,
        saturation=saturation,
        pH=ph,
        temperature=temperature,
    )
    import zipfile
    with zipfile.ZipFile(envelope_file_path, 'r') as zip_ref:
        zip_ref.extractall(os.path.dirname(os.path.normpath(envelope_file_path)))
    raw_spectra_path = os.path.join(os.path.dirname(os.path.normpath(envelope_file_path)),
                                    os.path.splitext(os.path.basename(envelope_file_path))[0])

    load_raw_ms_to_hdxms_data(
        hdxms_data,
        raw_spectra_path,
        protein_state,
        save_match=True,
        back_exchange_threshold=1.0, # no filtering
        refine_peptides=False
    )

    hdxms_data_list.append(hdxms_data)
    state_names = set([state.state_name for data in hdxms_data_list for state in data.states])
    if protein_name not in hxms.proteins:
        hxms.proteins.append(protein_name)
        hxms.state[protein_name] = []
        hxms.peptides[protein_name] = {}

    for state_name in state_names:
        protein_states = [state for data in hdxms_data_list for state in data.states if
                          state.state_name == state_name]
        all_peptides = [pep for state in protein_states for pep in state.peptides]
        all_timepoints = [tp for pep in all_peptides for tp in pep.timepoints]
        tp_dict = {}
        for tp_idx, tp in enumerate(all_timepoints):
            tp_dict[
                protein_name + "_" + state_name + "_" + str(tp.peptide.start) + "_" + str(tp.peptide.end) + "_" + str(
                    tp.deut_time)] = \
                tp.isotope_envelope is not None
        if state_name not in hxms.state[protein_name]:
            hxms.state[protein_name].append(state_name)
            hxms.peptides[protein_name][state_name] = {}

        for tp_idx, tp in enumerate(all_timepoints):
            peptide_name = f"{tp.peptide.start}~{tp.peptide.sequence}~N/A~{tp.peptide.end}~0"
            if peptide_name not in hxms.peptides[protein_name][state_name]:
                hxms.peptides[protein_name][state_name][peptide_name] = {}
            if tp.isotope_envelope is None and tp.deut_time != np.inf:
                print(f"No envelope for {tp.peptide.identifier} {tp.deut_time}")
                continue
            key = protein_name + "_" + state_name + "_" + str(tp.peptide.start) + "_" + str(tp.peptide.end) + "_" + str(0.0)
            if key in tp_dict:
                if not tp_dict[key]:
                    continue
            else:
                continue
            if str(tp.deut_time) not in hxms.peptides[protein_name][state_name][peptide_name]:
                hxms.peptides[protein_name][state_name][peptide_name][str(tp.deut_time)] = [
                 {
                "uptake": float(tp.num_d),
                "start": int(tp.peptide.start),
                "end": int(tp.peptide.end),
                "sequence": str(tp.peptide.sequence),
                "PTM": "0000",
                "time": float(tp.deut_time),
                "mod": "A",
                "envelope": f"{','.join(f'{val:.3f}' for val in custom_pad(tp.isotope_envelope[:50].copy(), 50)) if tp.deut_time != np.inf else ''}",
                "replicate": 0,
                "mono_mz": f"{tp.mono_mz:.6f}" if tp.deut_time != np.inf else '',
                "charge_state": tp.charge_state,
                "RT": f"{tp.peptide.RT:.3f}",
                "match": f"{','.join(f'{val[0]:.6f}:{val[1]:.6f}' for val in tp.match)}" if tp.deut_time != np.inf else '',
                "raw_ms_fine_structure": ",".join(";".join(f"{m:.6f}:{i:.6f}" for m, i in zip(*val)) for val in tp.raw_ms_fine_structure) if tp.deut_time != np.inf else '',
                "score": f"{tp.score:.2f}"
            }]  
                # if tp.deut_time != np.inf:
                # print(tp.peptide.identifier, tp.deut_time)
                # print(tp.raw_ms_fine_structure)

            else:
                hxms.peptides[protein_name][state_name][peptide_name][str(tp.deut_time)].append(
                    {
                        "uptake": float(tp.num_d),
                        "start": int(tp.peptide.start),
                        "end": int(tp.peptide.end),
                        "sequence": str(tp.peptide.sequence),
                        "PTM": "0000",
                        "time": float(tp.deut_time),
                        "mod": "A",
                        "envelope": f"{','.join(f'{val:.3f}' for val in custom_pad(tp.isotope_envelope[:50].copy(), 50)) if tp.deut_time != np.inf else ''}",
                        "replicate": 0,
                        "mono_mz": f"{tp.mono_mz:.6f}" if tp.deut_time != np.inf else '',
                        "charge_state": tp.charge_state,
                        "RT": f"{tp.peptide.RT:.3f}",
                        "match": f"{','.join(f'{val[0]:.6f}:{val[1]:.6f}' for val in tp.match)}" if tp.deut_time != np.inf else '',
                        "raw_ms_fine_structure": ",".join(";".join(f"{m:.6f}:{i:.6f}" for m, i in zip(*val)) for val in tp.raw_ms_fine_structure) if tp.deut_time != np.inf else '',
                        "score": f"{tp.score:.2f}"
                    })

    return hxms

def _check_peptide_list_format(file_path: str)->bool:
    if not os.path.exists(file_path) or not os.path.isfile(file_path):
        return False
    peptide_df = pd.read_csv(file_path)
    peptide_df.columns = peptide_df.columns.str.lower()
    if "start" not in peptide_df.columns:
        return False
    if "end" not in peptide_df.columns:
        return False
    return True



def _check_envelope_file(zip_file_path: str)->bool:
    if not os.path.exists(zip_file_path) or not os.path.isfile(zip_file_path):
        return False
    protein_states = []

    try:
        with zipfile.ZipFile(zip_file_path, 'r') as zip_ref:
            all_zip_entries = zip_ref.namelist()
            for entry in all_zip_entries:
                if entry.endswith('.csv'):
                    if ".csv" not in entry:
                        continue
                    folder_path = os.path.dirname(os.path.normpath(entry))
                    folder_name = os.path.basename(folder_path)
                    if len(folder_name.split("-")) != 3:
                        continue
                    state_path = os.path.dirname(os.path.normpath(folder_path))
                    state_name = os.path.basename(os.path.normpath(state_path))
                    if state_name not in protein_states:
                        protein_states.append(state_name)

    except zipfile.BadZipFile:
        print(f"Error: The file '{zip_file_path}' is not a valid zip file.")
        return False
    except Exception as e:
        print(f"An unexpected error occurred: {e}")
        return False
    if len(protein_states) != 0:
        return True
    return False
def _sort_hxms_data_cc(hx_data: HxmsData) -> HxmsData:

    hx_data.proteins.sort()
    for protein in hx_data.state:
        hx_data.state[protein].sort()

    for protein in hx_data.peptides:
        for state in hx_data.peptides[protein]:
            peptide_data = hx_data.peptides[protein][state]
            for peptide_name, time_data in peptide_data.items():
                for time_key in time_data:
                    time_data[time_key].sort(
                        key=lambda measurement: (
                            measurement['replicate'],
                            measurement['time']
                        )
                    )
                peptide_data[peptide_name] = dict(sorted(
                    time_data.items(),
                    key=lambda item: float('inf') if item[0] == 'inf' else float(item[0])
                ))
            hx_data.peptides[protein][state] = dict(sorted(
                peptide_data.items(),
                key=lambda item: (
                    int(float(item[0].split('~')[0])),
                    int(float(item[0].split('~')[-1])),
                    item[0].split('~')[2],
                )
            ))

    return hx_data



def _get_avg_zero_uptake(hxms_data: HxmsData, protein: str, state: str, peptide_name: str) -> Optional[float]:

    zero_uptakes = []
    base_peptide_name = "~".join(peptide_name.split("~")[:-1])
    for p_name, p_data in hxms_data.peptides[protein][state].items():
        if "~".join(p_name.split("~")[:-1]) == base_peptide_name:
            if "0" in p_data:
                zero_uptakes.append(p_data["0"][0]["uptake"])
            if "0.0" in p_data:
                zero_uptakes.append(p_data["0.0"][0]["uptake"])
    if zero_uptakes:
        return sum(zero_uptakes) / len(zero_uptakes)
    return None


def check_file_format(file_path: str) -> str:

    known_formats = {
        "HDXworkbench": [
            "peptide", "charge", "start", "end", "experiment",
            "sample", "timepoint", "replicate", "discarded_replicate",
            "percentd_replicate", "centroid"
        ],
        "DynamX": [
            "Protein", "Start", "End", "Sequence", "Modification", "Fragment",
            "MaxUptake", "MHP", "State", "Exposure", "Center", "Center SD",
            "Uptake", "Uptake SD", "RT", "RT SD"
        ],


        "DynamX_short": [
            "Protein", "Start", "End", "Sequence", "Modification", "Fragment", "MaxUptake", "MHP", "State",
            "Exposure", "File","z","RT","Inten","Center"
        ],

        "BioPharma": [
            "Retention Time", "Peptide", "Charge", "Residue From", "Residue To", "Ref", "0% Control",
            "100% Control", "stdev"
        ],
        "HDExaminer":[
            "State","Protein","Sequence","Search RT","Charge","Max D"
        ],
        "Custom": [
            "Protein","State","Start","Stop","Replicate","Modification","Bimodal_group","Exposure","Uptake","Envelope","RT","Conf","Z","Mono_M","M/Z"
        ]
    }

    if not os.path.exists(file_path):
        return "unknown"
    header_row_index = 0
    with open(file_path, 'r', newline='') as file:
        reader = csv.reader(file)
        for i, row in enumerate(reader):
            if any(h in row for h in ["peptide", "Protein", "Condition","State","Peptide"]):
                header_row_index = i
                break

    try:
        df = pd.read_csv(file_path, skiprows=header_row_index, nrows=0)
        file_headers = set(df.columns.str.strip())
    except Exception as e:
        return "unknown"
    for file_type, required_headers in known_formats.items():
        required_set = set(required_headers)
        if required_set.issubset(file_headers):
            return file_type
    return "unknown"


def validate_protein_sequence(sequence: str) -> bool:
    allowed_chars = "ACDEFGHIKLMNPQRSTVWY"

    for char in sequence:
        if not char.isalpha():
            return False
        if char.upper() not in allowed_chars:
            return False
    return True

def combine_hxms_data(hxms_data1:HxmsData,ptm_1_dic:dict,hxms_data2:HxmsData,ptm_2_dic:dict):
    largest_replicate = 0
    out_hxms = HxmsData()
    out_hxms.metadata = hxms_data1.metadata
    protein = hxms_data1.proteins[0]
    state = hxms_data1.state[protein][0]
    peptides = hxms_data1.peptides[protein][state]
    if protein not in out_hxms.proteins:
        out_hxms.proteins.append(protein)
        out_hxms.state[protein] = []
        out_hxms.peptides[protein] = {}
    if state not in out_hxms.state[protein]:
        out_hxms.state[protein].append(state)
        out_hxms.peptides[protein][state] = {}


    for peptide in peptides:
        if peptide not in out_hxms.peptides[protein][state]:
            out_hxms.peptides[protein][state][peptide] = {}
        for tp in peptides[peptide]:
            if tp not in out_hxms.peptides[protein][state][peptide]:
                out_hxms.peptides[protein][state][peptide][tp] = []
            for measurment in peptides[peptide][tp]:
                try:
                    rep_value = int(measurment["replicate"])
                    if rep_value > largest_replicate:
                        largest_replicate = rep_value
                    out_hxms.peptides[protein][state][peptide][tp].append(measurment)
                except:
                    return
    largest_replicate += 1
    peptides = hxms_data2.peptides[protein][state]
    for peptide in peptides:
        if peptide not in out_hxms.peptides[protein][state]:
            out_hxms.peptides[protein][state][peptide] = {}
        for tp in peptides[peptide]:
            if tp not in out_hxms.peptides[protein][state][peptide]:
                out_hxms.peptides[protein][state][peptide][tp] = []
            for measurment in peptides[peptide][tp]:
                try:
                    measurment["replicate"] = measurment["replicate"] + largest_replicate
                    out_hxms.peptides[protein][state][peptide][tp].append(measurment)
                except:
                    return

    out_hxms = _sort_hxms_data_cc(out_hxms)
    return out_hxms