/*************************************************************************
* *
* 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.Map;
import java.util.HashMap;
import java.util.Arrays;
import java.util.ArrayList;
import java.util.Iterator;
import java.util.Collections;
import java.io.PrintStream;
import java.io.OutputStream;
import java.io.FileOutputStream;
import org.apache.commons.lang3.tuple.Pair;
import java.io.File;
public class AnimateNetCDFAroundMass {
class XYZ {
Double time;
Double mz;
Double intensity;
public XYZ(Double time, Double mz, Double intensity) {
this.time = time;
this.mz = mz;
this.intensity = intensity;
}
}
private List<List<XYZ>> getSpectraInWindowAll(List<List<XYZ>> spectra, Double time, Double tWin,
Double mz, Double mzWin) {
List<List<XYZ>> windowOnMultiple = new ArrayList<>();
for (List<XYZ> one : spectra)
windowOnMultiple.add(getSpectraInWindow(one, time, tWin, mz, mzWin));
return windowOnMultiple;
}
final static int gridSize = 100;
private List<XYZ> getSpectraInWindow(List<XYZ> spectra, Double time, Double tWin, Double mz, Double mzWin) {
// tWin and mzWin are half windows around time and mz, respectively
// Range of mz is ( mz - mzWin, mz + mzWin )
// Range of time is ( time - timeWin, time + timeWin )
double mzLow = mz - mzWin;
double mzHigh = mz + mzWin;
double timeLow = time - tWin;
double timeHigh = time + tWin;
Map<Pair<Long, Long>, List<Double>> gridVals = new HashMap<>();
// initialize the grid to empty intensity data points
for (long i = 0; i <= gridSize; i++) {
for (long j = 0; j <= gridSize; j++) {
gridVals.put(Pair.of(i, j), new ArrayList<>());
}
}
// what is the step size between x -> x+1 and y -> y+1
double timeStep = (2 * tWin) / gridSize;
double mzStep = (2 * mzWin) / gridSize;
// for each x,y,z find the grid position it lies in
// (if it indeed is within the grid window we care)
// and add that to the data points we see within
for (XYZ xyz : spectra) {
// ignore if outside the window
if (xyz.mz < mzLow || xyz.mz > mzHigh)
continue;
if (xyz.time < timeLow || xyz.time > timeHigh)
continue;
// if inside the window find the grid position
long timeGridPos = Math.round((xyz.time - timeLow) / timeStep);
long mzGridPos = Math.round((xyz.mz - mzLow) / mzStep);
Pair<Long, Long> xy = Pair.of(timeGridPos, mzGridPos);
gridVals.get(xy).add(xyz.intensity);
}
// average out all the z values that appear at each x,y
Map<Pair<Long, Long>, Double> gridAvg = new HashMap<>();
for (Pair<Long, Long> xy : gridVals.keySet()) {
Double avg = 0.0;
List<Double> intensities = gridVals.get(xy);
for (Double intensity : intensities)
avg += intensity;
if (!intensities.isEmpty())
avg /= intensities.size();
gridAvg.put(xy, avg);
}
// convert back to list of x,y,z coordinates
List<XYZ> grid = new ArrayList<>();
for (long i = 0; i <= gridSize; i++) {
for (long j = 0; j <= gridSize; j++) {
Double t = timeLow + i * timeStep;
Double m = mzLow + j * mzStep;
Double it = gridAvg.get(Pair.of(i, j));
grid.add(new XYZ(t, m, it));
}
}
return grid;
}
private List<XYZ> getSpectra(Iterator<LCMSSpectrum> spectraIt, Double time, Double tWin, Double mz,
Double mzWin) {
double mzLow = mz - mzWin;
double mzHigh = mz + mzWin;
double timeLow = time - tWin;
double timeHigh = time + tWin;
List<XYZ> spectra = new ArrayList<>();
while (spectraIt.hasNext()) {
LCMSSpectrum timepoint = spectraIt.next();
if (timepoint.getTimeVal() < timeLow || timepoint.getTimeVal() > timeHigh)
continue;
// get all (mz, intensity) at this timepoint
for (Pair<Double, Double> mz_int : timepoint.getIntensities()) {
double mzHere = mz_int.getLeft();
double intensity = mz_int.getRight();
if (mzHere < mzLow || mzHere > mzHigh)
continue;
spectra.add(new XYZ(timepoint.getTimeVal(), mzHere, intensity));
}
}
return spectra;
}
public List<List<XYZ>> getSpectra(String[] fnames, Double time, Double tWin, Double mz,
Double mzWin) throws Exception {
List<List<XYZ>> extracted = new ArrayList<>();
LCMSParser parser = new LCMSNetCDFParser();
for (String fname : fnames) {
Iterator<LCMSSpectrum> iter = parser.getIterator(fname);
extracted.add(getSpectra(iter, time, tWin, mz, mzWin));
}
return extracted;
}
private static boolean areNCFiles(String[] fnames) {
for (String n : fnames) {
System.out.println(".nc file = " + n);
if (!n.endsWith(".nc"))
return false;
}
return true;
}
public static void main(String[] args) throws Exception {
if (args.length < 7 || !areNCFiles(Arrays.copyOfRange(args, 5, args.length))) {
throw new RuntimeException("Needs: \n" +
"(1) mass value, e.g., 132.0772 \n" +
"(2) time value, e.g., 39.2, (seconds), \n" +
"(3) minimum Mz Precision, 0.04 \n" +
"(4) max z axis, e.g., 20000 \n" +
"(5) prefix for .data and rendered .pdf \n" +
"(6..) 2 or more NetCDF .nc files"
);
}
Double mz = Double.parseDouble(args[0]);
Double time = Double.parseDouble(args[1]);
Double minMzPrecision = Double.parseDouble(args[2]);
Double maxZAxis = Double.parseDouble(args[3]);
String outPrefix = args[4];
// the mz values go from 50-950, we start with a big window and exponentially narrow down
double mzWin = 100;
// time values go from 0-450, we start with a big window and exponentially narrow down
double timeWin = 50;
// the factor by which to zoom in every step (has to be >1, a value of 2 is good)
double factor = 1.2;
// the animation frame count
int frame = 1;
AnimateNetCDFAroundMass c = new AnimateNetCDFAroundMass();
String[] netCDFFnames = Arrays.copyOfRange(args, 5, args.length);
List<List<XYZ>> spectra = c.getSpectra(netCDFFnames, time, timeWin, mz, mzWin);
for (List<XYZ> s : spectra) {
System.out.format("%d xyz datapoints in (initial narrowed) spectra\n", s.size());
}
String[] labels = new String[netCDFFnames.length];
for (int i=0; i<labels.length; i++)
labels[i] = "Dataset: " + i;
// you could set labels to netCDFFnames to get precise labels on the graphs
Gnuplotter plotter = new Gnuplotter();
String fmt = "png";
List<String> outImgFiles = new ArrayList<>(), outDataFiles = new ArrayList<>();
while (mzWin > minMzPrecision) {
// exponentially narrow windows down
mzWin /= factor;
timeWin /= factor;
List<List<XYZ>> windowedSpectra = c.getSpectraInWindowAll(spectra, time, timeWin, mz, mzWin);
String frameid = String.format("%03d", frame);
String outPDF = outPrefix + frameid + "." + fmt;
String outDATA = outPrefix + frameid + ".data";
outImgFiles.add(outPDF);
outDataFiles.add(outDATA);
frame++;
// Write data output to outfile
PrintStream out = new PrintStream(new FileOutputStream(outDATA));
// print out the spectra to outDATA
for (List<XYZ> windowOfSpectra : windowedSpectra) {
for (XYZ xyz : windowOfSpectra) {
out.format("%.4f\t%.4f\t%.4f\n", xyz.time, xyz.mz, xyz.intensity);
out.flush();
}
// delimit this dataset from the rest
out.print("\n\n");
}
// close the .data
out.close();
// render outDATA to outPDF using gnuplot
plotter.plotMulti3D(outDATA, outPDF, fmt, labels, maxZAxis);
}
String outImgs = outPrefix + "*." + fmt;
plotter.makeAnimatedGIF(outImgs, outPrefix + ".gif");
// all the frames are now in the animated gif, remove the intermediate files
for (String f: outDataFiles)
new File(f).delete();
for (String f: outImgFiles)
new File(f).delete();
}
}