MassIndexingEngine (smith-chem-wisc#886)

* Created IndexedMass, MassIndexingEngine, and mass indexing method in the engine

* Added GetXIC method in MassIndexingEngine and tests

* added some tests

* removed unnecessary override methods; Added optional charge in generic peak finding method so it can work with all kinds of indexed peaks

* ppmtol

* bug

* error

* converted all doubles into float in IndexedMass

* added comments and exception check

* .

* .

* comments

* moved exception handling to GetPeakFromBin

* .

* merge

* added override GetAllXics in massIndexingEngine

* added more tests

* .

* moved exception check

* changed RT to  zerobasedScanIndex in getallXIcs

* check if peak is indexedmass in GetAllXics

---------

Co-authored-by: trishorts <[email protected]>
Co-authored-by: Nic Bollis <[email protected]>
Co-authored-by: Sugar Ray <[email protected]>
Co-authored-by: Alexander-Sol <[email protected]>

Author: 5 people

mzLib/FlashLFQ/PeakIndexingEngine/IFlashLfqIndexingEngine.cs

@@ -21,7 +21,8 @@ public interface IFlashLfqIndexingEngine
         /// </summary>
         /// <param name="mz"> the m/z of the peak to be searched for </param>
         /// <param name="zeroBasedScanIndex"> the zero based index of the scan where the peak is to be found </param>
-        public IIndexedPeak GetIndexedPeak(double mz, int zeroBasedScanIndex, Tolerance tolerance);
+        /// <param name="charge"> an optional parameter used only for IIndexedMass and massIndexingEngine; must be null for mz peak indexing </param>
+        public IIndexedPeak GetIndexedPeak(double mz, int zeroBasedScanIndex, Tolerance tolerance, int? charge = null);
         /// <summary>
         /// A generic method of peak tracing across the retention time. Finds peaks with a given mz that occur on either side of a given
         /// retention time. Peak searching iterates backwards through the scans until the peak 
@@ -34,8 +35,9 @@ public interface IFlashLfqIndexingEngine
         /// <param name="missedScansAllowed"> the number of successive missed scans allowed before the xic is terminated </param>
         /// <param name="maxPeakHalfWidth"> the maximum distance from the apex RT of the XIC to both start RT and end RT </param>
         /// <param name="matchedPeaks"> the dictionary that stores all the peaks already matched to an xic </param>
+        /// <param name="charge"> an optional parameter used only for IIndexedMass and massIndexingEngine; must be null for mz peak indexing </param>
         /// <returns> A list of IIndexedPeak objects, ordered by retention time </returns>
-        public List<IIndexedPeak> GetXic(double mz, double retentionTime, Tolerance ppmTolerance, int missedScansAllowed, double maxPeakHalfWidth = int.MaxValue, Dictionary<IIndexedPeak, ExtractedIonChromatogram> matchedPeaks = null);
+        public List<IIndexedPeak> GetXic(double mz, double retentionTime, Tolerance ppmTolerance, int missedScansAllowed, double maxPeakHalfWidth = int.MaxValue, int? charge = null, Dictionary<IIndexedPeak, ExtractedIonChromatogram> matchedPeaks = null);
         /// <summary>
         /// Clear the indexed peaks and the jagged array of indexed peaks to free up memory
         /// </summary>

mzLib/MassSpectrometry/PeakIndexing/ExtractedIonChromatogram.cs

@@ -15,7 +15,7 @@ namespace MassSpectrometry
     /// </summary>
     public class ExtractedIonChromatogram
     {
-        public List<IIndexedPeak> Peaks { get; set; }
+        public virtual List<IIndexedPeak> Peaks { get; set; }
 
         public double ApexRT;
         public int ApexScanIndex;

mzLib/MassSpectrometry/PeakIndexing/IndexedMass.cs

@@ -0,0 +1,30 @@
+using System;
+using System.Collections.Generic;
+using System.Linq;
+using System.Text;
+using System.Threading.Tasks;
+
+namespace MassSpectrometry
+{
+    public class IndexedMass : IIndexedPeak
+    {
+        public float Intensity { get; set; }
+        public float RetentionTime { get; set; }
+        public int ZeroBasedScanIndex { get; set; }
+        public float M { get; set; }
+        public int Charge { get; set; } 
+        public int MsLevel { get; set; } 
+        public IsotopicEnvelope IsotopicEnvelope { get; set; }
+
+        public IndexedMass(IsotopicEnvelope envelope, double retentionTime, int zeroBasedScanIndex, int msLevel)
+        {
+            IsotopicEnvelope = envelope;
+            Intensity = envelope.Peaks.Max(p => (float)p.intensity);
+            RetentionTime = (float)retentionTime;
+            ZeroBasedScanIndex = zeroBasedScanIndex;
+            M = (float)envelope.MonoisotopicMass;
+            Charge = envelope.Charge;
+            MsLevel = msLevel;
+        }
+    }
+}

mzLib/MassSpectrometry/PeakIndexing/IndexingEngine.cs

@@ -1,5 +1,4 @@
-using Chemistry;
-using MzLibUtil;
+using MzLibUtil;
 using System;
 using System.Collections.Generic;
 using System.Linq;
@@ -19,7 +18,7 @@ public abstract class IndexingEngine<T> where T : IIndexedPeak
         /// </summary>
         protected List<T>[]? IndexedPeaks;
         protected virtual int BinsPerDalton => 100;
-        public ScanInfo[]? ScanInfoArray { get; private set; }
+        public ScanInfo[]? ScanInfoArray { get; protected set; }
 
         /// <summary>
         /// Read in all spectral peaks from scans, index the peaks and store them in a list ordered by m/z
@@ -56,26 +55,19 @@ public virtual bool IndexPeaks(MsDataScan[] scanArray)
                 return true;
         }
 
-        /// <summary>
-        /// A generic method for finding the closest peak with a specified mass and charge state and in a specified scan. Returns null if no peaks within tolerance are found.
-        /// </summary>
-        /// <param name="mz"> the m/z of the peak to be searched for </param>
-        /// <param name="zeroBasedScanIndex"> the zero based index of the scan where the peak is to be found </param>
-        public IIndexedPeak? GetIndexedPeak(double theorMass, int zeroBasedScanIndex, Tolerance ppmTolerance, int chargeState) =>
-            GetIndexedPeak(theorMass.ToMz(chargeState), zeroBasedScanIndex, ppmTolerance);
-
         /// <summary>
         /// A generic method for finding the closest peak with a specified m/z and in a specified scan. Returns null if no peaks within tolerance are found.
         /// </summary>
-        /// <param name="mz"> the m/z of the peak to be searched for </param>
+        /// <param name="m"> the m/z of the peak to be searched for </param>
         /// <param name="zeroBasedScanIndex"> the zero based index of the scan where the peak is to be found </param>
-        public IIndexedPeak? GetIndexedPeak(double mz, int zeroBasedScanIndex, Tolerance ppmTolerance)
+        /// <param name="charge"> an optional parameter used only for IIndexedMass and massIndexingEngine; must be null for mz peak indexing </param>
+        public IIndexedPeak? GetIndexedPeak(double m, int zeroBasedScanIndex, Tolerance ppmTolerance, int? charge = null)
         {
             if (IndexedPeaks == null) throw new MzLibException("Error: Attempt to retrieve indexed peak before peak indexing was performed");
-            var bins = GetBinsInRange(mz, ppmTolerance);
+            var bins = GetBinsInRange(m, ppmTolerance);
             if (bins.Count == 0) return default(T);
             List<int> peakIndicesInBins = bins.Select(b => BinarySearchForIndexedPeak(b, zeroBasedScanIndex)).ToList();
-            return GetBestPeakFromBins(bins, mz, zeroBasedScanIndex, peakIndicesInBins, ppmTolerance);
+            return GetBestPeakFromBins(bins, m, zeroBasedScanIndex, peakIndicesInBins, ppmTolerance, charge);
         }
 
         /// <summary>
@@ -85,17 +77,20 @@ public virtual bool IndexPeaks(MsDataScan[] scanArray)
         /// missedScansAllowed parameter. Missed scans don't have to be sequential. The same procedure
         /// is then repeated in the forward direction.
         /// </summary>
-        /// <param name="mz"> the m/z of the peak to be searched for </param>
+        /// <param name="m"> the m/z of the peak to be searched for </param>
         /// <param name="retentionTime"> the retention time where peak searching will begin </param>
         /// <param name="missedScansAllowed"> the number of successive missed scans allowed before the xic is terminated </param>
         /// <param name="maxPeakHalfWidth"> the maximum distance from the apex RT of the XIC to both start RT and end RT </param>
+        /// <param name="charge"> an optional parameter used only for IIndexedMass and massIndexingEngine; must be null for mz peak indexing </param>
         /// <param name="matchedPeaks"> the dictionary that stores all the peaks already matched to an xic </param>
         /// <returns> A list of IIndexedPeak objects, ordered by retention time </returns>
-        public List<IIndexedPeak> GetXic(double mz, double retentionTime, Tolerance ppmTolerance, int missedScansAllowed, double maxPeakHalfWidth = double.MaxValue, Dictionary<IIndexedPeak, ExtractedIonChromatogram> matchedPeaks = null)
+        public List<IIndexedPeak> GetXic(double m, double retentionTime, Tolerance ppmTolerance,
+            int missedScansAllowed, double maxPeakHalfWidth = double.MaxValue, int? charge = null, Dictionary<IIndexedPeak, ExtractedIonChromatogram> matchedPeaks = null)
         {
             // get precursor scan to start at
             int scanIndex = -1;
             if (ScanInfoArray == null) throw new MzLibException("Error: Attempt to retrieve XIC before peak indexing was performed");
+
             foreach (ScanInfo scan in ScanInfoArray)
             {
                 if (scan.RetentionTime < retentionTime)
@@ -108,7 +103,7 @@ public List<IIndexedPeak> GetXic(double mz, double retentionTime, Tolerance ppmT
                 }
             }
 
-            return GetXic(mz, scanIndex, ppmTolerance, missedScansAllowed, maxPeakHalfWidth, matchedPeaks);
+            return GetXic(m, scanIndex, ppmTolerance, missedScansAllowed, maxPeakHalfWidth, charge, matchedPeaks);
         }
 
         /// <summary>
@@ -118,23 +113,25 @@ public List<IIndexedPeak> GetXic(double mz, double retentionTime, Tolerance ppmT
         /// missedScansAllowed parameter. Missed scans don't have to be sequential. The same procedure
         /// is then repeated in the forward direction.
         /// </summary>
-        /// <param name="mz"> the m/z of the peak to be searched for </param>
+        /// <param name="m"> the m/z of the peak to be searched for </param>
         /// <param name="zeroBasedStartIndex"> the scan where peak searching begins </param>
         /// <param name="missedScansAllowed"> the number of successive missed scans allowed before the xic is terminated </param>
         /// <param name="maxPeakHalfWidth"> the maximum distance from the apex RT of the XIC to both start RT and end RT </param>
+        /// <param name="charge"> an optional parameter used only for IIndexedMass and massIndexingEngine; must be null for mz peak indexing </param>
         /// <param name="matchedPeaks"> the dictionary that stores all the peaks already matched to an xic </param>
         /// <returns> A list of IIndexedPeak objects, ordered by retention time </returns>
-        public List<IIndexedPeak> GetXic(double mz, int zeroBasedStartIndex, Tolerance ppmTolerance, int missedScansAllowed, double maxPeakHalfWidth = double.MaxValue, Dictionary<IIndexedPeak, ExtractedIonChromatogram> matchedPeaks = null)
+        public List<IIndexedPeak> GetXic(double m, int zeroBasedStartIndex, Tolerance ppmTolerance, int missedScansAllowed, double maxPeakHalfWidth = double.MaxValue, int? charge = null, Dictionary<IIndexedPeak, ExtractedIonChromatogram> matchedPeaks = null)
         {
             if (IndexedPeaks == null || ScanInfoArray == null) throw new MzLibException("Error: Attempt to retrieve XIC before peak indexing was performed");
+
             List<IIndexedPeak> xic = new List<IIndexedPeak>();
-            var allBins = GetBinsInRange(mz, ppmTolerance);
+            var allBins = GetBinsInRange(m, ppmTolerance);
             if (allBins.Count == 0)
                 return xic;
 
             // For each bin, find + store a pointer to the current index
             int[] peakPointerArray = allBins.Select(b => BinarySearchForIndexedPeak(b, zeroBasedStartIndex)).ToArray();
-            var initialPeak = GetBestPeakFromBins(allBins, mz, zeroBasedStartIndex, peakPointerArray, ppmTolerance);
+            var initialPeak = GetBestPeakFromBins(allBins, m, zeroBasedStartIndex, peakPointerArray, ppmTolerance, charge);
 
             if (initialPeak.IsNotDefaultOrNull())
                 xic.Add(initialPeak);
@@ -175,7 +172,7 @@ public List<IIndexedPeak> GetXic(double mz, int zeroBasedStartIndex, Tolerance p
                     }
 
                     // Search for the next peak
-                    var nextPeak = GetBestPeakFromBins(allBins, mz, currentZeroBasedScanIndex, pointerArrayCopy, ppmTolerance);
+                    var nextPeak = GetBestPeakFromBins(allBins, m, currentZeroBasedScanIndex, pointerArrayCopy, ppmTolerance, charge);
 
                     // Add the peak to the XIC or increment the missed peaks
                     if (nextPeak == null || (matchedPeaks != null && matchedPeaks.ContainsKey(nextPeak)))
@@ -197,7 +194,7 @@ public List<IIndexedPeak> GetXic(double mz, int zeroBasedStartIndex, Tolerance p
         /// <summary>
         /// A generic method performing peak tracing for all the peaks in an indexingEngine and trying to find all XICs.
         /// <returns> A list of ExtractedIonChromatogram objects representing all XICs that can be found in an indexingEngine </returns>
-        public List<ExtractedIonChromatogram> GetAllXics(PpmTolerance peakFindingTolerance, int maxMissedScanAllowed, double maxRTRange, int numPeakThreshold)
+        public virtual List<ExtractedIonChromatogram> GetAllXics(Tolerance peakFindingTolerance, int maxMissedScanAllowed, double maxRTRange, int numPeakThreshold)
         {
             var xics = new List<ExtractedIonChromatogram>();
             var matchedPeaks = new Dictionary<IIndexedPeak, ExtractedIonChromatogram>();
@@ -206,7 +203,8 @@ public List<ExtractedIonChromatogram> GetAllXics(PpmTolerance peakFindingToleran
             {
                 if (!matchedPeaks.ContainsKey(peak))
                 {
-                    var peakList = GetXic(peak.M, peak.RetentionTime, peakFindingTolerance, maxMissedScanAllowed, maxRTRange, matchedPeaks);
+                    int? charge = peak is IndexedMass indexedMass ? indexedMass.Charge : null;
+                    var peakList = GetXic(peak.M, peak.RetentionTime, peakFindingTolerance, maxMissedScanAllowed, maxRTRange, charge, matchedPeaks);
                     if (peakList.Count >= numPeakThreshold)
                     {
                         var newXIC = new ExtractedIonChromatogram(peakList);
@@ -240,12 +238,20 @@ internal List<List<T>> GetBinsInRange(double mz, Tolerance ppmTolerance)
             return allBins;
         }
 
-        internal static T? GetBestPeakFromBins(List<List<T>> allBins, double mz, int zeroBasedScanIndex, IList<int> peakIndicesInBins, Tolerance ppmTolerance)
+        /// <summary>
+        /// <param name="charge"> an optional parameter used only for IIndexedMass and massIndexingEngine; must be null for mz peak indexing </param>
+        /// Returns the peak that is closest to the target mz from all possible bins
+        /// </summary>
+        internal static T? GetBestPeakFromBins(List<List<T>> allBins, double mz, int zeroBasedScanIndex, IList<int> peakIndicesInBins, Tolerance ppmTolerance, int? charge = null)
         {
             T? bestPeak = default(T);
+            if (charge != null && typeof(T) != typeof(IndexedMass))
+            {
+                throw new MzLibException("Error: Attempted to access a peak using a charge parameter, but the peaks do not have charge information available.");
+            }
             for (int i = 0; i < allBins.Count; i++)
             {
-                var tempPeak = GetPeakFromBin(allBins[i], mz, zeroBasedScanIndex, peakIndicesInBins[i], ppmTolerance);
+                var tempPeak = GetPeakFromBin(allBins[i], mz, zeroBasedScanIndex, peakIndicesInBins[i], ppmTolerance, charge);
                 if (tempPeak.IsDefaultOrNull()) continue;
                 // Check if the peak is within the tolerance and if it is closer to the target M than the current peak
                 if (bestPeak.IsDefaultOrNull() || Math.Abs(tempPeak.M - mz) < Math.Abs(bestPeak.M - mz))
@@ -257,9 +263,10 @@ internal List<List<T>> GetBinsInRange(double mz, Tolerance ppmTolerance)
         }
 
         /// <summary>
-        /// Returns the peak that is closest to the target mz
+        /// <param name="charge"> an optional parameter used only for IIndexedMass and massIndexingEngine; must be null for mz peak indexing </param>
+        /// Returns the peak that is closest to the target mz from one bin
         /// </summary>
-        internal static T GetPeakFromBin(List<T> bin, double mz, int zeroBasedScanIndex, int peakIndexInBin, Tolerance ppmTolerance)
+        internal static T GetPeakFromBin(List<T> bin, double mz, int zeroBasedScanIndex, int peakIndexInBin, Tolerance ppmTolerance, int? charge = null)
         {
             T? bestPeak = default(T);
             if (peakIndexInBin < 0 || peakIndexInBin >= bin.Count) return bestPeak;
@@ -274,7 +281,8 @@ internal static T GetPeakFromBin(List<T> bin, double mz, int zeroBasedScanIndex,
 
                 if (ppmTolerance.Within(peak.M, mz)
                     && peak.ZeroBasedScanIndex == zeroBasedScanIndex
-                    && (bestPeak.IsDefaultOrNull() || Math.Abs(peak.M - mz) < Math.Abs(bestPeak.M - mz)))
+                    && (bestPeak.IsDefaultOrNull() || Math.Abs(peak.M - mz) < Math.Abs(bestPeak.M - mz))
+                    && (charge == null || (peak is IndexedMass mass && mass.Charge == charge)))
                 {
                     bestPeak = peak;
                 }

mzLib/MassSpectrometry/PeakIndexing/MassIndexingEngine.cs

@@ -0,0 +1,55 @@
+using MathNet.Numerics.RootFinding;
+using MzLibUtil;
+using System;
+using System.Collections.Generic;
+using System.Linq;
+using System.Runtime.CompilerServices;
+using System.Text;
+using System.Threading.Tasks;
+
+namespace MassSpectrometry
+{
+    public class MassIndexingEngine: IndexingEngine<IndexedMass>
+    {
+        protected override int BinsPerDalton => 1;
+        public int MaxMass { get; set; } = 30000;
+
+        public MassIndexingEngine()
+        {
+        }
+
+        public static MassIndexingEngine? InitializeMassIndexingEngine(MsDataScan[] scanArray, DeconvolutionParameters deconParameters)
+        {
+            MassIndexingEngine newEngine = new();
+            if (newEngine.IndexPeaks(scanArray, deconParameters))
+                return newEngine;
+            return null;
+        }
+
+        public bool IndexPeaks(MsDataScan[] scanArray, DeconvolutionParameters deconParameters, MzRange mzRange = null, double minMass = 0, int minCharge = 1)
+        {
+            if (scanArray.IsNullOrEmpty() || scanArray.All(p => p == null))
+                return false;
+
+            IndexedPeaks = new List<IndexedMass> [MaxMass];
+            ScanInfoArray = new ScanInfo[scanArray.Length];
+            for (int scanIndex = 0; scanIndex < scanArray.Length; scanIndex++)
+            {
+                ScanInfoArray[scanIndex] = new ScanInfo(scanArray[scanIndex].OneBasedScanNumber, scanIndex, scanArray[scanIndex].RetentionTime, scanArray[scanIndex].MsnOrder);
+                var envelopes = Deconvoluter.Deconvolute(scanArray[scanIndex].MassSpectrum, deconParameters, mzRange);
+                foreach (var envelope in envelopes)
+                {
+                    if (envelope.MonoisotopicMass < minMass || envelope.Charge < minCharge)
+                        continue;
+                    int roundedMass = (int)Math.Round(envelope.MonoisotopicMass * BinsPerDalton, 0);
+                    IndexedPeaks[roundedMass] ??= new List<IndexedMass>();
+                    IndexedPeaks[roundedMass].Add(new IndexedMass(envelope, scanArray[scanIndex].RetentionTime, scanIndex, scanArray[scanIndex].MsnOrder));
+                }
+            }
+            if (IndexedPeaks == null || IndexedPeaks.Length == 0)
+                return false;
+            else
+                return true;
+        }
+    }
+}

mzLib/MassSpectrometry/PeakIndexing/PeakSpline/XicSpline.cs

@@ -52,9 +52,9 @@ public void SetXicSplineXYData(ExtractedIonChromatogram xic, bool cycle = false,
             }
             else
             {
-                peakRts = xic.Peaks.Select(p => p.RetentionTime).ToArray();
+                peakRts = xic.Peaks.Select(p => (float)p.RetentionTime).ToArray();
             }
-            var peakIntensities = xic.Peaks.Select(p => p.Intensity).ToArray();
+            var peakIntensities = xic.Peaks.Select(p => (float)p.Intensity).ToArray();
             xic.XYData = GetXicSplineData(peakRts, peakIntensities, start, end);
         }