Thermo Stellar SRM spectra improvements

pwiz/data/vendor_readers/Thermo/ChromatogramList_Thermo.cpp

@@ -26,6 +26,8 @@
 
 #include "ChromatogramList_Thermo.hpp"
 
+#include <boost/range/algorithm/find.hpp>
+
 
 #ifdef PWIZ_READER_THERMO
 #include "Reader_Thermo_Detail.hpp"
@@ -193,10 +195,9 @@ PWIZ_API_DECL ChromatogramPtr ChromatogramList_Thermo::chromatogram(size_t index
 
                 string q1 = (format("%.10g", std::locale::classic()) % ci.q1).str();
                 string q3Range = (format("%.10g-%.10g", std::locale::classic())
-                                  % (ci.q3 - ci.q3Offset)
-                                  % (ci.q3 + ci.q3Offset)
+                                  % (ci.filterQ3 - ci.filterOffset) // offset must match filter
+                                  % (ci.filterQ3 + ci.filterOffset)
                                  ).str();
-
                 ChromatogramDataPtr cd = rawfile_->getChromatogramData(Type_MassRange,
                     polarity + "SRM ms2 " + q1 + " [" + q3Range + "]", ci.q3 - ci.q3Offset, ci.q3 + ci.q3Offset, 0,
                     rawfile_->getFirstScanTime(), rawfile_->getLastScanTime());
@@ -330,9 +331,45 @@ PWIZ_API_DECL ChromatogramList_Thermo::IndexEntry& ChromatogramList_Thermo::addC
     return ci;
 }
 
+template<typename MapT>
+typename MapT::const_iterator find_nearest(MapT const& m, typename MapT::key_type const& query, typename MapT::key_type const& tolerance)
+{
+    typename MapT::const_iterator cur, min, max, best;
+
+    min = m.lower_bound(query - tolerance);
+    max = m.lower_bound(query + tolerance);
+
+    if (min == m.end() || fabs(query - min->first) > tolerance)
+        return m.end();
+    else if (min == max)
+        return min;
+    else
+        best = min;
+
+    double minDiff = fabs(query - best->first);
+    for (cur = min; cur != max; ++cur)
+    {
+        double curDiff = fabs(query - cur->first);
+        if (curDiff < minDiff)
+        {
+            minDiff = curDiff;
+            best = cur;
+        }
+    }
+    return best;
+}
 
 PWIZ_API_DECL void ChromatogramList_Thermo::createIndex() const
 {
+    /*MassListTable massListTable;
+    for(const auto& method : rawfile_->getInstrumentMethods())
+    {
+        massListTable = parseMassListTables(method);
+        if (!massListTable.empty())
+            break;
+    }*/
+
+
     for (int controllerType = Controller_MS;
          controllerType < Controller_Count;
          ++controllerType)
@@ -390,22 +427,53 @@ PWIZ_API_DECL void ChromatogramList_Thermo::createIndex() const
                                 ci.q1 = filterParser.precursorMZs_[0];
                                 idMap_[ci.id] = ci.index;*/
 
+                                double q1 = scanInfo->precursorMZ(0);
+
                                 for (size_t j=0, jc=scanInfo->scanRangeCount(); j < jc; ++j)
                                 {
-                                    double q1 = scanInfo->precursorMZ(0);
-                                    double q3 = (scanInfo->scanRange(j).first + scanInfo->scanRange(j).second) / 2.0;
+                                    double scanRange = scanInfo->scanRange(j).second - scanInfo->scanRange(j).first;
+                                    double filterQ3 = (scanInfo->scanRange(j).first + scanInfo->scanRange(j).second) / 2.0;
                                     auto polarityType = translate(scanInfo->polarityType());
                                     string polarity = polarityStringForFilter(polarityType);
+
+                                    // if there's a mass list table entry for this q1/q3, use it instead of the scan range
+                                    /*auto findPrecursorItr = find_nearest(massListTable, q1, 0.001);
+                                    if (findPrecursorItr != massListTable.end())
+                                    {
+                                        auto findProductStartItr = findPrecursorItr->second.lower_bound(scanInfo->scanRange(j).first);
+                                        auto findProductEndItr = findPrecursorItr->second.lower_bound(scanInfo->scanRange(j).second);
+                                        for(const auto& itr : boost::iterator_range<MassListTable::mapped_type::const_iterator>(findProductStartItr, findProductEndItr))
+                                        {
+                                            double q3 = itr.second.product_mz;
+                                            string id = (format("%sSRM SIC %s,%.10g", std::locale::classic())
+                                                % polarity
+                                                % precursorMZ
+                                                % q3
+                                                ).str();
+                                            IndexEntry& ci = addChromatogram(id, (ControllerType)controllerType, n, MS_SRM_chromatogram, filterString);
+                                            ci.q1 = q1;
+                                            ci.q3 = q3;
+                                            ci.polarityType = polarityType;
+                                            ci.q3Offset = 0.35;
+                                            ci.filterQ3 = filterQ3;
+                                            ci.filterOffset = scanRange / 2.0;
+                                        }
+                                        continue;
+                                    }*/
+
+                                    if (scanRange > MAX_SRM_SCAN_RANGE)
+                                        continue; // skip this one, it's too wide
+
                                     string id = (format("%sSRM SIC %s,%.10g", std::locale::classic())
                                              % polarity
                                              % precursorMZ
-                                             % q3
+                                             % filterQ3
                                             ).str();
                                     IndexEntry& ci = addChromatogram(id, (ControllerType)controllerType, n, MS_SRM_chromatogram, filterString);
                                     ci.q1 = q1;
-                                    ci.q3 = q3;
+                                    ci.q3 = ci.filterQ3 = filterQ3;
                                     ci.polarityType = polarityType;
-                                    ci.q3Offset = (scanInfo->scanRange(j).second - scanInfo->scanRange(j).first) / 2.0;
+                                    ci.filterOffset = ci.q3Offset = scanRange / 2.0;
                                 }
                             }
                             break; // case ScanType_SRM

pwiz/data/vendor_readers/Thermo/ChromatogramList_Thermo.hpp

@@ -81,6 +81,7 @@ class PWIZ_API_DECL ChromatogramList_Thermo : public ChromatogramListBase
         string filter;
         double q1, q3;
         double q3Offset;
+        double filterQ3, filterOffset; // SRM Q3 and its range must match scan filter, even if we use mass table to get actual Q3
         CVID polarityType;
     };
 

pwiz/data/vendor_readers/Thermo/Jamfile.jam

@@ -119,7 +119,8 @@ lib pwiz_reader_thermo
         <conditional>@vendor-api-usage-requirements
     ;
 
-doctest SpectrumList_ThermoTest : SpectrumList_Thermo.cpp : <library>$(PWIZ_ROOT_PATH)/pwiz/data/msdata//pwiz_data_msdata <conditional>@vendor-api-requirements ;
+doctest SpectrumList_ThermoTest : SpectrumList_Thermo.cpp : <library>pwiz_reader_thermo ;
+doctest Reader_Thermo_DetailTest : Reader_Thermo_Detail.cpp : <library>pwiz_reader_thermo ;
 
 
 rule warn-once ( message )

pwiz/data/vendor_readers/Thermo/Reader_Thermo_Detail.cpp

@@ -27,6 +27,7 @@
 #include "pwiz/utility/misc/Container.hpp"
 #include "pwiz/utility/misc/String.hpp"
 #include <boost/range/algorithm/find_if.hpp>
+#include "pwiz/utility/misc/unit.hpp"
 
 namespace pwiz {
 namespace msdata {
@@ -584,6 +585,126 @@ PWIZ_API_DECL void setActivationType(ActivationType activationType, ActivationTy
     // ActivationType_Unknown:
 }
 
+
+MassListTable parseMassListTables(const std::string& instrumentMethod)
+{
+    MassListTable entries;
+    const std::string startToken = ">>>>>>>>>>>>> Mass List Table <<<<<<<<<<<<<<";
+    const std::string endToken = ">>>>>>>>>>>>> End Mass List Table <<<<<<<<<<<<<<";
+    size_t pos = 0;
+
+    while ((pos = instrumentMethod.find(startToken, pos)) != std::string::npos)
+    {
+        size_t tableStart = pos + startToken.size();
+        size_t tableEnd = instrumentMethod.find(endToken, tableStart);
+        if (tableEnd == std::string::npos) break;
+
+        std::istringstream iss(instrumentMethod.substr(tableStart, tableEnd - tableStart));
+        std::string line;
+        std::vector<std::string> headers;
+        std::vector<size_t> colIdx(5, std::string::npos);
+
+        // Find header line
+        while (std::getline(iss, line))
+        {
+            bal::trim(line);
+            if (line.empty()) continue;
+            if (line.find('|') != std::string::npos)
+            {
+                std::istringstream headerStream(line);
+                std::string header;
+                size_t idx = 0;
+                while (std::getline(headerStream, header, '|'))
+                {
+                    bal::trim(header);
+                    headers.push_back(header);
+                    if (header == "Compound") colIdx[0] = idx;
+                    else if (header == "Precursor (m/z)") colIdx[1] = idx;
+                    else if (header == "Product (m/z)") colIdx[2] = idx;
+                    else if (header == "Collision Energies (V)") colIdx[3] = idx;
+                    else if (header == "Polarity") colIdx[4] = idx;
+                    ++idx;
+                }
+                break;
+            }
+        }
+
+        // Only parse tables with all required columns
+        if (std::all_of(colIdx.begin(), colIdx.end(), [](size_t i) { return i != std::string::npos; }))
+        {
+            // Parse rows
+            while (std::getline(iss, line))
+            {
+                bal::trim(line);
+                if (line.empty() || line.find('|') == std::string::npos) continue;
+
+                std::istringstream rowStream(line);
+                std::string cell;
+                std::vector<std::string> cells;
+                while (std::getline(rowStream, cell, '|'))
+                    cells.push_back(bal::trim_copy(cell));
+
+                if (cells.size() < headers.size()) continue;
+
+                MassListTableEntry entry;
+                entry.compound = cells[colIdx[0]];
+                entry.precursor_mz = lexical_cast<double>(cells[colIdx[1]]);
+                entry.product_mz = lexical_cast<double>(cells[colIdx[2]]);
+                entry.collision_energy = lexical_cast<double>(cells[colIdx[3]]);
+                entry.polarity = cells[colIdx[4]];
+                entries[entry.precursor_mz][entry.product_mz] = entry;
+            }
+        }
+        pos = tableEnd + endToken.size();
+    }
+    return entries;
+}
+
+namespace {
+
+TEST_CASE("parseMassListTables()") {
+    MassListTable result = parseMassListTables(
+        R"(Stellar Method Summary	
+
+Creator: THERMO-RJM0O65J\Thermo Scientific          Last Modified: 09/11/2025 11:26:26 by THERMO-RJM0O65J\Thermo Scientific	
+
+Global Settings	
+	Use Ion Source Settings from Tune = Not Checked	
+	Method Duration (min) = 15	
+			Use Chromatographic Filter = False	
+>>>>>>>>>>>>> Mass List Table <<<<<<<<<<<<<<	
+                Center Mass (m/z)|           Isolation Window (m/z)|	
+                              125|                               50|	
+                              175|                               50|	
+                              225|                               50|	
+                              275|                               50|	
+                              325|                               50|	
+                              375|                               50|	
+>>>>>>>>>>>>> End Mass List Table <<<<<<<<<<<<<<	
+Experiment 2	
+>>>>>>>>>>>>> Mass List Table <<<<<<<<<<<<<<	
+Compound|  Formula|  Adduct|  Precursor (m/z)|  Precursor z|  Product (m/z)|  Retention Time (min)|  RT Window (min)|  Collision Energies (V)|  Polarity|	
+AcylCarnitine 14:2|  |  (no adduct)|  368.3|  1|  85.1|  1.4|  0.4|  27|  Positive|	
+AcylCarnitine 12:0|  |  (no adduct) | 344.3 | 1 | 84.1 | 1.5 | 0.4 | 24 | Negative |
+>>>>>>>>>>>>> End Mass List Table <<<<<<<<<<<<<<
+)");
+
+    CHECK(result.size() == 2);
+    CHECK(result[368.3][85.1].compound == "AcylCarnitine 14:2");
+    CHECK(result[368.3][85.1].precursor_mz == 368.3);
+    CHECK(result[368.3][85.1].product_mz == 85.1);
+    CHECK(result[368.3][85.1].collision_energy == 27);
+    CHECK(result[368.3][85.1].polarity == "Positive");
+
+    CHECK(result[344.3][84.1].compound == "AcylCarnitine 12:0");
+    CHECK(result[344.3][84.1].precursor_mz == 344.3);
+    CHECK(result[344.3][84.1].product_mz == 84.1);
+    CHECK(result[344.3][84.1].collision_energy == 24);
+    CHECK(result[344.3][84.1].polarity == "Negative");
+}
+
+} // namespace
+
 } // Thermo
 } // detail
 } // msdata

pwiz/data/vendor_readers/Thermo/Reader_Thermo_Detail.hpp

@@ -54,6 +54,23 @@ PWIZ_API_DECL CVID translateAsInletType(IonizationType ionizationType);
 PWIZ_API_DECL CVID translate(PolarityType polarityType);
 PWIZ_API_DECL void setActivationType(ActivationType activationType, ActivationType supplementalActivationType, Activation& activation);
 
+constexpr double MAX_SRM_SCAN_RANGE = 1.0; // if Q3 m/z range is larger than this, don't produce SRM chromatogram
+
+
+PWIZ_API_DECL struct MassListTableEntry
+{
+    std::string compound;
+    double precursor_mz;
+    double product_mz;
+    double collision_energy;
+    std::string polarity;
+};
+
+typedef std::map<double, std::map<double, MassListTableEntry>> MassListTable; /// keyed by precursor m/z, then product m/z
+
+PWIZ_API_DECL MassListTable parseMassListTables(const std::string& instrumentMethod);
+
+
 } // Thermo
 } // detail
 } // msdata

pwiz/data/vendor_readers/Thermo/SpectrumList_Thermo.cpp

@@ -823,6 +823,8 @@ PWIZ_API_DECL void SpectrumList_Thermo::createIndex()
     spectraByScanType.resize(ScanType_Count, 0);
     spectraByMSOrder.resize(MSOrder_Count+3, 0); // can't use negative index and a std::map would be inefficient
 
+    auto raw = rawfile_->getRawByThread(0);
+
     // calculate total spectra count from all controllers
     for (int controllerType = Controller_MS;
          controllerType < Controller_Count;
@@ -877,9 +879,28 @@ PWIZ_API_DECL void SpectrumList_Thermo::createIndex()
                                     break;  // break out of switch (scanType)
                                 continue;
                             case ScanType_SRM:
+                            {
                                 if (config_.srmAsSpectra)
                                     break;  // break out of switch (scanType)
+
+                                // if any scan range is bigger than MAX_SRM_SCAN_RANGE, we can't skip the SRM scan
+                                const auto scanInfo = raw->getScanInfo(scan);
+                                bool hasExcessiveSrmScanRange = false;
+                                for (size_t i = 0, end = scanInfo->scanRangeCount(); i < end; ++i)
+                                {
+                                    const double scanRange = scanInfo->scanRange(i).second - scanInfo->scanRange(i).first;
+                                    if (scanRange > MAX_SRM_SCAN_RANGE)
+                                    {
+                                        hasExcessiveSrmScanRange = true;
+                                        break;
+                                    }
+
+                                }
+                                if (hasExcessiveSrmScanRange)
+                                    break;  // break out of switch (scanType)
+
                                 continue;
+                            }
                         }
 
                         index_.push_back(IndexEntry());

pwiz_tools/Shared/ProteowizardWrapper/MsDataFileImpl.cs

@@ -745,7 +745,7 @@ protected virtual SpectrumList SpectrumList
                 {
                     var centroidLevel = new List<int>();
                     _spectrumList = _msDataFile.run.spectrumList;
-                    bool hasSrmSpectra = HasSrmSpectraInList(_spectrumList);
+                    bool hasSrmSpectra = HasSrmSpectraInList();
                     if (!hasSrmSpectra)
                     {
                         if (_requireVendorCentroidedMS1)
@@ -1164,10 +1164,14 @@ private double[] GetIonMobilityArray(Spectrum s)
                             if (data == null)
                             {
                                 data = TryGetIonMobilityData(s, CVID.MS_mean_ion_mobility_drift_time_array, ref _cvidIonMobility);
-                                if (data == null && HasCombinedIonMobilitySpectra && !s.id.Contains(MERGED_TAG))
+                                if (data == null)
                                 {
-                                    _cvidIonMobility = null; // We can't learn anything from a lockmass spectrum that has no IMS
-                                    return null;
+                                    data = TryGetIonMobilityData(s, CVID.MS_raw_ion_mobility_drift_time_array, ref _cvidIonMobility);
+                                    if (data == null && HasCombinedIonMobilitySpectra && !s.id.Contains(MERGED_TAG))
+                                    {
+                                        _cvidIonMobility = null; // We can't learn anything from a lockmass spectrum that has no IMS
+                                        return null;
+                                    }
                                 }
                             }
                         }
@@ -1402,7 +1406,7 @@ private SpectrumMetadata GetSpectrumMetadata(Spectrum spectrum)
 
         public bool HasSrmSpectra
         {
-            get { return HasSrmSpectraInList(SpectrumList); }
+            get { return HasSrmSpectraInList(); }
         }
 
         public bool HasIonMobilitySpectra
@@ -1427,13 +1431,16 @@ public bool HasChromatogramData
             }
         }
 
-        private static bool HasSrmSpectraInList(SpectrumList spectrumList)
+        private bool HasSrmSpectraInList()
         {
-            if (spectrumList == null || spectrumList.size() == 0)
+            if (_spectrumList == null || _spectrumList.size() == 0)
                 return false;
 
+            if (_msDataFile.fileDescription.fileContent.hasCVParam(CVID.MS_SRM_spectrum))
+                return true;
+
             // If the first spectrum is not SRM, the others will not be either
-            using (var spectrum = spectrumList.spectrum(0, false))
+            using (var spectrum = _spectrumList.spectrum(0, false))
             {
                 return IsSrmSpectrum(spectrum);
             }