IsoDec Deconvolution Algorithm

mzLib/Development/Deconvolution/TestCases/SinglePeakDeconvolutionTestCase.cs → mzLib/Development/DeconvolutionDevelopment/TestCases/SinglePeakDeconvolutionTestCase.cs

@@ -23,7 +23,7 @@ public class SinglePeakDeconvolutionTestCase
         public SinglePeakDeconvolutionTestCase(DeconvolutionParameters deconParameters, string sampleInformation, string spectrumPath, int scanNumber,
             double expectedMostAbundantObservedIsotopicMass, int expectedIonChargeState, double selectedIonMz, double precursorPpmMassTolerance)
         {
-
+            DeconvolutionParameters = deconParameters;
             SampleInformation = sampleInformation;
             ExpectedMostAbundantObservedIsotopicMass = expectedMostAbundantObservedIsotopicMass;
             ExpectedIonChargeState = expectedIonChargeState;

mzLib/Development/Development.csproj

@@ -23,13 +23,13 @@
   </ItemGroup>
 
   <ItemGroup>
-    <None Update="Deconvolution\TestData\Averaged_221110_CytoOnly.mzML">
+    <None Update="DeconvolutionDevelopment\TestData\Averaged_221110_CytoOnly.mzML">
       <CopyToOutputDirectory>PreserveNewest</CopyToOutputDirectory>
     </None>
-    <None Update="Deconvolution\TestData\Averaged_221110_HGHOnly.mzML">
+    <None Update="DeconvolutionDevelopment\TestData\Averaged_221110_HGHOnly.mzML">
       <CopyToOutputDirectory>PreserveNewest</CopyToOutputDirectory>
     </None>
-    <None Update="Deconvolution\TestData\Averaged_221110_UbiqOnly.mzML">
+    <None Update="DeconvolutionDevelopment\TestData\Averaged_221110_UbiqOnly.mzML">
       <CopyToOutputDirectory>PreserveNewest</CopyToOutputDirectory>
     </None>
   </ItemGroup>

mzLib/MassSpectrometry/Deconvolution/Algorithms/DeconvolutionAlgorithm.cs

@@ -1,8 +1,6 @@
 using System;
 using System.Collections.Generic;
 using System.Linq;
-using System.Text;
-using System.Threading.Tasks;
 using Chemistry;
 using MzLibUtil;
 

mzLib/MassSpectrometry/Deconvolution/Algorithms/IsoDecAlgorithm.cs

@@ -0,0 +1,159 @@
+using System;
+using System.Collections.Generic;
+using System.Linq;
+using System.Runtime.InteropServices;
+using MzLibUtil;
+
+namespace MassSpectrometry
+{
+    /// <summary>
+    /// Performs deconvolution on a single spectrum or region of spectrum using the Isodec algorithm
+    /// <remarks>
+    /// Isodec only needs to region of interest and does not use surrounding charge states as references.
+    /// Isodec can report multiple monoisotopic masses for a single peak if enabled by ReportMultipleMonoisos parameter
+    ///     In this case, the resulting isotopic envelopes will have the same precursor ID.
+    /// </remarks>
+    /// </summary>
+    internal class IsoDecAlgorithm : DeconvolutionAlgorithm
+    {
+        internal IsoDecAlgorithm(DeconvolutionParameters deconParameters) : base(deconParameters)
+        {
+
+        }
+
+        /// <summary>
+        /// Struct passed by pointer in memory to the Isodec.dll
+        /// </summary>
+        [StructLayout(LayoutKind.Sequential, Pack =1)]
+        public struct MatchedPeak
+        {
+            public float mz;
+            public int z;
+            public float monoiso;
+            public float peakmass;
+            public float avgmass;
+            public float area;
+            public float peakint;
+            [MarshalAs(UnmanagedType.ByValArray, SizeConst = 64)]
+            public float[] matchedindsiso;
+            [MarshalAs(UnmanagedType.ByValArray, SizeConst = 64)]
+            public float[] matchedindsexp;
+            [MarshalAs(UnmanagedType.ByValArray, SizeConst = 64)]
+            public float[] isomz;
+            [MarshalAs(UnmanagedType.ByValArray, SizeConst = 64)]
+            public float[] isodist;
+            [MarshalAs(UnmanagedType.ByValArray, SizeConst = 64)]
+            public float[] isomass;
+            [MarshalAs(UnmanagedType.ByValArray, SizeConst = 16)]
+            public float[] monoisos;
+            int startindex;
+            int endindex;
+            public float score;
+            public int realisolength;
+        }
+
+        /// <summary>
+        /// Calls the Isodec.dll to perform deconvolution on the given spectrum
+        /// The Isodec.dll requires three other dll's as dependencies: isogenmass.dll, libmmd.dll, scml_dispmd.dll
+        /// </summary>
+        /// <param name="cmz"></param>
+        /// <param name="cintensity"></param>
+        /// <param name="c"></param>
+        /// <param name="fname"></param>
+        /// <param name="matchedpeaks"></param>
+        /// <param name="settings"></param>
+        /// <returns></returns>
+
+        [DllImport("isodeclib.dll", EntryPoint = "process_spectrum", CallingConvention = CallingConvention.Cdecl)]
+        protected static extern int process_spectrum(double[] cmz, float[] cintensity, int c, string fname, IntPtr matchedpeaks, IsoDecDeconvolutionParameters.IsoSettings settings);
+
+        internal override IEnumerable<IsotopicEnvelope> Deconvolute(MzSpectrum spectrum, MzRange range)
+        {
+            var deconParams = DeconvolutionParameters as IsoDecDeconvolutionParameters ?? throw new MzLibException("Deconvolution params and algorithm do not match");
+
+            var firstIndex = spectrum.GetClosestPeakIndex(range.Minimum);
+            var lastIndex = spectrum.GetClosestPeakIndex(range.Maximum);
+
+            var mzs = spectrum.XArray[firstIndex..lastIndex]
+                .Select(p => p)
+                .ToArray();
+            var intensities = spectrum.YArray[firstIndex..lastIndex]
+                .Select(p => (float)p)
+                .ToArray();
+
+            var mpArray = new byte[intensities.Length * Marshal.SizeOf(typeof(MatchedPeak))];
+            GCHandle handle = GCHandle.Alloc(mpArray, GCHandleType.Pinned);
+            try
+            {
+                IntPtr matchedPeaksPtr = (IntPtr)handle.AddrOfPinnedObject();
+                IsoDecDeconvolutionParameters.IsoSettings settings = deconParams.ToIsoSettings();
+                int result = process_spectrum(mzs, intensities, intensities.Length, null, matchedPeaksPtr, settings);
+                if (result <= 0)
+                    return Enumerable.Empty<IsotopicEnvelope>();
+
+                // Handle results
+                MatchedPeak[] matchedpeaks = new MatchedPeak[result];
+                for (int i = 0; i < result; i++)
+                {
+                    matchedpeaks[i] = Marshal.PtrToStructure<MatchedPeak>(matchedPeaksPtr + i * Marshal.SizeOf(typeof(MatchedPeak)));
+                }
+
+                return ConvertToIsotopicEnvelopes(deconParams, matchedpeaks, spectrum);
+            }
+            finally
+            {
+                handle.Free();
+            }
+        }
+
+        /// <summary>
+        /// Converts the isodec output (MatchedPeak) to IsotopicEnvelope for return
+        /// </summary>
+        /// <param name="parameters"></param>
+        /// <param name="matchedpeaks"></param>
+        /// <param name="spectrum"></param>
+        /// <returns></returns>
+        private List<IsotopicEnvelope> ConvertToIsotopicEnvelopes(IsoDecDeconvolutionParameters parameters, MatchedPeak[] matchedpeaks, MzSpectrum spectrum)
+        {
+            List<IsotopicEnvelope> result = new List<IsotopicEnvelope>();
+            int currentId = 0;
+            var tolerance = new PpmTolerance(5);
+            foreach(MatchedPeak peak in matchedpeaks)
+            {
+                List<(double,double)> peaks = new List<(double,double)> ();
+                for (int i = 0; i < peak.realisolength; i++)
+                {
+
+                    List<int> indicesWithinTolerance = spectrum.GetPeakIndicesWithinTolerance(peak.isomz[i], tolerance);
+                    double maxIntensity = 0;
+                    int maxIndex = -1;
+                    foreach (int index in indicesWithinTolerance)
+                    {
+                        if (spectrum.YArray[index] > maxIntensity) { maxIntensity = spectrum.YArray[index]; maxIndex = index; }
+                    }
+                    if (maxIndex >= 0)
+                    {
+                        peaks.Add((spectrum.XArray[maxIndex], spectrum.YArray[maxIndex]));
+                    }
+                    else
+                    {
+                        peaks.Add((peak.isomz[i], 0));
+                    }
+
+                }
+                int charge = peak.z;
+                if(parameters.Polarity == Polarity.Negative) { charge = -peak.z; }
+                if(parameters.ReportMulitpleMonoisos)
+                {
+                    foreach (float monoiso in peak.monoisos)
+                    {
+                        if (monoiso > 0) { result.Add(new IsotopicEnvelope(currentId, peaks, (double)monoiso, charge, peak.peakint, peak.score)); }
+                    }
+                }
+                else { result.Add(new IsotopicEnvelope(currentId, peaks, (double)peak.monoiso, charge, peak.peakint, peak.score)); }
+                currentId++;
+            }
+            return result;
+        }
+    }
+}

mzLib/MassSpectrometry/Deconvolution/Deconvoluter.cs

@@ -4,12 +4,6 @@
 
 namespace MassSpectrometry
 {
-    public enum DeconvolutionType
-    {
-        ClassicDeconvolution,
-        ExampleNewDeconvolutionTemplate,
-    }
-
     /// <summary>
     /// Context class for all deconvolution
     /// </summary>
@@ -26,9 +20,7 @@ public static IEnumerable<IsotopicEnvelope> Deconvolute(MsDataScan scan,
             DeconvolutionParameters deconvolutionParameters, MzRange rangeToGetPeaksFrom = null)
         {
             // set any specific deconvolution parameters found only in the MsDataScan
-
-            foreach (var isotopicEnvelope in Deconvolute(scan.MassSpectrum, deconvolutionParameters, rangeToGetPeaksFrom)) 
-                yield return isotopicEnvelope;
+            return Deconvolute(scan.MassSpectrum, deconvolutionParameters, rangeToGetPeaksFrom);
         }
 
         /// <summary>
@@ -43,37 +35,53 @@ public static IEnumerable<IsotopicEnvelope> Deconvolute(MzSpectrum spectrum,
         {
             rangeToGetPeaksFrom ??= spectrum.Range;
 
+            // Short circuit deconvolution if it is called on a neutral mass spectrum
+            if (spectrum is NeutralMassSpectrum newt)
+                return DeconvoluteNeutralMassSpectrum(newt, rangeToGetPeaksFrom);
+
             // set deconvolution algorithm 
-            DeconvolutionAlgorithm deconAlgorithm;
-            switch (deconvolutionParameters.DeconvolutionType)
-            {
-                case DeconvolutionType.ClassicDeconvolution:
-                    deconAlgorithm = new ClassicDeconvolutionAlgorithm(deconvolutionParameters);
-                    break;
+            DeconvolutionAlgorithm deconAlgorithm = CreateAlgorithm(deconvolutionParameters);
 
-                case DeconvolutionType.ExampleNewDeconvolutionTemplate:
-                    deconAlgorithm = new ExampleNewDeconvolutionAlgorithmTemplate(deconvolutionParameters);
-                    break;
+            // Delegate deconvolution to the algorithm
+            return deconAlgorithm.Deconvolute(spectrum, rangeToGetPeaksFrom);
+        }
 
-                default: throw new MzLibException("DeconvolutionType not yet supported");
-            }
+        /// <summary>
+        /// Factory method to create the correct deconvolution algorithm from the parameters
+        /// </summary>
+        /// <param name="parameters"></param>
+        /// <returns></returns>
+        /// <exception cref="MzLibException"></exception>
+        private static DeconvolutionAlgorithm CreateAlgorithm(DeconvolutionParameters parameters)
+        {
+            return parameters.DeconvolutionType switch
+            {
+                DeconvolutionType.ClassicDeconvolution => new ClassicDeconvolutionAlgorithm(parameters),
+                DeconvolutionType.ExampleNewDeconvolutionTemplate => new ExampleNewDeconvolutionAlgorithmTemplate(parameters),
+                DeconvolutionType.IsoDecDeconvolution => new IsoDecAlgorithm(parameters),
+                _ => throw new MzLibException("DeconvolutionType not yet supported")
+            };
+        }
 
-            // Short circuit deconvolution if it is called on a neutral mass spectrum
-            if (spectrum is NeutralMassSpectrum newt)
+        /// <summary>
+        /// Returns all peaks in the neutral mass spectrum as an isotopic envelope with a single peak
+        /// </summary>
+        /// <param name="neutralSpectrum"></param>
+        /// <param name="range"></param>
+        /// <returns></returns>
+        private static IEnumerable<IsotopicEnvelope> DeconvoluteNeutralMassSpectrum(NeutralMassSpectrum neutralSpectrum, MzRange range)
+        {
+            for (int i = 0; i < neutralSpectrum.XArray.Length; i++)
             {
-                for (int i = 0; i < newt.XArray.Length; i++)
+                double neutralMass = neutralSpectrum.XArray[i];
+                double intensity = neutralSpectrum.YArray[i];
+                int chargeState = neutralSpectrum.Charges[i];
+
+                if (range.Contains(neutralMass.ToMz(chargeState)))
                 {
-                    // skip this peak if it's outside the range of interest (e.g. if we're only interested in deconvoluting a small m/z range)
-                    if (!rangeToGetPeaksFrom.Contains(newt.XArray[i].ToMz(newt.Charges[i])))
-                        continue; 
-                    yield return new IsotopicEnvelope(newt.XArray[i], newt.YArray[i], newt.Charges[i]);
+                    yield return new IsotopicEnvelope(neutralMass, intensity, chargeState);
                 }
             }
-            else
-            {
-                foreach (var isotopicEnvelope in deconAlgorithm.Deconvolute(spectrum, rangeToGetPeaksFrom)) 
-                    yield return isotopicEnvelope;
-            }
         }
     }
 }

mzLib/MassSpectrometry/Deconvolution/DeconvolutionType.cs

@@ -0,0 +1,15 @@
+using System;
+using System.Collections.Generic;
+using System.Linq;
+using System.Text;
+using System.Threading.Tasks;
+
+namespace MassSpectrometry
+{
+    public enum DeconvolutionType
+    {
+        ClassicDeconvolution,
+        ExampleNewDeconvolutionTemplate,
+        IsoDecDeconvolution,
+    }
+}