Merge pull request #57 from mwang87/cdf-gc-support

[CDF Support] CDF support from GC and single quad instruments

Author: mwang87

app.py

@@ -106,7 +106,7 @@
         ),
         dbc.Nav(
             [
-                dbc.NavItem(dbc.NavLink("GNPS LCMS Dashboard - Version 0.14", href="/")),
+                dbc.NavItem(dbc.NavLink("GNPS LCMS Dashboard - Version 0.15", href="/")),
             ],
         navbar=True)
     ],
@@ -726,6 +726,8 @@
             html.Br(),
             html.A("Thermo GCMS", href="/?usi=mzspec:MSV000086150:BA1.mzML"),
             html.Br(),
+            html.A("GCMS in CDF Format", href="/?usi=mzspec:GNPS:TASK-ce31b7fdd01244dbb31478147889de1e-f.aaksenov/GC_data/Sterols_data/Samples/0104006.cdf"),
+            html.Br(),
             html.A("Thermo LCMS from GNPS Analysis Classical Molecular Networking Task", href="/?usi=mzspec:GNPS:TASK-5ecfcf81cb3c471698995b194d8246a0-f.MSV000085444/ccms_peak/peak/Hui_N1_fe.mzML#%7B%7D"),
             html.Br(),
         ]
@@ -1147,8 +1149,11 @@ def update_output(search, filecontent, filename, filedate):
     usi2 = ""
 
     if filecontent is not None:
+        if len(filecontent) > 100000000:
+            raise Exception
+
         extension = os.path.splitext(filename)[1]
-        if extension == ".mzML" and len(filecontent) < 100000000:
+        if extension == ".mzML":
             temp_filename = os.path.join("temp", "{}.mzML".format(str(uuid.uuid4())))
             data = filecontent.encode("utf8").split(b";base64,")[1]
 
@@ -1159,7 +1164,7 @@ def update_output(search, filecontent, filename, filedate):
 
             return [usi, usi2, "FILE Uploaded {}".format(filename)]
 
-        if extension == ".mzXML" and len(filecontent) < 100000000:
+        if extension == ".mzXML":
             mangled_name = str(uuid.uuid4())
             temp_filename_mzXML = os.path.join("temp", "{}.mzXML".format(mangled_name))
             temp_filename = os.path.join("temp", "{}.mzXML".format(mangled_name))
@@ -1172,6 +1177,19 @@ def update_output(search, filecontent, filename, filedate):
 
             return [usi, usi2, "FILE Uploaded {}".format(filename)]
 
+        if extension.lower() == ".cdf":
+            mangled_name = str(uuid.uuid4())
+            temp_filename = os.path.join("temp", "{}.cdf".format(mangled_name))
+            data = filecontent.encode("utf8").split(b";base64,")[1]
+
+            with open(temp_filename, "wb") as temp_file:
+                temp_file.write(base64.decodebytes(data))
+
+            usi = "mzspec:LOCAL:{}".format(os.path.basename(temp_filename))
+
+            return [usi, usi2, "FILE Uploaded {}".format(filename)]
+
+
 
     # Resolving USI
     usi = _get_param_from_url(search, "usi", usi)

lcms_map.py

@@ -1,4 +1,5 @@
 
+import os
 import pymzml
 import numpy as np
 import datashader as ds
@@ -11,7 +12,6 @@
 import plotly.graph_objects as go 
 
 
-
 def _gather_lcms_data(filename, min_rt, max_rt, min_mz, max_mz, polarity_filter="None"):
     all_mz = []
     all_rt = []

requirements.txt

@@ -43,4 +43,5 @@ Flask-Caching
 scipy
 pyteomics
 psims
-molmass
+molmass
+netcdf4

test/bin

@@ -0,0 +1 @@
+../bin/

test/test.py

@@ -44,4 +44,3 @@ def test_resolve():
         print(record["usi"])
         remote_link, local_filename = utils._resolve_usi(record["usi"])
         lcms_map._create_map_fig(local_filename)
-

test/usi_list.tsv

@@ -15,4 +15,5 @@ mzspec:GNPS:GNPS-EMBL-MCF:accession:CCMSLIB00000479565
 mzspec:MTBLS138:QXA01PDNEG20181113_KGCO0418VWSS7000000126_RAT_SERUM1_03.mzML
 mzspec:MSV000083500:ccms_peak/9177.mzML
 mzspec:ST000763:20160411_MB_CS00000074-1_P.mzXML
-mzspec:MSV000086456:peak/Pilot_plate_5_NEG_mzXML/PLT5_B1_BB1_01_60017.mzXML
+mzspec:MSV000086456:peak/Pilot_plate_5_NEG_mzXML/PLT5_B1_BB1_01_60017.mzXML
+mzspec:GNPS:TASK-ce31b7fdd01244dbb31478147889de1e-f.aaksenov/GC_data/Sterols_data/Samples/0104006.cdf

utils.py

@@ -29,7 +29,10 @@ def _resolve_usi(usi, temp_folder="temp"):
         if not os.path.isfile(converted_local_filename):
             temp_filename = os.path.join(temp_folder, str(uuid.uuid4()) + ".mzML")
             # Lets do a conversion
-            _convert_mzML(local_filename, temp_filename)
+            if file_extension == ".cdf":
+                _convert_cdf_to_mzML(local_filename, temp_filename)
+            else:
+                _convert_mzML(local_filename, temp_filename)
 
             os.rename(temp_filename, converted_local_filename)
 
@@ -120,8 +123,9 @@ def _resolve_usi(usi, temp_folder="temp"):
     # Getting Data Local, TODO: likely should serialize it
     local_filename = os.path.join(temp_folder, werkzeug.utils.secure_filename(remote_link))
     filename, file_extension = os.path.splitext(local_filename)
+
     converted_local_filename = filename + ".mzML"
-    
+
     if not os.path.isfile(converted_local_filename):
         temp_filename = os.path.join(temp_folder, str(uuid.uuid4()) + file_extension)
         wget_cmd = "wget '{}' -O {}".format(remote_link, temp_filename)
@@ -130,12 +134,14 @@ def _resolve_usi(usi, temp_folder="temp"):
 
         temp_filename = os.path.join(temp_folder, str(uuid.uuid4()) + ".mzML")
         # Lets do a conversion
-        _convert_mzML(local_filename, temp_filename)
+        if file_extension == ".cdf":
+            _convert_cdf_to_mzML(local_filename, temp_filename)
+        else:
+            _convert_mzML(local_filename, temp_filename)
 
+        # Renaming the temp
         os.rename(temp_filename, converted_local_filename)
 
-        local_filename = converted_local_filename
-
     return remote_link, converted_local_filename
 
 # First try msconvert, if the output fails, then we will do pyteomics to mzML and then msconvert
@@ -220,8 +226,72 @@ def _convert_mzML(input_mzXML, output_mzML):
         except:
             pass
 
+# in python doing a conversion from cdf to mzML
+def _convert_cdf_to_mzML(input_cdf, output_mzML):
+    from netCDF4 import Dataset
+    from psims.mzml.writer import MzMLWriter
+    import numpy as np
+
+    temp_filename = os.path.join("temp", str(uuid.uuid4()) + ".mzML")
+
+    # lets put the cdf reader here
+    cdf_reader = Dataset(input_cdf, "r")
+    mass_values      = np.array(cdf_reader.variables["mass_values"][:])
+    intensity_values = np.array(cdf_reader.variables["intensity_values"][:])
+    time_values = np.array(cdf_reader.variables["scan_acquisition_time"][:])
+    scan_values = np.array(cdf_reader.variables["scan_index"][:])
+
+    #removing empty scans
+    dd = np.diff(scan_values) != 0
+    dd = np.append(dd, True)
+    ddi = np.arange(scan_values.shape[0], dtype = np.int64)[dd]
+    time_values = time_values[ddi]
+    scan_values = scan_values[ddi]
+
+    # getting scan boundaries
+    scan_end_values = np.append(scan_values[1:]-1, mass_values.shape[0]-1)
+    scan_indcs = zip(scan_values, scan_end_values)
+
+    # Writing everything out
+    with MzMLWriter(open(temp_filename, 'wb')) as out:
+        out.controlled_vocabularies()
+        with out.run(id="my_analysis"):
+            with out.spectrum_list(count=1000):
+                # Iterating through all scans in the reader
+                try:
+                    # reading through scans
+                    for i, scan_range in enumerate(scan_indcs):
+                        time_min_rt = time_values[i] / 60
+
+                        _mz_array = np.array(mass_values[scan_range[0]:scan_range[1]])
+                        _i_array = np.array(intensity_values[scan_range[0]:scan_range[1]])
+                        
+                        out.write_spectrum(
+                            _mz_array, _i_array,
+                            id=i, params=[
+                                "MS1 Spectrum",
+                                {"ms level": 1},
+                                {"total ion current": sum(_i_array)}
+                            ],
+                            scan_start_time=time_min_rt)
+                except:
+                    print("Reading Failed, skipping to end")
+                    pass
+
+    # # Round trip through MsConvert
+    # conversion_cmd = "export LC_ALL=C && ./bin/msconvert {} --mzML --32 --outfile {} --outdir {} --filter 'threshold count 500 most-intense'".format(temp_filename, output_mzML, os.path.dirname(output_mzML))
+    # conversion_ret_code = os.system(conversion_cmd)
+
+    # try:
+    #     os.remove(temp_filename)
+    # except:
+    #     pass
+
+    try:
+        os.rename(temp_filename, output_mzML)
+    except:
+        pass
 
-    
 
 import subprocess, io