private void cameraSearch(final RawDataFile rawFile) {
LOG.finest("Detecting peaks.");
// Get R engine.
final Rengine rEngine;
try {
rEngine = RUtilities.getREngine();
} catch (Throwable t) {
throw new IllegalStateException(
"CAMERA requires R but it couldn't be loaded ("
+ t.getMessage() + ')');
}
synchronized (RUtilities.R_SEMAPHORE) {
// Is R installed - load CAMERA library.
if (rEngine.eval("require(CAMERA)").asBool().isFALSE()) {
throw new IllegalStateException(
"The CAMERA R package couldn't be loaded - is it installed in R?");
}
// Check version of CAMERA.
if (rEngine
.eval("packageVersion('CAMERA') >= '" + CAMERA_VERSION
+ '\'').asBool().isFALSE()) {
throw new IllegalStateException(
"An old version of the CAMERA package is installed in R - please update CAMERA to version "
+ CAMERA_VERSION + " or later");
}
// Create empty peaks matrix.
rEngine.eval(
"columnHeadings <- c('mz','mzmin','mzmax','rt','rtmin','rtmax','into','intb','maxo','sn')",
false);
rEngine.eval(
"peaks <- matrix(nrow=0, ncol=length(columnHeadings))",
false);
rEngine.eval("colnames(peaks) <- columnHeadings", false);
// Initialize.
final ChromatographicPeak[] peaks = peakList.getPeaks(rawFile);
progress = 0.0;
// Initialize scan map.
final Map<Scan, Set<DataPoint>> peakDataPointsByScan = new HashMap<Scan, Set<DataPoint>>(
rawFile.getNumOfScans(MS_LEVEL));
int dataPointCount = 0;
for (final int scanNumber : rawFile.getScanNumbers(MS_LEVEL)) {
// Create a set to hold data points (sorted by m/z).
final Set<DataPoint> dataPoints = new TreeSet<DataPoint>(
ASCENDING_MASS_SORTER);
// Add a dummy data point.
dataPoints.add(new SimpleDataPoint(0.0, 0.0));
dataPointCount++;
// Map the set.
peakDataPointsByScan.put(rawFile.getScan(scanNumber),
dataPoints);
}
// Add peaks.
double progressInc = 1.0 / (double) peaks.length;
for (final ChromatographicPeak peak : peaks) {
// Get peak data.
Range rtRange = null;
Range intRange = null;
final double mz = peak.getMZ();
// Get the peak's data points per scan.
for (final int scanNumber : peak.getScanNumbers()) {
final Scan scan = rawFile.getScan(scanNumber);
if (scan.getMSLevel() != MS_LEVEL) {
throw new IllegalStateException(
"CAMERA can only process peak lists from MS-level "
+ MS_LEVEL);
}
// Copy the data point.
final DataPoint dataPoint = peak.getDataPoint(scanNumber);
if (dataPoint != null) {
final double intensity = dataPoint.getIntensity();
peakDataPointsByScan.get(scan).add(
new SimpleDataPoint(mz, intensity));
dataPointCount++;
// Update RT range.
final double rt = scan.getRetentionTime();
if (rtRange == null) {
rtRange = new Range(rt);
} else {
rtRange.extendRange(rt);
}
// Update intensity range.
if (intRange == null) {
intRange = new Range(intensity);
} else {
intRange.extendRange(intensity);
}
}
}
// Set peak values.
final double area = peak.getArea();
final double maxo = intRange == null
? peak.getHeight()
: intRange.getMax();
final double rtMin = (rtRange == null ? peak
.getRawDataPointsRTRange() : rtRange).getMin();
final double rtMax = (rtRange == null ? peak
.getRawDataPointsRTRange() : rtRange).getMax();
// Add peak row.
rEngine.eval("peaks <- rbind(peaks, c(" + mz + ", " // mz
+ mz + ", " // mzmin: use the same as mz.
+ mz + ", " // mzmax: use the same as mz.
+ peak.getRT() + ", " // rt
+ rtMin + ", " // rtmin
+ rtMax + ", " // rtmax
+ area + ", " // into: peak area.
+ area + ", " // intb: doesn't affect result, use area.
+ maxo + ", " // maxo
+ SIGNAL_TO_NOISE + "))", false);
progress += progressInc;
}
progress = 0.0;
progressInc = 0.25;
// Create R vectors.
final int scanCount = peakDataPointsByScan.size();
final double[] scanTimes = new double[scanCount];
final int[] scanIndices = new int[scanCount];
final double[] masses = new double[dataPointCount];
final double[] intensities = new double[dataPointCount];
// Fill vectors.
int scanIndex = 0;
int pointIndex = 0;
for (final int scanNumber : rawFile.getScanNumbers(MS_LEVEL)) {
final Scan scan = rawFile.getScan(scanNumber);
scanTimes[scanIndex] = scan.getRetentionTime();
scanIndices[scanIndex] = pointIndex + 1;
scanIndex++;
for (final DataPoint dataPoint : peakDataPointsByScan.get(scan)) {
masses[pointIndex] = dataPoint.getMZ();
intensities[pointIndex] = dataPoint.getIntensity();
pointIndex++;
}
}
// Set vectors.
rEngine.assign("scantime", scanTimes);
rEngine.assign("scanindex", scanIndices);
rEngine.assign("mass", masses);
rEngine.assign("intensity", intensities);
// Construct xcmsRaw object
rEngine.eval("xRaw <- new(\"xcmsRaw\")", false);
rEngine.eval("xRaw@tic <- intensity", false);
rEngine.eval("xRaw@scantime <- scantime * " + SECONDS_PER_MINUTE,
false);
rEngine.eval("xRaw@scanindex <- scanindex", false);
rEngine.eval("xRaw@env$mz <- mass", false);
rEngine.eval("xRaw@env$intensity <- intensity", false);
// Create the xcmsSet object.
rEngine.eval("xs <- new('xcmsSet')", false);
// Set peaks.
rEngine.eval("xs@peaks <- peaks", false);
// Set file (dummy) file path.
rEngine.eval("xs@filepaths <- ''", false);
// Set sample name.
rEngine.assign("sampleName", peakList.getName());
rEngine.eval("sampnames(xs) <- sampleName", false);
// Create an empty xsAnnotate.
rEngine.eval("an <- xsAnnotate(xs, sample=1)", false);
// Group by RT.
rEngine.eval("an <- groupFWHM(an, sigma=" + fwhmSigma
+ ", perfwhm=" + fwhmPercentage + ')', false);
progress += progressInc;
// Identify isotopes.
rEngine.eval(
"an <- findIsotopes(an, maxcharge=" + isoMaxCharge
+ ", maxiso=" + isoMaxCount + ", ppm="
+ isoMassTolerance.getPpmTolerance() + ", mzabs="
+ isoMassTolerance.getMzTolerance() + ')', false);
progress += progressInc;
// Split groups by correlating peak shape (need to set xraw to raw
// data).
rEngine.eval(
"an <- groupCorr(an, calcIso=TRUE, xraw=xRaw, cor_eic_th="
+ corrThreshold + ", pval=" + corrPValue + ')',
false);
progress += progressInc;
// Get the peak list.
rEngine.eval("peakList <- getPeaklist(an)", false);
// Extract the pseudo-spectra and isotope annotations from the peak
// list.
final REXP spectraExp = rEngine
.eval("as.integer(peakList$pcgroup)");
final REXP isotopeExp = rEngine.eval("peakList$isotopes");
// Add identities.
if (spectraExp != null) {
addPseudoSpectraIdentities(peaks, spectraExp, isotopeExp);