Merge branch '151_coremshdf5importer2' into 'master'

Fix error on loading CoreMS hdf5 molecular formula embedded in mass spectrum objects

See merge request mass-spectrometry/corems!120

Author: corilo

corems/mass_spectrum/input/coremsHDF5.py

@@ -199,6 +199,11 @@ def get_dataframe(self, scan_index=0, time_index=-1):
             data_dict = {}
             for data_index, data in enumerate(row):
                 label = columnsLabels[data_index]
+                # if data starts with a b' it is a byte string, so decode it
+                if isinstance(data, bytes):
+                    data = data.decode("utf-8")
+                if data == "nan":
+                    data = None
                 data_dict[label] = data
 
             list_dict.append(data_dict)

corems/mass_spectrum/input/massList.py

@@ -1,9 +1,12 @@
 __author__ = "Yuri E. Corilo"
 __date__ = "Jun 12, 2019"
 
+import numpy as np
+
+from corems.encapsulation.constant import Atoms
 from corems.mass_spectrum.input.baseClass import MassListBaseClass
 from corems.mass_spectrum.factory.MassSpectrumClasses import MassSpecProfile, MassSpecCentroid
-from corems.molecular_formula.factory.MolecularFormulaFactory import MolecularFormula
+from corems.molecular_formula.factory.MolecularFormulaFactory import MolecularFormula, MolecularFormulaIsotopologue
 from corems.encapsulation.constant import Labels, Atoms
 from corems.encapsulation.factory.processingSetting  import DataInputSetting
 
@@ -101,17 +104,81 @@ def add_molecular_formula(self, mass_spec_obj, dataframe):
                 atoms = list(formula_df.columns.astype(str))
                 counts = list(formula_df.iloc[df_index].astype(int))
 
-                formula_list = [sub[item] for item in range(len(atoms))
-                                for sub in [atoms, counts]]
+                formula_dict = dict(zip(atoms, counts))
             if sum(counts) > 0:
 
                 ion_type = str(Labels.ion_type_translate.get(ion_type_df[df_index]))
                 if adduct_df is not None:
                     adduct_atom = str(adduct_df[df_index])
+                    if adduct_atom == 'None':
+                        adduct_atom = None
                 else:
                     adduct_atom = None
-                mfobj = MolecularFormula(formula_list, int(ion_charge_df[df_index]), mspeak_parent=mass_spec_obj[ms_peak_index] , ion_type=ion_type, adduct_atom=adduct_atom)
-                mfobj.is_isotopologue = bool(is_isotopologue_df[df_index])
+
+                # If not isotopologue, cast as MolecularFormula
+                if not bool(int(is_isotopologue_df[df_index])):
+                    mfobj = MolecularFormula(
+                        formula_dict, int(ion_charge_df[df_index]), 
+                        mspeak_parent=mass_spec_obj[ms_peak_index] , 
+                        ion_type=ion_type, adduct_atom=adduct_atom
+                        )
+                    
+                # if is isotopologue, recast as MolecularFormulaIsotopologue
+                if bool(int(is_isotopologue_df[df_index])):
+
+                    # First make a MolecularFormula object for the parent so we can get probabilities etc
+                    formula_list_parent = {}
+                    for atom in formula_dict:
+                        if atom in Atoms.isotopes.keys():
+                            formula_list_parent[atom] = formula_dict[atom]
+                        else:
+                            # remove any numbers from the atom name to cast as a mono-isotopic atom
+                            atom_mono = atom.strip('0123456789')
+                            if atom_mono in Atoms.isotopes.keys():
+                                formula_list_parent[atom_mono] = formula_list_parent[atom_mono]+formula_dict[atom]
+                            else:
+                                print(f"Atom {atom} not in Atoms.atoms_order")
+                    mono_index = int(dataframe.iloc[df_index]['Mono Isotopic Index'])
+                    mono_mfobj = MolecularFormula(
+                        formula_list_parent, 
+                        int(ion_charge_df[df_index]), 
+                        mspeak_parent=mass_spec_obj[mono_index], 
+                        ion_type=ion_type, 
+                        adduct_atom=adduct_atom
+                        )
+                    
+                    # Next, generate isotopologues from the parent
+                    isos = list(
+                        mono_mfobj.isotopologues(
+                        min_abundance = mass_spec_obj[df_index].abundance*0.1, 
+                        current_mono_abundance = mass_spec_obj[mono_index].abundance, 
+                        dynamic_range = mass_spec_obj.dynamic_range
+                         )
+                    )
+
+                    # Finally, find the isotopologue that matches the formula_dict
+                    matched_isos = isos
+                    for iso in isos:
+                        if set(iso.atoms) == set(formula_dict.keys()):
+                            # Check the values of the atoms match
+                            if all([iso[atom] == formula_dict[atom] for atom in formula_dict]):
+                                matched_isos = [iso]
+                    if len(matched_isos) > 1:
+                        raise ValueError("More than one isotopologue matched the formula_dict: {matched_isos}")
+                    if len(matched_isos) == 0:
+                        raise ValueError("No isotopologue matched the formula_dict")
+                    mfobj = matched_isos[0]        
+
+                    # Add the mono isotopic index, confidence score and isotopologue similarity    
+                    mfobj.mspeak_index_mono_isotopic = int(dataframe.iloc[df_index]['Mono Isotopic Index'])
+                
+                # Add the confidence score and isotopologue similarity and average MZ error score
+                if 'm/z Error Score' in dataframe:
+                    mfobj._mass_error_average_score = float(dataframe.iloc[df_index]['m/z Error Score'])
+                if 'Confidence Score' in dataframe:
+                    mfobj._confidence_score = float(dataframe.iloc[df_index]['Confidence Score'])
+                if 'Isotopologue Similarity' in dataframe:
+                    mfobj._isotopologue_similarity = float(dataframe.iloc[df_index]['Isotopologue Similarity'])
                 mass_spec_obj[ms_peak_index].add_molecular_formula(mfobj)
 
 

tests/test_input.py

@@ -169,6 +169,7 @@ def test_import_corems_hdf5():
     mass_list_reader = ReadCoreMSHDF_MassSpectrum(file_location)
 
     mass_spectrum = mass_list_reader.get_mass_spectrum()
+    mass_spectrum.to_dataframe()
 
     for mspeak in mass_spectrum:
 
@@ -362,11 +363,11 @@ def test_import_thermo_average():
     #test_import_lcms_from_transient()
     #test_import_thermo_profile_mass_list()
     # test_import_transient()
-    #test_import_corems_hdf5()
+    test_import_corems_hdf5()
     #test_import_corems_mass_list()
     #test_import_mass_list()
     #test_import_maglab_pks()
     #test_andi_netcdf_gcms()
-    test_import_corems_mass_list()
+    #test_import_corems_mass_list()
     #test_import_thermo_average()
 

tests/test_mass_spectra_export_import.py

@@ -0,0 +1,74 @@
+import shutil
+
+from corems.mass_spectrum.input.numpyArray import ms_from_array_centroid
+from corems.molecular_id.search.molecularFormulaSearch import SearchMolecularFormulas
+from corems.encapsulation.factory.parameters import MSParameters
+from corems.mass_spectrum.output.export import HighResMassSpecExport
+from corems.mass_spectrum.input.coremsHDF5 import ReadCoreMSHDF_MassSpectrum
+
+
+def prep_mass_spec_obj():
+    # Test for generating accurate molecular formula from a single mass using the local sql database
+    # Now also tests that it is handling isotopes correctly (for non-adducts)
+    mz = [760.58156938877, 761.58548]
+    abundance = [1000, 400]
+    rp, s2n = [[1, 1], [10, 10]]
+
+    MSParameters.mass_spectrum.noise_threshold_method = "relative_abundance"
+    MSParameters.mass_spectrum.noise_threshold_absolute_abundance = 0
+
+    MSParameters.molecular_search.url_database = ""
+    MSParameters.molecular_search.error_method = "None"
+    MSParameters.molecular_search.min_ppm_error = -5
+    MSParameters.molecular_search.max_ppm_error = 5
+    MSParameters.molecular_search.mz_error_range = 1
+    MSParameters.molecular_search.isProtonated = True
+    MSParameters.molecular_search.isRadical = False
+    MSParameters.molecular_search.isAdduct = False
+
+    usedatoms = {"C": (1, 57), "H": (4, 200), "N": (0, 1)}
+    MSParameters.molecular_search.usedAtoms = usedatoms
+    mass_spectrum_obj = ms_from_array_centroid(
+        mz, abundance, rp, s2n, "single mf search", polarity=1, auto_process=True
+    )
+    return mass_spectrum_obj
+
+
+def run_molecular_formula_search(mass_spectrum_obj):
+    mass_spectrum_obj.molecular_search_settings.use_min_peaks_filter = False
+    mass_spectrum_obj.molecular_search_settings.use_isotopologue_filter = False
+    SearchMolecularFormulas(
+        mass_spectrum_obj, find_isotopologues=True
+    ).run_worker_ms_peaks([mass_spectrum_obj[0]])
+    return mass_spectrum_obj
+
+
+def test_mass_spec_export_import_with_annote():
+    mass_spectrum_obj = prep_mass_spec_obj()
+    mass_spectrum_obj = run_molecular_formula_search(mass_spectrum_obj)
+    ms_df1 = mass_spectrum_obj.to_dataframe()
+    assert mass_spectrum_obj[0][0].string == "C56 H73 N1"
+    assert ms_df1.shape == (2, 26)
+    assert mass_spectrum_obj[1][0].string == "C55 H73 N1 13C1"
+
+    exportMS = HighResMassSpecExport("my_mass_spec", mass_spectrum_obj)
+    exportMS._output_type = "hdf5"
+    exportMS.save()
+
+    parser = ReadCoreMSHDF_MassSpectrum("my_mass_spec.hdf5")
+    mass_spectrum_obj2 = parser.get_mass_spectrum(auto_process=True, load_settings=True)
+
+    ms_df2 = mass_spectrum_obj2.to_dataframe()
+    assert mass_spectrum_obj2[0][0].string == "C56 H73 N1"
+    assert ms_df2.shape == (2, 26)
+    assert mass_spectrum_obj2[1][0].string == "C55 H73 N1 13C1"
+
+    # Remove the file
+    shutil.rmtree(
+        "my_mass_spec.hdf5",
+        ignore_errors=True,
+    )
+
+
+if __name__ == "__main__":
+    test_mass_spec_export_import_with_annote()