Java tutorial
/************************************************************************* * * * This file is part of the 20n/act project. * * 20n/act enables DNA prediction for synthetic biology/bioengineering. * * Copyright (C) 2017 20n Labs, Inc. * * * * Please direct all queries to act@20n.com. * * * * This program is free software: you can redistribute it and/or modify * * it under the terms of the GNU General Public License as published by * * the Free Software Foundation, either version 3 of the License, or * * (at your option) any later version. * * * * This program is distributed in the hope that it will be useful, * * but WITHOUT ANY WARRANTY; without even the implied warranty of * * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * * GNU General Public License for more details. * * * * You should have received a copy of the GNU General Public License * * along with this program. If not, see <http://www.gnu.org/licenses/>. * * * *************************************************************************/ package com.act.lcms; import java.util.List; import java.util.ArrayList; import java.util.Iterator; import org.apache.commons.lang3.tuple.Pair; import org.apache.commons.lang3.tuple.Triple; import java.io.PrintStream; import java.io.OutputStream; import java.io.FileOutputStream; public class ExtractFromNetCDFAroundMass { /** * Extracts a window around a particular m/z value within the spectra * @param file The netCDF file with the full LCMS (lockmass corrected) spectra * @param mz The m/z value +- 1 around which to extract data for * @param numTimepointsToExamine Since we stream the netCDF in, we can choose * to look at the first few timepoints or -1 for all */ public List<Triple<Double, Double, Double>> get2DWindow(String file, double mz, int numTimepointsToExamine) throws Exception { LCMSParser parser = new LCMSNetCDFParser(); Iterator<LCMSSpectrum> iter = parser.getIterator(file); List<Triple<Double, Double, Double>> window = new ArrayList<>(); int pulled = 0; Double mzTightL = mz - 0.1; Double mzTightR = mz + 0.1; Double mzMinus1Da = mz - 1; Double mzPlus1Da = mz + 1; // iterate through first few, or all if -1 given while (iter.hasNext() && (numTimepointsToExamine == -1 || pulled < numTimepointsToExamine)) { LCMSSpectrum timepoint = iter.next(); List<Pair<Double, Double>> intensities = timepoint.getIntensities(); // this time point is valid to look at if its max intensity is around // the mass we care about. So lets first get the max peak location Pair<Double, Double> max_peak = findMaxPeak(intensities); Double max_mz = max_peak.getLeft(); // If the max_mz value is pretty close to our target mass, ie in [mzTightL, mzTightR] // Then this timepoint is a good timepoint to output... proceed, shall we. if (max_mz >= mzTightL && max_mz <= mzTightR) { // For this timepoint, output a window for (int k = 0; k < intensities.size(); k++) { Double mz_here = intensities.get(k).getLeft(); Double intensity = intensities.get(k).getRight(); // The window not as tight, but +-1 Da around the target mass if (mz_here > mzMinus1Da && mz_here < mzPlus1Da) { window.add(Triple.of(timepoint.getTimeVal(), mz_here, intensity)); } } } pulled++; } return window; } private Pair<Double, Double> findMaxPeak(List<Pair<Double, Double>> raw) { // the intensity is the second value in the pairs // this finds the pair with the max intensity Pair<Double, Double> max = null; for (Pair<Double, Double> mz_int : raw) { if (max == null || max.getRight() < mz_int.getRight()) max = mz_int; } return max; } public static void main(String[] args) throws Exception { if (args.length != 4 || !args[0].endsWith(".nc")) { throw new RuntimeException( "Needs (1) NetCDF .nc file, " + "(2) mass value, e.g., 132.0772 for debugging, " + "(3) how many timepoints to process (-1 for all), " + "(4) prefix for .data and rendered .pdf, '-' for stdout"); } String netCDF = args[0]; Double mz = Double.parseDouble(args[1]); Integer numSpectraToProcess = Integer.parseInt(args[2]); String outPrefix = args[3]; String outPDF = outPrefix.equals("-") ? null : outPrefix + ".pdf"; String outDATA = outPrefix.equals("-") ? null : outPrefix + ".data"; ExtractFromNetCDFAroundMass e = new ExtractFromNetCDFAroundMass(); List<Triple<Double, Double, Double>> window = e.get2DWindow(netCDF, mz, numSpectraToProcess); // Write data output to outfile PrintStream whereTo = outDATA == null ? System.out : new PrintStream(new FileOutputStream(outDATA)); for (Triple<Double, Double, Double> xyz : window) { whereTo.format("%.4f\t%.4f\t%.4f\n", xyz.getLeft(), xyz.getMiddle(), xyz.getRight()); whereTo.flush(); } if (outDATA != null) { // if outDATA is != null, then we have written to .data file // now render the .data to the corresponding .pdf file // first close the .data whereTo.close(); // render outDATA to outPDF using gnuplo Gnuplotter plotter = new Gnuplotter(); plotter.plot3D(outDATA, outPDF, netCDF, mz); } } }