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