Skip to content
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.

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

Fixed Gaussian curve Peak Transfer overfitting #574

Merged
merged 1 commit into from
Dec 3, 2024
Merged
Changes from all commits
Commits
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
Original file line number Diff line number Diff line change
Expand Up @@ -8,6 +8,7 @@
*
* Contributors:
* Philip Wenig - initial API and implementation
* Matthias Mailänder - optimize sigma estimation
*******************************************************************************/
package net.openchrom.xxd.process.supplier.templates.peaks;

Expand Down Expand Up @@ -120,31 +121,31 @@ private IProcessingInfo<R> applyDetector(IChromatogramSelection<? extends IPeak,
IProcessingInfo<R> processingInfo = super.validate(chromatogramSelection, settings, monitor);
if(!processingInfo.hasErrorMessages()) {
if(settings instanceof PeakTransferSettings peakTransferSettings) {
transferPeaks(chromatogramSelection, peakTransferSettings);
transferPeaks(chromatogramSelection, peakTransferSettings, monitor);
} else {
processingInfo.addErrorMessage(PeakDetectorSettings.DETECTOR_DESCRIPTION, "The settings instance is wrong.");
}
}
return processingInfo;
}

private void transferPeaks(IChromatogramSelection<? extends IPeak, ?> chromatogramSelection, PeakTransferSettings peakTransferSettings) {
private void transferPeaks(IChromatogramSelection<? extends IPeak, ?> chromatogramSelection, PeakTransferSettings peakTransferSettings, IProgressMonitor monitor) {

IChromatogram<?> chromatogram = chromatogramSelection.getChromatogram();
List<? extends IPeak> peaks = chromatogram.getPeaks(chromatogramSelection);
List<IChromatogram<?>> referencedChromatograms = chromatogram.getReferencedChromatograms();
for(IChromatogram<?> referencedChromatogram : referencedChromatograms) {
transferPeaks(peaks, referencedChromatogram, peakTransferSettings);
transferPeaks(peaks, referencedChromatogram, peakTransferSettings, monitor);
}
}

private void transferPeaks(List<? extends IPeak> peaks, IChromatogram<?> chromatogramSink, PeakTransferSettings peakTransferSettings) {
private void transferPeaks(List<? extends IPeak> peaks, IChromatogram<?> chromatogramSink, PeakTransferSettings peakTransferSettings, IProgressMonitor monitor) {

Map<Integer, List<IPeak>> peakGroups = extractPeakGroups(peaks, peakTransferSettings);
List<Integer> groups = new ArrayList<>();
groups.addAll(peakGroups.keySet());
Collections.sort(groups);
//
monitor.beginTask("Transfer Peaks", peaks.size());
for(int group : groups) {
List<IPeak> peakz = peakGroups.get(group);
if(peakz.size() == 1) {
Expand All @@ -154,6 +155,7 @@ private void transferPeaks(List<? extends IPeak> peaks, IChromatogram<?> chromat
IPeak peak = peakz.get(0);
double percentageIntensity = getPercentageIntensity(peak);
transfer(peak, percentageIntensity, chromatogramSink, peakTransferSettings);
monitor.worked(1);
} else {
/*
* Peak Group
Expand All @@ -164,6 +166,7 @@ private void transferPeaks(List<? extends IPeak> peaks, IChromatogram<?> chromat
} catch(Exception e) {
logger.warn(e);
}
monitor.worked(1);
}
}
}
Expand Down Expand Up @@ -272,7 +275,7 @@ private void transfer(IPeak peakSource, int startRetentionTime, int stopRetentio
//
IPeak peakSink = peakSupport.extractPeakByRetentionTime(chromatogramSink, startRetentionTime, stopRetentionTime, includeBackground, optimizeRange, traces);
if(peakSink != null) {
adjustPeakIntensity(peakSource, peakSink, percentageIntensity, peakTransferSettings);
adjustPeakIntensity(peakSink, percentageIntensity, peakTransferSettings);
transferTargets(peakSource, peakSink, peakTransferSettings);
addPeak(chromatogramSink, peakSink);
}
Expand Down Expand Up @@ -306,19 +309,7 @@ private void transfer(IPeak peakSource, IChromatogram<? extends IPeak> chromatog
Point maxPosition = getMaxPosition(chromatogramCSD, peakModel.getRetentionTimeAtPeakMaximum(), offsetRetentionTime);
//
if(maxPosition.getX() > 0 && maxPosition.getY() > 0) {
/*
* Parameter
*/
double sigma = calculateSigma(chromatogramCSD, peakSource);
int centerRetentionTime = (int)maxPosition.getX();
float intensity = (float)(maxPosition.getY() * percentageIntensity);
int centerScan = chromatogramCSD.getScanNumber(centerRetentionTime);
double[] parameterStandard = new double[3];
parameterStandard[0] = intensity; // Norm
parameterStandard[1] = centerScan; // Mean
parameterStandard[2] = sigma; // Sigma
IPeak peakSink = createDefaultGaussPeakNormal(chromatogramCSD, startRetentionTime, centerRetentionTime, stopRetentionTime, parameterStandard);
//
IPeak peakSink = fitPeak(chromatogramCSD, maxPosition, percentageIntensity, startRetentionTime, stopRetentionTime);
if(peakSink != null) {
transferTargets(peakSource, peakSink, peakTransferSettings);
addPeak(chromatogramSink, peakSink);
Expand All @@ -333,7 +324,41 @@ private void transfer(IPeak peakSource, IChromatogram<? extends IPeak> chromatog
}
}

private double calculateSigma(IChromatogramCSD chromatogramCSD, IPeak peak) {
private IPeak fitPeak(IChromatogramCSD chromatogramCSD, Point maxPosition, double percentageIntensity, int startRetentionTime, int stopRetentionTime) {

int centerRetentionTime = (int)maxPosition.getX();
float intensity = (float)(maxPosition.getY() * percentageIntensity);
int centerScan = chromatogramCSD.getScanNumber(centerRetentionTime);
double sigma = calculateSigma(chromatogramCSD);
Gaussian gaussian = new Gaussian(intensity, centerScan, sigma);
IPeak peakSink = createDefaultGaussPeakNormal(chromatogramCSD, startRetentionTime, stopRetentionTime, gaussian, intensity);
if(peakSink == null) {
return null;
}
while(!isWithinBounds(peakSink, chromatogramCSD)) {
sigma = sigma / 2;
gaussian = new Gaussian(intensity, centerScan, sigma);
peakSink = createDefaultGaussPeakNormal(chromatogramCSD, startRetentionTime, stopRetentionTime, gaussian, intensity);
if(peakSink == null) {
return null;
}
}
return peakSink;
}

private boolean isWithinBounds(IPeak peakSink, IChromatogramCSD chromatogramCSD) {

IPeakModel peakModel = peakSink.getPeakModel();
for(int retentionTime : peakModel.getRetentionTimes()) {
int scanNumber = chromatogramCSD.getScanNumber(retentionTime);
if(peakModel.getPeakAbundance(retentionTime) > chromatogramCSD.getScan(scanNumber).getTotalSignal()) {
return false;
}
}
return true;
}

private double calculateSigma(IChromatogramCSD chromatogramCSD) {

int scanInterval = chromatogramCSD.getScanInterval();
return 1000.0d / scanInterval;
Expand All @@ -344,6 +369,7 @@ private Set<Integer> getTraces(IPeak peakSource, int numberTraces) {
Set<Integer> traces = new HashSet<>();
if(peakSource instanceof IChromatogramPeakMSD peakMSD) {
if(peakMSD.getPurity() < 1.0f && numberTraces > 0) {
// probably a deconvoluted peak
IScanMSD scanMSD = peakMSD.getExtractedMassSpectrum();
if(scanMSD.getIons().size() <= numberTraces) {
IExtractedIonSignal extractedIonSignal = peakMSD.getExtractedMassSpectrum().getExtractedIonSignal();
Expand All @@ -358,7 +384,7 @@ private Set<Integer> getTraces(IPeak peakSource, int numberTraces) {
return traces;
}

private void adjustPeakIntensity(IPeak peakSource, IPeak peakSink, double percentageIntensity, PeakTransferSettings peakTransferSettings) {
private void adjustPeakIntensity(IPeak peakSink, double percentageIntensity, PeakTransferSettings peakTransferSettings) {

if(peakTransferSettings.isAdjustPeakHeight()) {
if(percentageIntensity > 0.0d && percentageIntensity < 1.0d) {
Expand Down Expand Up @@ -436,20 +462,15 @@ private Point getMaxPosition(IChromatogram<?> chromatogram, int centerRetentionT
return new Point(retentionTime, maxIntensity);
}

public IChromatogramPeakCSD createDefaultGaussPeakNormal(IChromatogramCSD chromatogram, int startRetentionTime, int centerRetentionTime, int stopRetentionTime, double[] parameter) {
public IChromatogramPeakCSD createDefaultGaussPeakNormal(IChromatogramCSD chromatogram, int startRetentionTime, int stopRetentionTime, Gaussian gaussian, float norm) {

int startScan = chromatogram.getScanNumber(startRetentionTime);
int stopScan = chromatogram.getScanNumber(stopRetentionTime);
/*
* Intensity profile
*/
IScanCSD peakMaximum = new ScanCSD((float)parameter[0]);
IPeakIntensityValues peakIntensityValues = new PeakIntensityValues((float)parameter[0]);
//
double norm = parameter[0];
double mean = parameter[1];
double sigma = parameter[2];
Gaussian gaussian = new Gaussian(norm, mean, sigma);
IScanCSD peakMaximum = new ScanCSD(norm);
IPeakIntensityValues peakIntensityValues = new PeakIntensityValues(norm);
//
for(int i = startScan; i <= stopScan; i++) {
IScan scan = chromatogram.getScan(i);
Expand Down
Loading