Merge branch 'master' into StefanBranch

Author: stefanks

EngineLayer/Analysis/AnalysisEngine.cs

@@ -33,6 +33,9 @@ public class AnalysisEngine : MyEngine
         private readonly bool doParsimony;
         private readonly bool noOneHitWonders;
         private readonly bool doHistogramAnalysis;
+        private readonly bool quantify;
+        private readonly double quantifyRtTol;
+        private readonly double quantifyPpmTol;
         private readonly List<ProductType> lp;
         private readonly InitiatorMethionineBehavior initiatorMethionineBehavior;
         private readonly List<string> nestedIds;
@@ -42,7 +45,7 @@ public class AnalysisEngine : MyEngine
 
         #region Public Constructors
 
-        public AnalysisEngine(PsmParent[][] newPsms, Dictionary<CompactPeptide, HashSet<PeptideWithSetModifications>> compactPeptideToProteinPeptideMatching, List<Protein> proteinList, List<ModificationWithMass> variableModifications, List<ModificationWithMass> fixedModifications, List<ModificationWithMass> localizeableModifications, Protease protease, List<SearchMode> searchModes, IMsDataFile<IMsDataScan<IMzSpectrum<IMzPeak>>> myMSDataFile, Tolerance fragmentTolerance, Action<BinTreeStructure, string> action1, Action<List<NewPsmWithFdr>, string> action2, Action<List<ProteinGroup>, string> action3, bool doParsimony, bool noOneHitWonders, int maximumMissedCleavages, int maxModIsoforms, bool doHistogramAnalysis, List<ProductType> lp, double binTol, InitiatorMethionineBehavior initiatorMethionineBehavior, List<string> nestedIds)
+        public AnalysisEngine(PsmParent[][] newPsms, Dictionary<CompactPeptide, HashSet<PeptideWithSetModifications>> compactPeptideToProteinPeptideMatching, List<Protein> proteinList, List<ModificationWithMass> variableModifications, List<ModificationWithMass> fixedModifications, List<ModificationWithMass> localizeableModifications, Protease protease, List<SearchMode> searchModes, IMsDataFile<IMsDataScan<IMzSpectrum<IMzPeak>>> myMSDataFile, Tolerance fragmentTolerance, Action<BinTreeStructure, string> action1, Action<List<NewPsmWithFdr>, string> action2, Action<List<ProteinGroup>, string> action3, bool doParsimony, bool noOneHitWonders, int maximumMissedCleavages, int maxModIsoforms, bool doHistogramAnalysis, List<ProductType> lp, double binTol, InitiatorMethionineBehavior initiatorMethionineBehavior, List<string> nestedIds, bool Quantify, double QuantifyRtTol, double QuantifyPpmTol)
         {
             this.doParsimony = doParsimony;
             this.noOneHitWonders = noOneHitWonders;
@@ -65,6 +68,9 @@ public AnalysisEngine(PsmParent[][] newPsms, Dictionary<CompactPeptide, HashSet<
             this.lp = lp;
             this.binTol = binTol;
             this.initiatorMethionineBehavior = initiatorMethionineBehavior;
+            this.quantify = Quantify;
+            this.quantifyRtTol = QuantifyRtTol;
+            this.quantifyPpmTol = QuantifyPpmTol;
             this.nestedIds = nestedIds;
         }
 
@@ -74,8 +80,6 @@ public AnalysisEngine(PsmParent[][] newPsms, Dictionary<CompactPeptide, HashSet<
 
         public Dictionary<CompactPeptide, HashSet<PeptideWithSetModifications>> ApplyProteinParsimony(out List<ProteinGroup> proteinGroups)
         {
-            Status("Applying protein parsimony...", nestedIds);
-
             var proteinToPeptidesMatching = new Dictionary<Protein, HashSet<CompactPeptide>>();
             var parsimonyDict = new Dictionary<Protein, HashSet<CompactPeptide>>();
             var proteinsWithUniquePeptides = new Dictionary<Protein, HashSet<CompactPeptide>>();
@@ -483,6 +487,101 @@ public List<ProteinGroup> DoProteinFdr(List<ProteinGroup> proteinGroups)
             return sortedProteinGroups;
         }
 
+        public void RunQuantification(List<NewPsmWithFdr> psms, double rtTolerance, double ppmTolerance)
+        {
+            var baseSeqToPsmMatching = new Dictionary<string, HashSet<NewPsmWithFdr>>();
+
+            foreach (var psm in psms)
+            {
+                if (!baseSeqToPsmMatching.ContainsKey(psm.thisPSM.FullSequence))
+                {
+                    var psmList = new HashSet<NewPsmWithFdr>();
+                    psmList.Add(psm);
+                    baseSeqToPsmMatching.Add(psm.thisPSM.FullSequence, psmList);
+                }
+                else
+                {
+                    HashSet<NewPsmWithFdr> psmList;
+                    baseSeqToPsmMatching.TryGetValue(psm.thisPSM.FullSequence, out psmList);
+                    psmList.Add(psm);
+                }
+            }
+
+            foreach (var kvp in baseSeqToPsmMatching)
+            {
+                // calculate apex intensity
+                var rt1 = kvp.Value.Select(r => r.thisPSM.newPsm.scanRetentionTime).Min();
+                var rt2 = kvp.Value.Select(r => r.thisPSM.newPsm.scanRetentionTime).Max();
+                
+                double theoreticalMz = Chemistry.ClassExtensions.ToMz(kvp.Value.First().thisPSM.PeptideMonoisotopicMass, kvp.Value.First().thisPSM.newPsm.scanPrecursorCharge);
+
+                double mzTol = ((ppmTolerance / 1e6) * kvp.Value.First().thisPSM.PeptideMonoisotopicMass) / kvp.Value.First().thisPSM.newPsm.scanPrecursorCharge;
+
+                var spectraInThisWindow = myMsDataFile.GetMsScansInTimeRange(rt1 - rtTolerance, rt2 + rtTolerance).ToList();
+                var ms1SpectraInThisWindow = spectraInThisWindow.Where(s => s.MsnOrder == 1).ToList();
+
+                var intensities = new List<double>();
+                var retentionTimes = new List<double>();
+
+                foreach (var spectrum in ms1SpectraInThisWindow)
+                {
+                    var i = spectrum.MassSpectrum.Where(s => (Math.Abs(s.Mz - theoreticalMz) < mzTol)).ToList();
+                    foreach (var p in i)
+                    {
+                        intensities.Add(p.Intensity);
+                        retentionTimes.Add(spectrum.RetentionTime);
+                    }
+                }
+
+                double apexIntensity = 0;
+                if(intensities.Any())
+                    apexIntensity = intensities.Max();
+
+                foreach(var p in kvp.Value)
+                    p.thisPSM.newPsm.apexIntensity = apexIntensity;
+
+                // calculate full width half max (peak quality)
+                var apexIntensityIndex = intensities.IndexOf(apexIntensity);
+
+                double leftHalfMax = double.NaN;
+                double rightHalfMax = double.NaN;
+                double fullWidthHalfMax = double.NaN;
+
+                // left width half max
+                for (int i = apexIntensityIndex; i >= 0; i--)
+                {
+                    if (intensities[i] < (apexIntensity / 2))
+                    {
+                        leftHalfMax = retentionTimes[i];
+                        break;
+                    }
+                }
+
+                // right width half max
+                if (apexIntensityIndex >= 0)
+                {
+                    for (int i = apexIntensityIndex; i < intensities.Count; i++)
+                    {
+                        if (intensities[i] < (apexIntensity / 2))
+                        {
+                            rightHalfMax = retentionTimes[i];
+                            break;
+                        }
+                    }
+                }
+
+                if (!double.IsNaN(leftHalfMax) && !double.IsNaN(rightHalfMax))
+                {
+                    fullWidthHalfMax = rightHalfMax - leftHalfMax;
+                }
+
+                foreach (var p in kvp.Value)
+                    p.thisPSM.newPsm.fullWidthHalfMax = fullWidthHalfMax;
+
+                // calculate SNR (TODO**)
+            }
+        }
+
         #endregion Public Methods
 
         #region Protected Methods
@@ -551,6 +650,7 @@ protected override MyResults RunSpecific()
             List<ProteinGroup>[] proteinGroups = null;
             if (doParsimony)
             {
+                Status("Applying protein parsimony...", nestedIds);
                 proteinGroups = new List<ProteinGroup>[searchModes.Count];
                 // TODO**: make this faster (only apply parsimony once but make multiple instances of the same ProteinGroups
                 for (int i = 0; i < searchModes.Count; i++)
@@ -575,6 +675,13 @@ protected override MyResults RunSpecific()
 
                     Status("Running FDR analysis...", nestedIds);
                     var orderedPsmsWithFDR = DoFalseDiscoveryRateAnalysis(orderedPsmsWithPeptides, searchModes[j]);
+
+                    if (quantify)
+                    {
+                        Status("Quantifying peptides...", nestedIds);
+                        RunQuantification(orderedPsmsWithFDR, quantifyRtTol, quantifyPpmTol);
+                    }
+
                     writePsmsAction(orderedPsmsWithFDR, searchModes[j].FileNameAddition);
 
                     if (doHistogramAnalysis)
@@ -588,6 +695,7 @@ protected override MyResults RunSpecific()
                             writeHistogramPeaksAction(myTreeStructure, searchModes[j].FileNameAddition);
                         }
                     }
+
                     else
                     {
                         Status("Running FDR analysis on unique peptides...", nestedIds);

EngineLayer/ClassicSearch/ClassicSearchEngine.cs

@@ -150,7 +150,7 @@ protected override MyResults RunSpecific()
                                     //}
 
                                     var score = PsmWithMultiplePossiblePeptides.MatchIons(scan.TheScan, productMassTolerance, sortedProductMasses, matchedIonsArray);
-                                    var psm = new PsmClassic(yyy, fileName, scan.RetentionTime, scan.MonoisotopicPrecursorIntensity, scan.MonoisotopicPrecursorMass, scan.OneBasedScanNumber, scan.PrecursorCharge, scan.NumPeaks, scan.TotalIonCurrent, scan.MonoisotopicPrecursorMZ, score, theTuple.Item2);
+                                    var psm = new PsmClassic(yyy, fileName, scan.RetentionTime, scan.MonoisotopicPrecursorIntensity, scan.MonoisotopicPrecursorMass, scan.OneBasedScanNumber, scan.OneBasedPrecursorScanNumber, scan.PrecursorCharge, scan.NumPeaks, scan.TotalIonCurrent, scan.MonoisotopicPrecursorMZ, score, theTuple.Item2);
                                     if (psm.score > 1)
                                     {
                                         PsmClassic current_best_psm = psms[aede][scan.OneBasedScanNumber - 1];

EngineLayer/ClassicSearch/PsmClassic.cs

@@ -22,7 +22,7 @@ public class PsmClassic : PsmParent
 
         #region Public Constructors
 
-        public PsmClassic(PeptideWithSetModifications ps, string fileName, double scanRetentionTime, double scanPrecursorIntensity, double scanPrecursorMass, int scanNumber, int scanPrecursorCharge, int scanExperimentalPeaks, double totalIonCurrent, double scanPrecursorMZ, double score, int notch) : base(fileName, scanRetentionTime, scanPrecursorIntensity, scanPrecursorMass, scanNumber, scanPrecursorCharge, scanExperimentalPeaks, totalIonCurrent, scanPrecursorMZ, score, notch)
+        public PsmClassic(PeptideWithSetModifications ps, string fileName, double scanRetentionTime, double scanPrecursorIntensity, double scanPrecursorMass, int scanNumber, int scanPrecursorNumber, int scanPrecursorCharge, int scanExperimentalPeaks, double totalIonCurrent, double scanPrecursorMZ, double score, int notch) : base(fileName, scanRetentionTime, scanPrecursorIntensity, scanPrecursorMass, scanNumber, scanPrecursorNumber, scanPrecursorCharge, scanExperimentalPeaks, totalIonCurrent, scanPrecursorMZ, score, notch)
         {
             this.ps = ps;
         }

EngineLayer/LocalMs2Scan.cs

@@ -22,6 +22,7 @@ public LocalMS2Scan(IMsDataScanWithPrecursor<IMzSpectrum<IMzPeak>> b)
             MonoisotopicPrecursorMass = MonoisotopicPrecursorMZ.ToMass(PrecursorCharge);
 
             OneBasedScanNumber = b.OneBasedScanNumber;
+            OneBasedPrecursorScanNumber = b.OneBasedPrecursorScanNumber;
 
             RetentionTime = b.RetentionTime;
 
@@ -49,6 +50,7 @@ public LocalMS2Scan(double precursorMass)
 
         internal IMsDataScanWithPrecursor<IMzSpectrum<IMzPeak>> TheScan { get; private set; }
         internal int OneBasedScanNumber { get; private set; }
+        internal int OneBasedPrecursorScanNumber { get; private set; }
         internal double RetentionTime { get; private set; }
         internal int PrecursorCharge { get; private set; }
         internal int NumPeaks { get; private set; }

EngineLayer/ModernSearch/ModernSearchEngine.cs

@@ -137,7 +137,7 @@ protected override MyResults RunSpecific()
                         CompactPeptide theBestPeptide = bestPeptides[j];
                         if (theBestPeptide != null)
                         {
-                            newPsms[j][thisScan.OneBasedScanNumber - 1] = new PsmModern(theBestPeptide, Path.GetFileNameWithoutExtension(fileToSearch), thisScan.RetentionTime, thisScan.MonoisotopicPrecursorIntensity, thisScanprecursorMass, thisScan.OneBasedScanNumber, thisScan.PrecursorCharge, thisScan.NumPeaks, thisScan.TotalIonCurrent, thisScan.MonoisotopicPrecursorMZ, bestScores[j], bestNotches[j]);
+                            newPsms[j][thisScan.OneBasedScanNumber - 1] = new PsmModern(theBestPeptide, Path.GetFileNameWithoutExtension(fileToSearch), thisScan.RetentionTime, thisScan.MonoisotopicPrecursorIntensity, thisScanprecursorMass, thisScan.OneBasedScanNumber, thisScan.OneBasedPrecursorScanNumber, thisScan.PrecursorCharge, thisScan.NumPeaks, thisScan.TotalIonCurrent, thisScan.MonoisotopicPrecursorMZ, bestScores[j], bestNotches[j]);
                         }
                     }
                 }

EngineLayer/ModernSearch/PsmModern.cs

@@ -14,7 +14,7 @@ public class PsmModern : PsmParent
 
         #region Public Constructors
 
-        public PsmModern(CompactPeptide theBestPeptide, string fileName, double scanRetentionTime, double scanPrecursorIntensity, double scanPrecursorMass, int scanNumber, int scanPrecursorCharge, int scanExperimentalPeaks, double totalIonCurrent, double scanPrecursorMZ, double score, int notch) : base(fileName, scanRetentionTime, scanPrecursorIntensity, scanPrecursorMass, scanNumber, scanPrecursorCharge, scanExperimentalPeaks, totalIonCurrent, scanPrecursorMZ, score, notch)
+        public PsmModern(CompactPeptide theBestPeptide, string fileName, double scanRetentionTime, double scanPrecursorIntensity, double scanPrecursorMass, int scanNumber, int scanPrecursorNumber, int scanPrecursorCharge, int scanExperimentalPeaks, double totalIonCurrent, double scanPrecursorMZ, double score, int notch) : base(fileName, scanRetentionTime, scanPrecursorIntensity, scanPrecursorMass, scanNumber, scanPrecursorNumber, scanPrecursorCharge, scanExperimentalPeaks, totalIonCurrent, scanPrecursorMZ, score, notch)
         {
             compactPeptide = theBestPeptide;
         }

EngineLayer/ProteinGroup.cs

@@ -277,10 +277,6 @@ public void ScoreThisProteinGroup(List<ModificationWithMass> variableModificatio
             {
                 proteinGroupScore += psm.thisPSM.Score;
             }
-
-            //TODO**: give option to ignore protein groups with one peptide
-            //if (StrictPeptideList.Count == 1)
-            //    proteinGroupScore = 0;
         }
 
         public void CalculateSequenceCoverage()

EngineLayer/PsmParent.cs

@@ -14,6 +14,7 @@ public abstract class PsmParent
 
         public readonly string fileName;
         public readonly int scanNumber;
+        public readonly int precursorScanNumber;
         public readonly double scanRetentionTime;
         public readonly int scanExperimentalPeaks;
         public readonly double totalIonCurrent;
@@ -22,6 +23,8 @@ public abstract class PsmParent
         public readonly double scanPrecursorMZ;
         public readonly double scanPrecursorMass;
         public readonly double score;
+        public double apexIntensity;
+        public double fullWidthHalfMax;
 
         public Dictionary<ProductType, double[]> matchedIonsList;
         public List<double> LocalizedScores;
@@ -36,10 +39,11 @@ public abstract class PsmParent
 
         #region Protected Constructors
 
-        protected PsmParent(string fileName, double scanRetentionTime, double scanPrecursorIntensity, double scanPrecursorMass, int scanNumber, int scanPrecursorCharge, int scanExperimentalPeaks, double totalIonCurrent, double scanPrecursorMZ, double score, int notch)
+        protected PsmParent(string fileName, double scanRetentionTime, double scanPrecursorIntensity, double scanPrecursorMass, int scanNumber, int precursorScanNumber, int scanPrecursorCharge, int scanExperimentalPeaks, double totalIonCurrent, double scanPrecursorMZ, double score, int notch)
         {
             this.fileName = fileName;
             this.scanNumber = scanNumber;
+            this.precursorScanNumber = precursorScanNumber;
             this.scanRetentionTime = scanRetentionTime;
             this.scanExperimentalPeaks = scanExperimentalPeaks;
             this.totalIonCurrent = totalIonCurrent;
@@ -72,6 +76,8 @@ public override string ToString()
             sb.Append(scanPrecursorMass.ToString("F5", CultureInfo.InvariantCulture) + '\t');
             sb.Append(score.ToString("F3", CultureInfo.InvariantCulture) + '\t');
             sb.Append(notch.ToString("F3", CultureInfo.InvariantCulture) + '\t');
+            sb.Append(apexIntensity.ToString("F3", CultureInfo.InvariantCulture) + '\t');
+            sb.Append(fullWidthHalfMax.ToString("F3", CultureInfo.InvariantCulture) + '\t');
 
             sb.Append("[");
             foreach (var kvp in matchedIonsList)
@@ -112,6 +118,8 @@ internal static string GetTabSeparatedHeader()
             sb.Append("scanPrecursorMass" + '\t');
             sb.Append("score" + '\t');
             sb.Append("notch" + '\t');
+            sb.Append("apexIntensity" + '\t');
+            sb.Append("fullWidthHalfMax" + '\t');
 
             sb.Append("matched ions" + '\t');
             sb.Append("matched ion counts" + '\t');