Merge pull request #51 from bigbio/dev

Update github actions, more examples added, remove filtering of features based on percentile tic.

Author: ypriverol

.github/workflows/python-package.yml

@@ -44,3 +44,15 @@ jobs:
     - name: Test commandline tool
       run: |
         quantmsutilsc --help
+    - name: Download test files.
+      run: |
+        wget https://ftp.pride.ebi.ac.uk/pub/databases/pride/resources/proteomes/quantms-ci-github/quantms-utils/TMT_Erwinia_1uLSike_Top10HCD_isol2_45stepped_60min_01.mzML
+        wget https://ftp.pride.ebi.ac.uk/pub/databases/pride/resources/proteomes/quantms-ci-github/quantms-utils/RD139_Narrow_UPS1_0_1fmol_inj1.mzML
+    - name: Test package online TMT
+      run: |
+        quantmsutilsc mzmlstats --ms_path TMT_Erwinia_1uLSike_Top10HCD_isol2_45stepped_60min_01.mzML --ms2_file --feature_detection
+    - name: Text package online DIA
+      run: |
+        quantmsutilsc mzmlstats --ms_path RD139_Narrow_UPS1_0_1fmol_inj1.mzML --ms2_file --feature_detection
+        
+        

quantmsutils/mzml/ms1_feature_finder.py

@@ -4,10 +4,13 @@
  previous algorithm by Andy Lim https://github.com/bmx8177/MS1Connect
  published in https://doi.org/10.1093/bioinformatics/btad058.
 
- We improved the original algorithm by using FeatureFinderMultiplexAlgorithm instead of FeatureFinder as originally
- implemented by Andy Lim. Also, we annotated additional features such as min and max retention time and mz values.
-
- This algorithm is used to detect MS1 features from mzML files and save them to parquet format.
+ We improved the original algorithm with the following ideas:
+  - Using FeatureFinderMultiplexAlgorithm instead of FeatureFinder as originally implemented by Andy Lim. This will provide
+    a more robust way to perform FeatureFinding.
+  - Remove the filtering of percentile TIC for features, we leave this step to future consuming tools of the data to perform
+    extra curation of the features based on percentile_tic, or quality of the feature, etc.
+  - We annotated additional features such as min and max retention time and mz values.
+  - This algorithm is used to detect MS1 features from mzML files and save them to parquet format.
 """
 
 import bisect
@@ -28,7 +31,7 @@ class MS1FeatureDetector:
     Class for detecting MS1 features from mzML files and saving to parquet format.
     """
 
-    def __init__(self, min_ptic: float = 0.05, max_ptic: float = 0.95, ms_level: int = 1):
+    def __init__(self, ms_level: int = 1):
         """
         Initialize the MS1 feature detector.
 
@@ -43,10 +46,7 @@ def __init__(self, min_ptic: float = 0.05, max_ptic: float = 0.95, ms_level: int
         """
         # Configure logging
 
-        self.min_ptic = min_ptic
-        self.max_ptic = max_ptic
         self.ms_level = ms_level
-
         # Initialize options for file loading
         self.options = oms.PeakFileOptions()
         self.options.setMSLevels([self.ms_level])
@@ -142,7 +142,11 @@ def _find_ptic_for_rt(rt: float, rt_list: List[float], ptic_list: List[float]) -
         return ptic_left + rt_frac * (ptic_right - ptic_left)
 
     def _extract_features(
-        self, features: oms.FeatureMap, rt_list: List[float], ptic_list: List[float], scans: List[str]
+        self,
+        features: oms.FeatureMap,
+        rt_list: List[float],
+        ptic_list: List[float],
+        scans: List[str],
     ) -> List[Dict[str, Any]]:
         """
         Extract feature information and filter by pTIC.
@@ -182,30 +186,23 @@ def _extract_features(
             select_scans = self._get_selected_scans(scans, rt_list, minRT, maxRT)
             num_scans = len(select_scans)
 
-
-
-            # Filter by pTIC
-            if self.min_ptic <= ptic <= self.max_ptic:
-                feature_list.append(
-                    {
-                        "feature_mz": mz,
-                        "feature_intensity": intensity,
-                        "feature_rt": rt,
-                        "feature_charge": charge,
-                        "feature_percentile_tic": ptic,
-                        "feature_quality": quality,
-                        "feature_id": feature.getUniqueId(),
-                        "feature_min_rt": minRT,
-                        "feature_min_mz": minMZ,
-                        "feature_max_rt": maxRT,
-                        "feature_max_mz": maxMZ,
-                        "feature_num_scans": num_scans,
-                        "feature_scans": select_scans
-
-                    }
-                )
-            else:
-                logger.debug(f"Skipping feature at RT {rt} due to pTIC {ptic}")
+            feature_list.append(
+                {
+                    "feature_mz": mz,
+                    "feature_intensity": intensity,
+                    "feature_rt": rt,
+                    "feature_charge": charge,
+                    "feature_percentile_tic": ptic,
+                    "feature_quality": quality,
+                    "feature_id": feature.getUniqueId(),
+                    "feature_min_rt": minRT,
+                    "feature_min_mz": minMZ,
+                    "feature_max_rt": maxRT,
+                    "feature_max_mz": maxMZ,
+                    "feature_num_scans": num_scans,
+                    "feature_scans": select_scans,
+                }
+            )
 
         return feature_list
 
@@ -321,4 +318,3 @@ def _get_selected_scans(scans: List[str], rt_list: List[float], min_rt: float, m
             if min_rt <= rt <= max_rt:
                 selected_scans.append(scans[i])
         return selected_scans
-

quantmsutils/mzml/mzml_statistics.py

@@ -27,7 +27,7 @@
     SCAN,
     SUMMED_PEAK_INTENSITY,
     PRECURSOR_RT,
-    PRECURSOR_TOTAL_INTENSITY
+    PRECURSOR_TOTAL_INTENSITY,
 )
 
 logging.basicConfig(format="%(asctime)s [%(funcName)s] - %(message)s", level=logging.INFO)
@@ -607,7 +607,7 @@ def mzml_statistics(
                 batch_size=batch_size,
             )
             if feature_detection:
-                feature_detector = MS1FeatureDetector(min_ptic=0.05, max_ptic=0.95)
+                feature_detector = MS1FeatureDetector()
                 feature_detector.process_file(input_file=ms_path, output_file=feature_output_path)
             logger.info("The file {} has been processed".format(ms_path))
         else:

quantmsutils/openms.py

@@ -21,4 +21,4 @@ def extract_scan_id(spectrum: oms.MSSpectrum) -> str:
     match = re.search(SCAN_PATTERN, spectrum.getNativeID())
     if match:
         return match.group(1)
-    return spectrum.getNativeID()
+    return spectrum.getNativeID()

tests/test_commands.py

@@ -204,7 +204,7 @@ class TestFeatureFinder:
     def test_feature_finder(self):
         """Test feature finder on TMT data"""
 
-        detector = MS1FeatureDetector(min_ptic=0.05, max_ptic=0.95)
+        detector = MS1FeatureDetector()
         result = detector.process_file(input_file=TMT_MZML_FILE, output_file=TMT_MS1_FEAURES)
 
         parquet_df = pd.read_parquet(result)