Merge branch 'master' into SpectralDecon

mzLib/FlashLFQ/ChromatographicPeak.cs

@@ -1,4 +1,5 @@
 using Chemistry;
+using MathNet.Numerics.Statistics;
 using System;
 using System.Collections.Generic;
 using System.Linq;
@@ -14,7 +15,6 @@ public class ChromatographicPeak
         public double SplitRT;
         public readonly bool IsMbrPeak;
         public double MbrScore;
-        public double PosteriorErrorProbability { get { return NumIdentificationsByFullSeq > 1 ? 1 : Identifications.Min(p => p.PosteriorErrorProbability); } }
 
         public ChromatographicPeak(Identification id, bool isMbrPeak, SpectraFileInfo fileInfo)
         {
@@ -35,6 +35,18 @@ public ChromatographicPeak(Identification id, bool isMbrPeak, SpectraFileInfo fi
         public int NumIdentificationsByBaseSeq { get; private set; }
         public int NumIdentificationsByFullSeq { get; private set; }
         public double MassError { get; private set; }
+        /// <summary>
+        /// Expected retention time for MBR acceptor peaks (mean)
+        /// </summary>
+        public double? RtHypothesis { get; private set; }
+        /// <summary>
+        /// Std. Dev of retention time differences between MBR acceptor file and donor file, used if # calibration points < 6
+        /// </summary>
+        public double? RtStdDev { get; private set;  }
+        /// <summary>
+        /// Interquartile range of retention time differences between MBR acceptor file and donor file, used if # calibration points >= 6
+        /// </summary>
+        public double? RtInterquartileRange { get; private set; }
 
         public static string TabSeparatedHeader
         {
@@ -68,6 +80,19 @@ public static string TabSeparatedHeader
             }
         }
 
+        /// <summary>
+        /// Sets retention time information for a given peak. Used for MBR peaks
+        /// </summary>
+        /// <param name="rtHypothesis"> Expected retention time for peak, based on alignment between a donor and acceptor file </param>
+        /// <param name="rtStdDev"> Standard deviation in the retention time differences between aligned peaks </param>
+        /// <param name="rtInterquartileRange"> Interquartile range og the retention time differences between aligned peaks</param>
+        internal void SetRtWindow(double rtHypothesis, double? rtStdDev, double? rtInterquartileRange)
+        {
+            RtHypothesis = rtHypothesis;
+            RtStdDev = rtStdDev;
+            RtInterquartileRange = rtInterquartileRange;
+        }
+
         public void CalculateIntensityForThisFeature(bool integrate)
         {
             if (IsotopicEnvelopes.Any())

mzLib/FlashLFQ/FlashLfqEngine.cs

@@ -676,23 +676,27 @@ private void QuantifyMatchBetweenRunsPeaks(SpectraFileInfo idAcceptorFile)
                                 .Select(p => p.AcceptorFilePeak.Apex.IndexedPeak.RetentionTime - p.DonorFilePeak.Apex.IndexedPeak.RetentionTime)
                                 .ToList();
 
-                            double medianIntensityDiff = nearbyCalibrationPoints
-                                .Select(p => Math.Log(p.AcceptorFilePeak.Intensity, 2) - Math.Log(p.DonorFilePeak.Intensity, 2))
-                                .Median();
+                            //double medianIntensityDiff = nearbyCalibrationPoints
+                            //    .Select(p => Math.Log(p.AcceptorFilePeak.Intensity, 2) - Math.Log(p.DonorFilePeak.Intensity, 2))
+                            //    .Median();
 
                             // figure out the range of RT differences between the files that are "reasonable", centered around the median difference
                             double median = rtDiffs.Median();
 
                             // default range (if only 1 datapoint, or SD is 0, range is very high, etc)
                             double rtRange = MbrRtWindow;
+                            double? rtStdDev = null;
+                            double? rtInterquartileRange = null;
 
                             if (nearbyCalibrationPoints.Count < 6 && nearbyCalibrationPoints.Count > 1 && rtDiffs.StandardDeviation() > 0)
                             {
-                                rtRange = rtDiffs.StandardDeviation() * 6;
+                                rtStdDev = rtDiffs.StandardDeviation();
+                                rtRange = (double)rtStdDev * 6.0; // Multiplication inherited from legacy code, unsure of reason for 6
                             }
                             else if (nearbyCalibrationPoints.Count >= 6 && rtDiffs.InterquartileRange() > 0)
                             {
-                                rtRange = rtDiffs.InterquartileRange() * 4.5;
+                                rtInterquartileRange = rtDiffs.InterquartileRange();
+                                rtRange = (double)rtInterquartileRange * 4.5; // Multiplication inherited from legacy code, unsure of reason for 4.5
                             }
 
                             rtRange = Math.Min(rtRange, MbrRtWindow);
@@ -761,9 +765,10 @@ private void QuantifyMatchBetweenRunsPeaks(SpectraFileInfo idAcceptorFile)
 
                                     var xic = Peakfind(seedEnv.IndexedPeak.RetentionTime, donorIdentification.PeakfindingMass, z, idAcceptorFile, mbrTol);
                                     List<IsotopicEnvelope> bestChargeEnvelopes = GetIsotopicEnvelopes(xic, donorIdentification, z);
-
                                     acceptorPeak.IsotopicEnvelopes.AddRange(bestChargeEnvelopes);
                                     acceptorPeak.CalculateIntensityForThisFeature(Integrate);
+                                    acceptorPeak.SetRtWindow(acceptorFileRtHypothesis, rtStdDev, rtInterquartileRange);
+
                                     CutPeak(acceptorPeak, seedEnv.IndexedPeak.RetentionTime);
 
                                     var claimedPeaks = new HashSet<IndexedMassSpectralPeak>(acceptorPeak.IsotopicEnvelopes.Select(p => p.IndexedPeak));

mzLib/MassSpectrometry/Deconvolution/Algorithms/ClassicDeconvolutionAlgorithm.cs

@@ -23,8 +23,9 @@ public ClassicDeconvolutionAlgorithm(DeconvolutionParameters deconParameters) :
         /// Override to deconvolute the spectra using the Classic Deconvolution algorithm
         /// </summary>
         /// <param name="spectrumToDeconvolute">spectrum to deconvolute</param>
+        /// <param name="range">Range of peaks to deconvolute</param>
         /// <returns></returns>
-        public override IEnumerable<IsotopicEnvelope> Deconvolute(MzSpectrum spectrumToDeconvolute)
+        public override IEnumerable<IsotopicEnvelope> Deconvolute(MzSpectrum spectrumToDeconvolute, MzRange range)
         {
             var deconParams = DeconvolutionParameters as ClassicDeconvolutionParameters;
             if (deconParams == null)
@@ -41,7 +42,7 @@ public override IEnumerable<IsotopicEnvelope> Deconvolute(MzSpectrum spectrumToD
 
             var isolatedMassesAndCharges = new List<IsotopicEnvelope>();
 
-            (int start, int end) indexes = ExtractIndices(deconParams.Range.Minimum, deconParams.Range.Maximum);
+            (int start, int end) indexes = ExtractIndices(range.Minimum, range.Maximum);
 
             //find the most intense peak in the range
             double maxIntensity = 0;

mzLib/MassSpectrometry/Deconvolution/Algorithms/DeconvolutionAlgorithm.cs

@@ -77,8 +77,8 @@ protected DeconvolutionAlgorithm(DeconvolutionParameters deconParameters)
         /// Deconvolutes a mass spectrum
         /// </summary>
         /// <param name="spectrum">spectrum to be deconvoluted</param>
-        /// <param name="deconvolutionParams">parameters for the deconvolution</param>
+        /// <param name="range">Range of peaks to deconvolute</param>
         /// <returns></returns>
-        public abstract IEnumerable<IsotopicEnvelope> Deconvolute(MzSpectrum spectrum);
+        public abstract IEnumerable<IsotopicEnvelope> Deconvolute(MzSpectrum spectrum, MzRange range);
     }
 }

mzLib/MassSpectrometry/Deconvolution/Algorithms/ExampleNewDeconvolutionAlgorithm.cs

@@ -16,9 +16,11 @@ public ExampleNewDeconvolutionAlgorithm(DeconvolutionParameters deconParameters)
 
         }
 
-        public override IEnumerable<IsotopicEnvelope> Deconvolute(MzSpectrum spectrum)
+        public override IEnumerable<IsotopicEnvelope> Deconvolute(MzSpectrum spectrum, MzRange range = null)
         {
             var deconParams = DeconvolutionParameters as ExampleNewDeconvolutionParameters ?? throw new MzLibException("Deconvolution params and algorithm do not match");
+            range ??= spectrum.Range;
+
             throw new NotImplementedException();
         }
 

mzLib/MassSpectrometry/Deconvolution/Deconvoluter.cs

@@ -34,10 +34,13 @@ public Deconvoluter(DeconvolutionTypes deconType, DeconvolutionParameters deconP
         /// <summary>
         /// Deconvolute a MsDataScan
         /// </summary>
-        /// <param name="scan"></param>
+        /// <param name="scan">scan to deconvolute</param>
+        /// <param name="rangeToGetPeaksFrom">Range of peaks to deconvolute, if null, will deconvolute entire spectra</param>
         /// <returns></returns>
-        public IEnumerable<IsotopicEnvelope> Deconvolute(MsDataScan scan)
+        public IEnumerable<IsotopicEnvelope> Deconvolute(MsDataScan scan, MzRange rangeToGetPeaksFrom = null)
         {
+            rangeToGetPeaksFrom ??= scan.MassSpectrum.Range;
+
             // set deconvolution parameters that are only present in the MsDataScan
             switch (DeconvolutionType)
             {
@@ -52,7 +55,7 @@ public IEnumerable<IsotopicEnvelope> Deconvolute(MsDataScan scan)
                     break;
             }
 
-            return DeconvolutionAlgorithm.Deconvolute(scan.MassSpectrum);
+            return DeconvolutionAlgorithm.Deconvolute(scan.MassSpectrum, rangeToGetPeaksFrom);
         }
 
 

mzLib/MassSpectrometry/Deconvolution/DeconvoluterExtensions.cs

@@ -13,17 +13,25 @@ namespace MassSpectrometry
     /// </summary>
     public static class DeconvoluterExtensions
     {
+        /// <summary>
+        /// Deconvolutes a MzSpectrum object
+        /// </summary>
+        /// <param name="deconvoluter">performs deconvolution</param>
+        /// <param name="spectrum">spectrum to deconvolute</param>
+        /// <param name="range">mz range of returned peaks</param>
+        /// <returns></returns>
+        /// <exception cref="MzLibException"></exception>
         public static IEnumerable<IsotopicEnvelope> ClassicDeconvoluteMzSpectra(this Deconvoluter deconvoluter,
-            MzSpectrum spectrum, MzRange range)
+            MzSpectrum spectrum, MzRange range = null)
         {
             if (deconvoluter.DeconvolutionType != DeconvolutionTypes.ClassicDeconvolution)
             {
                 throw new MzLibException("Deconvoluter is not of correct type for this extension method");
             }
             else
             {
-                ((ClassicDeconvolutionParameters)deconvoluter.DeconvolutionParameters).Range = range;
-                return deconvoluter.DeconvolutionAlgorithm.Deconvolute(spectrum);
+                range ??= spectrum.Range;
+                return deconvoluter.DeconvolutionAlgorithm.Deconvolute(spectrum, range);
             }
         }
 

mzLib/MassSpectrometry/Deconvolution/Parameters/ClassicDeconvolutionParameters.cs

@@ -13,6 +13,9 @@ namespace MassSpectrometry
     public class ClassicDeconvolutionParameters : DeconvolutionParameters
     {
         public MzRange Range { get; set; }
+        public int MinAssumedChargeState { get; set; }
+        public int MaxAssumedChargeState { get; set; }
+        public double DeconvolutionTolerancePpm { get; set; }
         public double IntensityRatioLimit { get; set; }
 
         /// <summary>
@@ -26,7 +29,6 @@ public class ClassicDeconvolutionParameters : DeconvolutionParameters
         public ClassicDeconvolutionParameters(int minCharge, int maxCharge, double deconPpm, double intensityRatio, MzRange range = null) : 
             base (minCharge, maxCharge, deconPpm)
         {
-            Range = range;
             IntensityRatioLimit = intensityRatio;
         }
     }

mzLib/MassSpectrometry/MsDataScan.cs

@@ -129,8 +129,10 @@ public byte[] Get64BitNoiseDataBaseline()
         /// <returns></returns>
         public IEnumerable<IsotopicEnvelope> GetIsolatedMassesAndCharges(Deconvoluter deconvoluter)
         {
-            var allDeconvolutedEnvelopes = deconvoluter.Deconvolute(this);
-            return allDeconvolutedEnvelopes.Where(b => b.Peaks.Any(cc => isolationRange.Contains(cc.mz)));
+            // +- 8.5 are magic numbers from the original implementation of Classic Deconvolution
+            // it is believe that they represent the maximum mz space an isotopic envelopes can exist within
+            // This ensure that a peak is not cut in half by the mz isolation range
+            return deconvoluter.Deconvolute(this, new MzRange(IsolationRange.Minimum - 8.5, IsolationRange.Maximum + 8.5));
         }
 
         /// <summary>