Skip to content

IsoTracker performance improvement #942

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 9 commits into
base: master
Choose a base branch
from
Open

IsoTracker performance improvement #942

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
9 commits
Select commit Hold shift + click to select a range
4fcf16d
Changed parallelization in IsoTracker
Alexander-Sol Aug 21, 2025
cc6fa44
implemented binary search in BuildXIC
Alexander-Sol Aug 22, 2025
766d6ff
Decreased resolution in XicGroups
Alexander-Sol Aug 22, 2025
15da18f
reverted change to nuspec
Alexander-Sol Aug 22, 2025
233a000
Added new method for efficient XIC retrieval, test
Alexander-Sol Aug 22, 2025
248e75c
Merge branch 'master' into IsoTrackerSpeedup
trishorts Aug 22, 2025
e518ea8
Added additional logging
Alexander-Sol Aug 22, 2025
362570c
Merge branch 'IsoTrackerSpeedup' of https://github.com/Alexander-Sol/…
Alexander-Sol Aug 22, 2025
b62c1af
Merge branch 'master' into IsoTrackerSpeedup
trishorts Aug 26, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
236 changes: 138 additions & 98 deletions mzLib/FlashLFQ/FlashLfqEngine.cs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1821,97 +1821,131 @@ private void CutPeak(ChromatographicPeak peak, double identificationTime)
private void QuantifyIsobaricPeaks()
{
if(!FlashParams.Silent)
Console.WriteLine("Quantifying isobaric species...");
Console.WriteLine("Quantifying isobaric species for {0} isobar groups...", PeptideGroupsForIsoTracker.Count);
int isoGroupsSearched = 0;
double lastReportedProgress = 0;
double currentProgress = 0;


Parallel.ForEach(Partitioner.Create(0, PeptideGroupsForIsoTracker.Count),
new ParallelOptions { MaxDegreeOfParallelism = FlashParams.MaxThreads },
(range, loopState) =>
foreach(var idGroup in PeptideGroupsForIsoTracker)
{
Stopwatch stopwatch = null;
if(currentProgress == lastReportedProgress)
{
for (int i = range.Item1; i < range.Item2; i++)
{
var idGroup = PeptideGroupsForIsoTracker[i];
List<XIC> xicGroup = new List<XIC>();
var mostCommonChargeIdGroup = idGroup.Identifications.GroupBy(p => p.PrecursorChargeState).OrderBy(p => p.Count()).Last();
var id = mostCommonChargeIdGroup.First();
stopwatch = Stopwatch.StartNew();
}

// Try to get the primitive window for the XIC , from the (firstOne - 2) -> (lastOne + 2)
double rightWindow = mostCommonChargeIdGroup.Select(p => p.Ms2RetentionTimeInMinutes).Max() + 2;
double leftWindow = mostCommonChargeIdGroup.Select(p => p.Ms2RetentionTimeInMinutes).Min() - 2;
List<XIC> xicGroup = new List<XIC>();
var mostCommonChargeIdGroup = idGroup.Identifications.GroupBy(p => p.PrecursorChargeState).MaxBy(p => p.Count());
var id = mostCommonChargeIdGroup.First();

//generate XIC from each file
foreach (var spectraFile in SpectraFileInfoList)
{
var xicIds = mostCommonChargeIdGroup.Where(p => p.FileInfo.Equals(spectraFile)).ToList();
// Try to get the primitive window for the XIC , from the (firstOne - 2) -> (lastOne + 2)
double rightWindow = mostCommonChargeIdGroup.Select(p => p.Ms2RetentionTimeInMinutes).Max() + 2;
double leftWindow = mostCommonChargeIdGroup.Select(p => p.Ms2RetentionTimeInMinutes).Min() - 2;

//If in this run, there is no any ID, we will borrow the ID from other run to build the XIC.
if (xicIds.Count == 0)
{
xicIds = mostCommonChargeIdGroup.Where(p => !p.IsDecoy).ToList();
}
//generate XIC from each file
foreach (var spectraFile in SpectraFileInfoList)
{
var xicIds = mostCommonChargeIdGroup.Where(p => p.FileInfo.Equals(spectraFile)).ToList();

XIC xic = BuildXIC(xicIds, spectraFile, false, leftWindow, rightWindow);
if (xic != null && xic.Ms1Peaks.Count() > 5) // each XIC should have at least 5 points
{
xicGroup.Add(xic);
}
//If in this run, there is no any ID, we will borrow the ID from other run to build the XIC.
if (xicIds.Count == 0)
{
xicIds = mostCommonChargeIdGroup.Where(p => !p.IsDecoy).ToList();
}
Stopwatch xicStopwatch = null;
if (currentProgress == lastReportedProgress)
{
xicStopwatch = Stopwatch.StartNew();
}
XIC xic = BuildXIC(xicIds, spectraFile, false, leftWindow, rightWindow);
if (xicStopwatch != null)
{
xicStopwatch.Stop();
if (!FlashParams.Silent)
{
Console.WriteLine("XIC for {0} took {1} seconds", id.ModifiedSequence, xicStopwatch.Elapsed.TotalSeconds);
}
}
if (xic != null && xic.Ms1Peaks.Count() > 5) // each XIC should have at least 5 points
{
xicGroup.Add(xic);
}
}

// In order to eliminate the bad XIC case that only one id for each file.
// We need to check if the XICGroup has more than one ID in one file.
if ( (!FlashParams.RequireMultipleIdsInOneFiles || MoreThanOneID(xicGroup)) && xicGroup.Count > 1)
{
// Step 1: Find the XIC with most IDs then, set as reference XIC
xicGroup.OrderByDescending(p => p.Ids.Count()).First().Reference = true;
// In order to eliminate the bad XIC case that only one id for each file.
// We need to check if the XICGroup has more than one ID in one file.
if ( (!FlashParams.RequireMultipleIdsInOneFiles || MoreThanOneID(xicGroup)) && xicGroup.Count > 1)
{
// Step 1: Find the XIC with most IDs then, set as reference XIC
xicGroup.OrderByDescending(p => p.Ids.Count()).First().Reference = true;

//Step 2: Build the XICGroups
XICGroups xICGroups = new XICGroups(xicGroup);
Stopwatch xicGroupContstructorStopwatch = null;
if (currentProgress == lastReportedProgress)
{
xicGroupContstructorStopwatch = Stopwatch.StartNew();
}
//Step 2: Build the XICGroups
XICGroups xicGroups = new XICGroups(xicGroup, maxThreads: FlashParams.MaxThreads);
if (xicGroupContstructorStopwatch != null)
{
xicGroupContstructorStopwatch.Stop();
if (!FlashParams.Silent)
{
Console.WriteLine("XICGroup construction for {0}, including peak Alignment, took {1} seconds", id.ModifiedSequence, xicGroupContstructorStopwatch.Elapsed.TotalSeconds);
}
}

//Step 3: Build the peakPairChrom dictionary
//The shared peak dictionary, each sharedPeak included serval ChromatographicPeak for each run [SharedPeak <-> ChromatographicPeaks]
Dictionary<PeakRegion, List<ChromatographicPeak>> sharedPeaksDict = new Dictionary<PeakRegion, List<ChromatographicPeak>>();
//Step 3: Build the peakPairChrom dictionary
//The shared peak dictionary, each sharedPeak included serval ChromatographicPeak for each run [SharedPeak <-> ChromatographicPeaks]
Dictionary<PeakRegion, List<ChromatographicPeak>> sharedPeaksDict = new Dictionary<PeakRegion, List<ChromatographicPeak>>();

//Step 4: Iterate the shared peaks in the XICGroups
foreach (var peak in xICGroups.SharedPeaks)
{
double peakWindow = peak.Width;
if (peakWindow > 0.001) //make sure we have enough length of the window (0.001 min) for peak searching
{
List<ChromatographicPeak> chromPeaksInSharedPeak = new List<ChromatographicPeak>();
CollectChromPeakInRuns(peak, chromPeaksInSharedPeak, xICGroups);
//Step 4: Iterate the shared peaks in the XICGroups
foreach (var peak in xicGroups.SharedPeaks)
{
double peakWindow = peak.Width;
if (peakWindow > 0.001) //make sure we have enough length of the window (0.001 min) for peak searching
{
List<ChromatographicPeak> chromPeaksInSharedPeak = new List<ChromatographicPeak>();
CollectChromPeakInRuns(peak, chromPeaksInSharedPeak, xicGroups);

var allChromPeakInDict = sharedPeaksDict.SelectMany(pair => pair.Value).ToList();
var allChromPeakInDict = sharedPeaksDict.SelectMany(pair => pair.Value).ToList();

// ensuring no duplicates are added and only non-null peaks are considered.
if (chromPeaksInSharedPeak.Where(p => p != null).Any() &&
!chromPeaksInSharedPeak.Any(chromPeak => chromPeak != null && allChromPeakInDict.Any(peakInDict => peakInDict != null && peakInDict.Equals(chromPeak))))
{
sharedPeaksDict[peak] = chromPeaksInSharedPeak;
}
}
// ensuring no duplicates are added and only non-null peaks are considered.
if (chromPeaksInSharedPeak.Where(p => p != null).Any() &&
!chromPeaksInSharedPeak.Any(chromPeak => chromPeak != null && allChromPeakInDict.Any(peakInDict => peakInDict != null && peakInDict.Equals(chromPeak))))
{
sharedPeaksDict[peak] = chromPeaksInSharedPeak;
}
IsobaricPeptideDict.TryAdd(id.ModifiedSequence, sharedPeaksDict);
}
}
IsobaricPeptideDict.TryAdd(id.ModifiedSequence, sharedPeaksDict);
}

// report search progress (proteins searched so far out of total proteins in database)
if (!FlashParams.Silent)
{
Interlocked.Increment(ref isoGroupsSearched);
if (stopwatch != null)
{
stopwatch.Stop();
if (!FlashParams.Silent)
{
Console.WriteLine("IsoTracker for {0} took {1} seconds", id.ModifiedSequence, stopwatch.Elapsed.TotalSeconds);
}
}

double percentProgress = ((double)isoGroupsSearched / PeptideGroupsForIsoTracker.Count * 100);
currentProgress = Math.Max(percentProgress, currentProgress);
// report search progress (proteins searched so far out of total proteins in database)
if (!FlashParams.Silent)
{
Interlocked.Increment(ref isoGroupsSearched);

if (currentProgress > lastReportedProgress + 10)
{
Console.WriteLine("{0:0.}% of isobaric species quantified", currentProgress);
lastReportedProgress = currentProgress;
}
}
double percentProgress = ((double)isoGroupsSearched / PeptideGroupsForIsoTracker.Count * 100);
currentProgress = Math.Max(percentProgress, currentProgress);

if (currentProgress > lastReportedProgress + 5)
{
Console.WriteLine("{0:0.}% of isobaric species quantified", currentProgress);
lastReportedProgress = currentProgress;
}
});
}
}
if (!FlashParams.Silent)
Console.WriteLine("Finished quantifying isobaric species!");
}
Expand All @@ -1938,40 +1972,40 @@ internal bool MoreThanOneID(List<XIC> xics)
internal XIC BuildXIC(List<Identification> ids, SpectraFileInfo spectraFile, bool isReference, double start, double end)
{
Identification id = ids.FirstOrDefault();
var peakIndexingEngine = IndexingEngineDictionary[spectraFile];
var peakIndexingEngine = (PeakIndexingEngine)IndexingEngineDictionary[spectraFile];
PpmTolerance isotopeTolerance = new PpmTolerance(FlashParams.PpmTolerance);
ScanInfo[] ms1ScanInfos = peakIndexingEngine.ScanInfoArray;

ScanInfo startScan = ms1ScanInfos
.Where(p => p.RetentionTime < start)
.OrderBy(p => p.RetentionTime)
.LastOrDefault()
?? ms1ScanInfos.OrderBy(p => p.RetentionTime).First(); // If the start time is before the first scan, use the first scan
// Binary search for start and end scans using a custom comparer for RetentionTime
var scanInfoComparer = Comparer<ScanInfo>.Create((a, b) => a.RetentionTime.CompareTo(b.RetentionTime));

// Create search keys (dummy ScanInfo objects with target retention times)
var startKey = new ScanInfo(0, 0, start, 0);
var endKey = new ScanInfo(0, 0, end, 0);

// Binary search for start scan (scan with RT <= start)
int startIndex = Array.BinarySearch(ms1ScanInfos, startKey, scanInfoComparer);
if (startIndex < 0) startIndex = ~startIndex - 1; // If exact match not found, get largest index with RT <= start
if (startIndex < 0) startIndex = 0; // Handle case where all scans have RT > start

// Binary search for end scan (scan with RT >= end)
int endIndex = Array.BinarySearch(ms1ScanInfos, endKey, scanInfoComparer);
if (endIndex < 0) endIndex = ~endIndex; // If exact match not found, get smallest index with RT >= end
if (endIndex >= ms1ScanInfos.Length) endIndex = ms1ScanInfos.Length - 1; // Handle case where all scans have RT < end

ScanInfo endScan = ms1ScanInfos
.Where(p => p.RetentionTime > end)
.OrderBy(p => p.RetentionTime)
.FirstOrDefault()
?? ms1ScanInfos.OrderBy(p => p.RetentionTime).Last(); // If the end time is after the last scan, use the last scan
ScanInfo startScan = ms1ScanInfos[startIndex];
ScanInfo endScan = ms1ScanInfos[endIndex];

// Collect all peaks from the Ms1 scans in the given time window, then build the XIC
List<IIndexedPeak> peaks = new List<IIndexedPeak>();
for (int j = startScan.ZeroBasedScanIndex; j <= endScan.ZeroBasedScanIndex; j++)
{
double mz = id.PeakfindingMass.ToMz(id.PrecursorChargeState);
IIndexedPeak peak = peakIndexingEngine.GetIndexedPeak(mz , j, isotopeTolerance);
if (peak != null)
{
peaks.Add(peak);
}
}
double mz = id.PeakfindingMass.ToMz(id.PrecursorChargeState);
List<IIndexedPeak> peaks = peakIndexingEngine.GetXicFromScanRange(startScan.ZeroBasedScanIndex, endScan.ZeroBasedScanIndex, mz, isotopeTolerance);

XIC xIC = null;
if (peaks.Count > 0)
{
xIC = new XIC(peaks, id.PeakfindingMass, spectraFile, isReference, ids);
return new XIC(peaks, id.PeakfindingMass, spectraFile, isReference, ids);
}
return xIC;
else
return null;
}

/// <summary>
Expand Down Expand Up @@ -2109,7 +2143,7 @@ internal void AddIsoPeaks()
//remove the repeated peaks from FlashLFQ with the same identification list, the priority is IsoTrack > MSMS
foreach (var peak in allIsoTrackerPeaksInFile)
{
_results.Peaks[fileInfo].RemoveAll(p => IDsEqual(p.Identifications,peak.Identifications));
_results.Peaks[fileInfo].RemoveAll(p => IDsEqual(p.Identifications, peak.Identifications));
}

// Add the peaks into the result dictionary, and remove the duplicated peaks.
Expand Down Expand Up @@ -2139,19 +2173,25 @@ internal void AddIsoPeaks()
/// <param name="idList1"></param>
/// <param name="idList2"></param>
/// <returns></returns>
private bool IDsEqual(List<Identification> idList1, List<Identification> idList2)
public bool IDsEqual(List<Identification> idList1, List<Identification> idList2)
Copy link

Copilot AI Aug 22, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Changing the method visibility from private to public and sorting the input lists in-place can cause side effects for callers. The method modifies the original lists, which could break other code that expects the lists to remain unchanged.

Copilot uses AI. Check for mistakes.
{
if (idList1.Count != idList2.Count)
{
return false;
}

var sortedIdList1 = idList1.OrderBy(id => id.ModifiedSequence).ToList();
var sortedIdList2 = idList2.OrderBy(id => id.ModifiedSequence).ToList();
if (idList1.Count == 1)
{
return idList1[0].ModifiedSequence.Equals(idList2[0].ModifiedSequence);
}

// Sort both lists in place by modified sequence
idList1.Sort((a, b) => a.ModifiedSequence.CompareTo(b.ModifiedSequence));
idList2.Sort((a, b) => a.ModifiedSequence.CompareTo(b.ModifiedSequence));

for (int i = 0; i < sortedIdList1.Count; i++)
for (int i = 0; i < idList1.Count; i++)
{
if (!sortedIdList1[i].ModifiedSequence.Equals(sortedIdList2[i].ModifiedSequence))
if (!idList1[i].ModifiedSequence.Equals(idList2[i].ModifiedSequence))
{
return false;
}
Expand Down
Loading
Loading