Add fix for 'empty' scans after filtering

Author: RJMW

dimspy/experiment.py

@@ -175,25 +175,33 @@ def check_metadata(fn_tsv):
     return fm_dict
 
 
-def update_metadata_and_labels(peaklists, fl):
+def update_metadata_and_labels(peaklists, fl, pl_id=""):
 
     if not isinstance(peaklists[0], PeakList):
         raise IOError("PeakList object required")
 
-    for k in fl.keys():
-        for pl in peaklists:
-            if pl.ID not in fl[fl.keys()[0]]:
-                raise IOError("filelist and peaklist do not match {}".format(pl.ID))
+    if len(fl) == 0:
+        return peaklists
 
-            index = fl[fl.keys()[0]].index(pl.ID)
+    for pl in peaklists:
+
+        if pl_id == "":
+            pl_ID = pl_id
+        else:
+            pl_ID = pl.ID
+
+        if pl_ID not in fl[fl.keys()[0]]:
+            raise IOError("filelist and peaklist do not match {}".format(pl_ID))
+
+        index = fl[fl.keys()[0]].index(pl_ID)
+        for k in fl.keys():
             pl.metadata[k] = fl[k][index]
-            #pl.metadata["filelist"] = {k:fl[k][index] for k in fl.keys()}
 
-            for tag_name in ["replicate", "replicates", "batch", "injectionOrder", "classLabel"]:
-                if tag_name in fl.keys():
-                    if pl.tags.has_tag_type(tag_name):
-                        pl.tags.drop_tag_type(tag_name)
-                    pl.tags.add_tag(fl[tag_name][index], tag_name)
+        for tag_name in ["replicate", "replicates", "batch", "injectionOrder", "classLabel"]:
+            if tag_name in fl.keys():
+                if pl.tags.has_tag_type(tag_name):
+                    pl.tags.drop_tag_type(tag_name)
+                pl.tags.add_tag(fl[tag_name][index], tag_name)
 
     return peaklists
 

dimspy/tools.py

@@ -105,6 +105,9 @@ def process_scans(source, function_noise, snr_thres, ppm, min_fraction=None, rsd
             if len(pls_scans[h]) >= 1:
                 if sum(pl.shape[0] for pl in pls_scans[h]) == 0:
                     logging.warning("No scan data available for {}".format(h))
+                    d = dict((k, [pls_scans[h][0].metadata[k]]) for k in pls_scans[h][0].metadata)
+                    pls_avg.append(PeakList(ID=h, mz=[], intensity=[], **d))
+                    n_peaks, median_rsd = 0, "NA"
                 else:
                     pl_avg = average_replicate_scans(h, pls_scans[h], ppm, min_fraction, rsd_thres, "intensity", block_size, ncpus)
                     pls_avg.append(pl_avg)
@@ -115,21 +118,26 @@ def process_scans(source, function_noise, snr_thres, ppm, min_fraction=None, rsd
             if report is not None:
                 out.write("{}\t{}\t{}\t{}\t{}\n".format(os.path.basename(filenames[i]), h, nscans, n_peaks, median_rsd))
 
-        if len(pls_avg) == 0:
-            raise IOError("No peaks remaining after filtering. Remove file from Study (filelist).")
+        if sum(pl.shape[0] for pl in pls_avg) == 0:
+            logging.warning("No peaks remaining after filtering. Remove file from Study (filelist).")
 
         if not skip_stitching or len(pls_scans.keys()) == 1:
             pl = join_peaklists(os.path.basename(filenames[i]), pls_avg)
             pl = update_metadata_and_labels([pl], fl)
             pls.extend(pl)
+
+            if hasattr(pl[0], 'rsd'):
+                median_rsd = np.nanmedian(pl[0].rsd)
+            else:
+                median_rsd = "NA"
+
             if len(pls_scans.keys()) > 1 and report is not None:
-                out.write("{}\t{}\t{}\t{}\t{}\n".format(os.path.basename(filenames[i]), "SIM-Stitch", "NA", pl[0].shape[0], np.nanmedian(pl[0].rsd)))
+                out.write("{}\t{}\t{}\t{}\t{}\n".format(os.path.basename(filenames[i]), "SIM-Stitch", "NA", pl[0].shape[0], median_rsd))
         else:
             for pl in pls_avg:
-                pl = update_metadata_and_labels([pl], fl)
-                pl = join_peaklists("{}#{}".format(os.path.basename(filenames[i]), pl[0].metadata["header"][0]), pl)
-                pls.append(pl)
-
+                pl = join_peaklists("{}#{}".format(os.path.basename(filenames[i]), pl.metadata["header"][0]), [pl])
+                pl = update_metadata_and_labels([pl], fl, os.path.basename(filenames[i]))
+                pls.extend(pl)
     return pls
 
 
@@ -244,9 +252,10 @@ def replicate_filter(source, ppm, replicates, min_peaks, rsd_thres=None, filelis
 
         if sum([comb[-1] for comb in temp]) == 0.0:
             logging.warning("insufficient data available to calculate scores for {}".format(str([comb[0].ID for comb in temp])))
-
-        # sort the scores from high to low
-        temp.sort(key=operator.itemgetter(-1), reverse=True)
+            temp.sort(key=operator.itemgetter(1), reverse=True)
+        else:
+            # sort the scores from high to low
+            temp.sort(key=operator.itemgetter(-1), reverse=True)
         # select the replicate filtered peaklist that is ranked first
         pls_rep_filt.append(temp[0][0])