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/ClassicDeconvolutionAlgorithm.cs

@@ -10,15 +10,11 @@
 
 namespace MassSpectrometry
 {
-    public class ClassicDeconvolutionAlgorithm : DeconvolutionAlgorithm
+    public class ClassicDeconvolutionAlgorithm(DeconvolutionParameters deconParameters)
+        : DeconvolutionAlgorithm(deconParameters)
     {
         private MzSpectrum spectrum;
 
-        public ClassicDeconvolutionAlgorithm(DeconvolutionParameters deconParameters) : base(deconParameters)
-        {
-
-        }
-
         /// <summary>
         /// Override to deconvolute the spectra using the Classic Deconvolution algorithm
         /// </summary>

mzLib/MassSpectrometry/Deconvolution/Algorithms/IsoDecAlgorithm.cs

@@ -0,0 +1,182 @@
+using System;
+using System.Collections.Generic;
+using System.IO;
+using System.Linq;
+using System.Numerics;
+using System.Runtime.InteropServices;
+using System.Text;
+using System.Threading.Tasks;
+using MassSpectrometry.MzSpectra;
+using MzLibUtil;
+using MathNet.Numerics.Statistics;
+using System.Reflection.Metadata.Ecma335;
+
+namespace MassSpectrometry
+{
+    public class IsoDecAlgorithm(DeconvolutionParameters deconParameters)
+        : DeconvolutionAlgorithm(deconParameters)
+    {
+        private static string _phaseModelPath;
+        static IsoDecAlgorithm()
+        {
+            _phaseModelPath = Path.Combine(AppDomain.CurrentDomain.BaseDirectory, "Deconvolution", "Algorithms", "IsoDecResources",  "phase_model.bin");
+        }
+
+        [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 startmz;
+            public float endmz;
+            public float score;
+            public int realisolength;
+        }
+
+        public struct IsoSettings
+        {
+            public int phaseres; // Precision of encoding matrix
+            public int verbose; // Verbose output
+            public int peakwindow; // Peak Detection Window
+            public float peakthresh; // Peak Detection Threshold
+            public int minpeaks; // Minimum Peaks for an allowed peak
+            public float css_thresh; // Minimum cosine similarity score for isotope distribution
+            public float matchtol; // Match Tolerance for peak detection in ppm
+            public int maxshift; // Maximum shift allowed for isotope distribution
+            [MarshalAs(UnmanagedType.ByValArray, SizeConst =2)]
+            public float[] mzwindow; // MZ Window for isotope distribution
+            [MarshalAs(UnmanagedType.ByValArray, SizeConst = 2)]
+            public float[] plusoneintwindow; // Plus One Intensity range. Will be used for charge state 1
+            public int knockdown_rounds; // Number of knockdown rounds
+            public float min_score_diff; // Minimum score difference for isotope distribution to allow missed monoisotopic peaks
+            public float minareacovered; // Minimum area covered by isotope distribution. Use in or with css_thresh
+            public int isolength; // Isotope Distribution Length
+            public double mass_diff_c; // Mass difference between isotopes
+            public float adductmass; // Adduct Mass
+            public int minusoneaszero; // Use set the -1 isotope as 0 to help force better alignments
+            public float isotopethreshold; // Threshold for isotope distribution. Will remove relative intensities below this.
+            public float datathreshold; // Threshold for data. Will remove relative intensities below this relative to max intensity in each cluster
+            public float zscore_threshold; //Ratio above which a secondary charge state prediction will be returned.
+        }
+
+
+
+        [DllImport("Deconvolution/Algorithms/IsoDecResources/isodeclib.dll", CallingConvention = CallingConvention.Cdecl)]
+        public static extern int process_spectrum(double[] mz, float[] intensity, int len, string modelpath, IntPtr matchedpeaks, IsoSettings settings);
+
+        public override IEnumerable<IsotopicEnvelope> Deconvolute(MzSpectrum spectrum, MzRange range)
+        {
+            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();
+
+            IntPtr matchedPeaksPtr = Marshal.AllocHGlobal(Marshal.SizeOf(typeof(MatchedPeak)) * intensities.Length);
+            IsoSettings settings = DeconParametersToIsoSettings(deconParameters);
+            int result = process_spectrum(mzs, intensities, intensities.Length, _phaseModelPath , matchedPeaksPtr, settings);
+            if(result > 0)
+            {
+                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(matchedpeaks, spectrum);
+            }
+
+            else return new List<IsotopicEnvelope>();
+        }
+
+        private List<IsotopicEnvelope> ConvertToIsotopicEnvelopes(MatchedPeak[] matchedpeaks, MzSpectrum spectrum)
+        {
+            List<IsotopicEnvelope> result = new List<IsotopicEnvelope>();
+            foreach(MatchedPeak peak in matchedpeaks)
+            {
+                List<(double,double)> peaks = new List<(double,double)> ();
+                List<double> listofratios = new List<double>();
+                for (int i = 0; i < peak.realisolength; i++)
+                {
+
+                    List<int> indicesWithinTolerance = spectrum.GetPeakIndicesWithinTolerance(peak.isomz[i], new PpmTolerance(5));
+                    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)
+                    {
+                        listofratios.Add(peak.isodist[i] / spectrum.YArray[maxIndex]);
+                        peaks.Add((spectrum.XArray[maxIndex], spectrum.YArray[maxIndex]));
+                    }
+                    else
+                    {
+                        listofratios.Add(0);
+                        peaks.Add((peak.isomz[i], 0));
+                    }
+
+                }
+                //One of the features of our algorithm is that it reports multiple possible monoisotopics if they score highly enough.
+                //This should do away with having to allow missed monoisotopics randomly.
+                foreach (float monoiso in peak.monoisos)
+                {
+                    if(monoiso > 0) { result.Add(new IsotopicEnvelope(peaks, (double)monoiso, peak.z, peak.peakint, Statistics.StandardDeviation(listofratios), 0)); }
+
+                }                
+            }
+            return result;
+        }
+
+        public IsoSettings DeconParametersToIsoSettings(DeconvolutionParameters parameters)
+        {
+            IsoSettings result = new IsoSettings();
+            result.phaseres = 8;
+            result.verbose = 1;
+            result.peakwindow = 80;
+            result.peakthresh = (float)0.0001;
+            result.minpeaks = 3;
+            result.css_thresh = (float)0.7;
+            result.matchtol = 5;
+            result.maxshift = 3;
+            result.mzwindow = [(float)-1.05, (float)2.05];
+            result.plusoneintwindow = [(float)0.1, (float)0.6];
+            result.knockdown_rounds = 5;
+            result.min_score_diff = (float)0.1;
+            result.minareacovered = (float)0.15;
+            result.isolength = 64;
+            result.mass_diff_c = 1.0033;
+            //If polarity is positive, adduct is a proton, if negative, it's the loss of a proton.
+            if(parameters.Polarity == Polarity.Positive) { result.adductmass = (float)1.007276467; }
+            else { result.adductmass = (float)-1.007276467; }
+            result.minusoneaszero = 1;
+            result.isotopethreshold = (float)0.01;
+            result.datathreshold = (float)0.05;
+            result.zscore_threshold = (float)0.95;
+            return result;
+        }
+    }
+}

mzLib/MassSpectrometry/Deconvolution/Deconvoluter.cs

@@ -9,11 +9,7 @@
 
 namespace MassSpectrometry
 {
-    public enum DeconvolutionType
-    {
-        ClassicDeconvolution,
-        ExampleNewDeconvolutionTemplate,
-    }
+
 
     /// <summary>
     /// Context class for all deconvolution
@@ -44,6 +40,10 @@ public static IEnumerable<IsotopicEnvelope> Deconvolute(MsDataScan scan,
                     deconAlgorithm = new ExampleNewDeconvolutionAlgorithmTemplate(deconvolutionParameters);
                     break;
 
+                case DeconvolutionType.IsoDecDeconvolution:
+                    deconAlgorithm = new IsoDecAlgorithm(deconvolutionParameters);
+                    break;
+
                 default: throw new MzLibException("DeconvolutionType not yet supported");
             }
 
@@ -65,19 +65,13 @@ public static IEnumerable<IsotopicEnvelope> Deconvolute(MzSpectrum spectrum,
             rangeToGetPeaksFrom ??= spectrum.Range;
 
             // set deconvolution algorithm 
-            DeconvolutionAlgorithm deconAlgorithm;
-            switch (deconvolutionParameters.DeconvolutionType)
+            DeconvolutionAlgorithm deconAlgorithm = deconvolutionParameters.DeconvolutionType switch
             {
-                case DeconvolutionType.ClassicDeconvolution:
-                    deconAlgorithm = new ClassicDeconvolutionAlgorithm(deconvolutionParameters);
-                    break;
-
-                case DeconvolutionType.ExampleNewDeconvolutionTemplate:
-                    deconAlgorithm = new ExampleNewDeconvolutionAlgorithmTemplate(deconvolutionParameters);
-                    break;
-
-                default: throw new MzLibException("DeconvolutionType not yet supported");
-            }
+                DeconvolutionType.ClassicDeconvolution => new ClassicDeconvolutionAlgorithm(deconvolutionParameters),
+                DeconvolutionType.ExampleNewDeconvolutionTemplate => new ExampleNewDeconvolutionAlgorithmTemplate(deconvolutionParameters),
+                DeconvolutionType.IsoDecDeconvolution => new IsoDecAlgorithm(deconvolutionParameters),
+                _ => throw new MzLibException("DeconvolutionType not yet supported")
+            };
 
             return deconAlgorithm.Deconvolute(spectrum, rangeToGetPeaksFrom);
         }

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,
+    }
+}

mzLib/MassSpectrometry/Deconvolution/Parameters/IsoDecDeconvolutionParameters.cs

@@ -0,0 +1,14 @@
+using System;
+using System.Collections.Generic;
+using System.Linq;
+using System.Text;
+using System.Threading.Tasks;
+
+namespace MassSpectrometry
+{
+    public class IsoDecDeconvolutionParameters(Polarity polarity = Polarity.Positive)
+        : DeconvolutionParameters(1, 100, polarity)
+    {
+        public override DeconvolutionType DeconvolutionType { get; protected set; } = DeconvolutionType.IsoDecDeconvolution;
+    }
+}

mzLib/MassSpectrometry/MassSpectrometry.csproj

@@ -9,6 +9,8 @@
     <DebugType>full</DebugType>
     <DebugSymbols>true</DebugSymbols>
   </PropertyGroup>
+
+
 
   <ItemGroup>
     <PackageReference Include="CsvHelper" Version="32.0.3" />
@@ -20,4 +22,26 @@
     <ProjectReference Include="..\MzLibUtil\MzLibUtil.csproj" />
   </ItemGroup>
 
+	<ItemGroup>
+		<None Remove="Deconvolution\Algorithms\isodeclib.dll" />
+		<None Remove="Deconvolution\Algorithms\libmmd.dll" />
+		<None Remove="Deconvolution\Algorithms\phase_model.bin" />
+		<None Remove="Deconvolution\Algorithms\svml_dispmd.dll" />
+	</ItemGroup>
+
+	<ItemGroup>
+		<Content Include="Deconvolution\Algorithms\IsoDecResources\isodeclib.dll">
+		  <CopyToOutputDirectory>Always</CopyToOutputDirectory>
+		</Content>
+		<Content Include="Deconvolution\Algorithms\IsoDecResources\libmmd.dll">
+		  <CopyToOutputDirectory>Always</CopyToOutputDirectory>
+		</Content>
+		<Content Include="Deconvolution\Algorithms\IsoDecResources\phase_model.bin">
+			<CopyToOutputDirectory>Always</CopyToOutputDirectory>
+		</Content>
+		<Content Include="Deconvolution\Algorithms\IsoDecResources\svml_dispmd.dll">
+		  <CopyToOutputDirectory>Always</CopyToOutputDirectory>
+		</Content>
+	</ItemGroup>
+
 </Project>

mzLib/MassSpectrometry/MzSpectra/MzSpectrum.cs

@@ -611,6 +611,36 @@
             return XArray.GetClosestIndex(x);
         }
 
+        public List<int> GetPeakIndicesWithinTolerance(double x, PpmTolerance tolerance)
+        {
+            List<int> indices = new List<int>();
+            int nearestIdx = XArray.GetClosestIndex(x);
+            if (tolerance.Within(XArray[nearestIdx], x)) { indices.Add(nearestIdx); }
+            int shift = 1;
+
+            bool findingUpper = true;
+            bool findingLower = true;
+
+            while(findingLower || findingLower)
+            {
+                int upperIdx = nearestIdx+ shift;
+                int lowerIdx = nearestIdx-shift;
+                if(findingUpper && upperIdx < XArray.Length)
+                {
+                    if(tolerance.Within(XArray[upperIdx], x)) {indices.Add(upperIdx); }
+                    else { findingUpper = false; }
+                }
+
+                if(findingLower && lowerIdx >= 0)
+                {
+                    if (tolerance.Within(XArray[lowerIdx], x)) { indices.Add(lowerIdx); }
+                    else { findingLower = false; }
+                }
+                shift++;
+            }
+            return indices;
+        }
+
         public void ReplaceXbyApplyingFunction(Func<MzPeak, double> convertor)
         {
             for (int i = 0; i < Size; i++)
@@ -764,14 +794,14 @@
             return numerator / denominator;
         }
 
         public bool Equals(MzSpectrum? other)
         {
             if (ReferenceEquals(null, other)) return false;
             if (ReferenceEquals(this, other)) return true;
             return (XArray.SequenceEqual(other.XArray) && YArray.SequenceEqual(other.YArray));
         }
 
         public override bool Equals(object? spectrumToCompare)
         {
             if (ReferenceEquals(null, spectrumToCompare)) return false;
             if (ReferenceEquals(this, spectrumToCompare)) return true;