/*************************************************************************
* *
* 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.plotter;
import com.act.utils.TSVWriter;
import com.act.lcms.Gnuplotter;
import com.act.lcms.MS1;
import com.act.lcms.XZ;
import com.act.lcms.db.analysis.PathwayProductAnalysis;
import com.act.lcms.db.analysis.ScanData;
import com.act.lcms.db.model.LCMSWell;
import com.act.lcms.db.model.MS1ScanForWellAndMassCharge;
import org.apache.commons.lang3.tuple.Pair;
import java.io.FileOutputStream;
import java.io.IOException;
import java.io.OutputStream;
import java.io.PrintStream;
import java.util.ArrayList;
import java.util.HashMap;
import java.util.List;
import java.util.Map;
import java.util.Set;
public class WriteAndPlotMS1Results {
public void plotTIC(List<XZ> tic, String outPrefix, String fmt) throws IOException {
String outImg = outPrefix + "." + fmt;
String outData = outPrefix + ".data";
// Write data output to outfile
try (PrintStream out = new PrintStream(new FileOutputStream(outData))) {
// print each time point + intensity to outDATA
for (XZ xz : tic) {
out.format("%.4f\t%.4f\n", xz.getTime(), xz.getIntensity());
out.flush();
}
}
// render outDATA to outPDF using gnuplot
new Gnuplotter().plot2D(outData, outImg, new String[] { "TIC" }, "time", null, "intensity",
fmt);
}
public void plotScan(List<MS1.YZ> scan, String outPrefix, String fmt) throws IOException {
String outPDF = outPrefix + "." + fmt;
String outDATA = outPrefix + ".data";
// Write data output to outfile
PrintStream out = new PrintStream(new FileOutputStream(outDATA));
// print out the spectra to outDATA
for (MS1.YZ yz : scan) {
out.format("%.4f\t%.4f\n", yz.getMZ(), yz.getIntensity());
out.flush();
}
// close the .data
out.close();
// render outDATA to outPDF using gnuplot
new Gnuplotter().plot2DImpulsesWithLabels(outDATA, outPDF, new String[] { "mz distribution at TIC max" },
null, "mz", null, "intensity", fmt);
}
private List<String> writeFeedMS1Values(List<Pair<Double, List<XZ>>> ms1s, Double maxIntensity,
OutputStream os) throws IOException {
// Write data output to outfile
PrintStream out = new PrintStream(os);
List<String> plotID = new ArrayList<>(ms1s.size());
for (Pair<Double, List<XZ>> ms1ForFeed : ms1s) {
Double feedingConcentration = ms1ForFeed.getLeft();
List<XZ> ms1 = ms1ForFeed.getRight();
plotID.add(String.format("concentration: %5e", feedingConcentration));
// print out the spectra to outDATA
for (XZ xz : ms1) {
out.format("%.4f\t%.4f\n", xz.getTime(), xz.getIntensity());
out.flush();
}
// delimit this dataset from the rest
out.print("\n\n");
}
return plotID;
}
private void writeFeedMS1Values(List<Pair<Double, Double>> concentrationIntensity, OutputStream os)
throws IOException {
PrintStream out = new PrintStream(os);
for (Pair<Double, Double> ci : concentrationIntensity) {
out.format("%f\t%f\n", ci.getLeft(), ci.getRight());
}
out.flush();
}
// input: list sorted on first field of pair of (concentration, ms1 spectra)
// the ion of relevance to compare across different spectra
// outPrefix for pdfs and data, and fmt (pdf or png) of output
public void plotFeedings(List<Pair<Double, MS1ScanForWellAndMassCharge>> feedings, String ion, String outPrefix,
String fmt, String gnuplotFile)
throws IOException {
String outSpectraImg = outPrefix + "." + fmt;
String outSpectraData = outPrefix + ".data";
String outFeedingImg = outPrefix + ".fed." + fmt;
String outFeedingData = outPrefix + ".fed.data";
String feedingGnuplotFile = gnuplotFile + ".fed";
boolean useMaxPeak = true;
// maps that hold the values for across different concentrations
List<Pair<Double, List<XZ>>> concSpectra = new ArrayList<>();
List<Pair<Double, Double>> concAreaUnderSpectra = new ArrayList<>();
List<Pair<Double, Double>> concMaxPeak = new ArrayList<>();
// we will compute a running max of the intensity in the plot, and integral
Double maxIntensity = 0.0d, maxAreaUnder = 0.0d;
// now compute the maps { conc -> spectra } and { conc -> area under spectra }
for (Pair<Double, MS1ScanForWellAndMassCharge> feedExpr : feedings) {
Double concentration = feedExpr.getLeft();
MS1ScanForWellAndMassCharge scan = feedExpr.getRight();
// get the ms1 spectra for the selected ion, and the max for it as well
List<XZ> ms1 = scan.getIonsToSpectra().get(ion);
Double maxInThisSpectra = scan.getMaxIntensityForIon(ion);
Double areaUnderSpectra = scan.getIntegralForIon(ion);
// update the max intensity over all different spectra
maxIntensity = Math.max(maxIntensity, maxInThisSpectra);
maxAreaUnder = Math.max(maxAreaUnder, areaUnderSpectra);
// install this concentration and spectra in map, to be dumped to file later
concSpectra.add(Pair.of(concentration, ms1));
concAreaUnderSpectra.add(Pair.of(concentration, areaUnderSpectra));
concMaxPeak.add(Pair.of(concentration, maxInThisSpectra));
}
// Write data output to outfiles
List<String> plotID = null;
try (FileOutputStream outSpectra = new FileOutputStream(outSpectraData)) {
plotID = writeFeedMS1Values(concSpectra, maxIntensity, outSpectra);
}
try (FileOutputStream outFeeding = new FileOutputStream(outFeedingData)) {
writeFeedMS1Values(useMaxPeak ? concMaxPeak : concAreaUnderSpectra, outFeeding);
}
// render outDATA to outPDF using gnuplot
Gnuplotter gp = new Gnuplotter();
String[] plotNames = plotID.toArray(new String[plotID.size()]);
gp.plotOverlayed2D(outSpectraData, outSpectraImg, plotNames, "time", maxIntensity, "intensity", fmt, gnuplotFile);
gp.plot2D(outFeedingData, outFeedingImg, new String[] { "feeding ramp" }, "concentration",
useMaxPeak ? maxIntensity : maxAreaUnder, "integrated area under spectra", fmt, null, null, null,
feedingGnuplotFile);
}
private List<Pair<String, String>> writeMS1Values(Map<String, List<XZ>> ms1s, Double maxIntensity,
Map<String, Double> metlinMzs, OutputStream os,
boolean heatmap) throws IOException {
return writeMS1Values(ms1s, maxIntensity, metlinMzs, os, heatmap, true, null);
}
public List<Pair<String, String>> writeMS1Values(Map<String, List<XZ>> ms1s, Double maxIntensity,
Map<String, Double> metlinMzs, OutputStream os, boolean heatmap,
boolean applyThreshold, Set<String> ionsToWrite) throws IOException {
// Write data output to outfile
PrintStream out = new PrintStream(os);
List<Pair<String, String>> plotID = new ArrayList<>(ms1s.size());
for (Map.Entry<String, List<XZ>> ms1ForIon : ms1s.entrySet()) {
String ion = ms1ForIon.getKey();
// Skip ions not in the ionsToWrite set if that set is defined.
if (ionsToWrite != null && !ionsToWrite.contains(ion)) {
continue;
}
List<XZ> ms1 = ms1ForIon.getValue();
String plotName = String.format("ion: %s, mz: %.5f", ion, metlinMzs.get(ion));
plotID.add(Pair.of(ion, plotName));
// print out the spectra to outDATA
for (XZ xz : ms1) {
if (heatmap) {
/*
* When we are building heatmaps, we use gnuplots pm3d package
* along with `dgrid3d 2000,2` (which averages data into grids
* that are 2000 on the time axis and 2 in the y axis), and
* `view map` that flattens a 3D graphs into a 2D view.
* We want time to be on the x-axis and intensity on the z-axis
* (because that is the one that is mapped to heat colors)
* but then we need an artificial y-axis. We create proxy y=1
* and y=2 datapoints, and then dgrid3d averaging over 2 creates
* a vertical "strip".
*/
out.format("%.4f\t1\t%.4f\n", xz.getTime(), xz.getIntensity());
out.format("%.4f\t2\t%.4f\n", xz.getTime(), xz.getIntensity());
} else {
out.format("%.4f\t%.4f\n", xz.getTime(), xz.getIntensity());
}
out.flush();
}
// delimit this dataset from the rest
out.print("\n\n");
}
return plotID;
}
public void plotSpectra(Map<String, List<XZ>> ms1s, Double maxIntensity,
Map<String, Double> individualMaxIntensities, Map<String, Double> metlinMzs,
String outPrefix, String fmt, boolean makeHeatmap, boolean overlayPlots)
throws IOException {
String outImg = outPrefix + "." + fmt;
String outData = outPrefix + ".data";
// Write data output to outfile
try (FileOutputStream out = new FileOutputStream(outData)) {
List<Pair<String, String>> ionAndplotID = writeMS1Values(ms1s, maxIntensity, metlinMzs, out, makeHeatmap);
// writeMS1Values picks an ordering of the plots.
// create two new sets plotID and yMaxes that have the matching ordering
// and contain plotNames, and yRanges respectively
List<Double> yMaxesInSameOrderAsPlots = new ArrayList<>();
List<String> plotID = new ArrayList<>();
for (Pair<String, String> plot : ionAndplotID) {
String ion = plot.getLeft();
Double yMax = individualMaxIntensities.get(ion);
yMaxesInSameOrderAsPlots.add(yMax);
plotID.add(plot.getRight());
}
Double[] yMaxes = yMaxesInSameOrderAsPlots.toArray(new Double[yMaxesInSameOrderAsPlots.size()]);
// render outDATA to outPDF using gnuplot
Gnuplotter gp = new Gnuplotter();
String[] plotNames = plotID.toArray(new String[plotID.size()]);
if (makeHeatmap) {
gp.plotHeatmap(outData, outImg, plotNames, maxIntensity, fmt);
} else {
if (!overlayPlots) {
gp.plot2D(outData, outImg, plotNames, "time", maxIntensity, "intensity", fmt,
null, null, yMaxes, outImg + ".gnuplot");
} else {
gp.plotOverlayed2D(outData, outImg, plotNames, "time", maxIntensity, "intensity", fmt, outImg + ".gnuplot");
}
}
}
}
/**
* This function writes the pathway product results to an analysis file
* @param writer This is a handle to the tsv writer
* @param chemicalName The name of the primary chemical of inspection
* @param positiveAndNegativeWells This is a list of positive and negative control well samples
* @param pathwayStepIon This is the metlin ion which usually is the best metlin ion based on standard ion analysis
* @param wellsToBestPeaks This is a map of well to the best peak positions in the spectral charts
* @throws IOException
*/
public static void writePathwayProductOutput(TSVWriter<String, String> writer, String chemicalName,
List<ScanData<LCMSWell>> positiveAndNegativeWells,
String pathwayStepIon,
Map<ScanData<LCMSWell>, XZ> wellsToBestPeaks) throws IOException {
for (ScanData<LCMSWell> well : positiveAndNegativeWells) {
String fedChemical = well.getWell().getChemical() == null ||
well.getWell().getChemical().isEmpty() ? "nothing" : well.getWell().getChemical();
String pelletOrSupernatant = well.getPlate().getDescription().contains("pellet") ? "pellet" : "supernatant";
String detected;
String intensity;
String time;
if (wellsToBestPeaks.get(well) != null) {
detected = "YES";
intensity = String.format("%.4f", wellsToBestPeaks.get(well).getIntensity());
time = String.format("%.4f", wellsToBestPeaks.get(well).getTime());
} else {
detected = "NO";
intensity = "-";
time = "-";
}
Map<String, String> row = new HashMap<>();
row.put(PathwayProductAnalysis.PATHWAY_PRODUCT_HEADER_FIELDS.TARGET_CHEMICAL.name(), chemicalName);
row.put(PathwayProductAnalysis.PATHWAY_PRODUCT_HEADER_FIELDS.TYPE.name(), pelletOrSupernatant);
row.put(PathwayProductAnalysis.PATHWAY_PRODUCT_HEADER_FIELDS.FED_CHEMICAL.name(), fedChemical);
row.put(PathwayProductAnalysis.PATHWAY_PRODUCT_HEADER_FIELDS.DETECTED.name(), detected);
row.put(PathwayProductAnalysis.PATHWAY_PRODUCT_HEADER_FIELDS.INTENSITY.name(), intensity);
row.put(PathwayProductAnalysis.PATHWAY_PRODUCT_HEADER_FIELDS.TIME.name(), time);
row.put(PathwayProductAnalysis.PATHWAY_PRODUCT_HEADER_FIELDS.PLATE_BARCODE.name(), well.getPlate().getBarcode());
row.put(PathwayProductAnalysis.PATHWAY_PRODUCT_HEADER_FIELDS.MODE.name(), well.getScanFile().getMode().toString().toLowerCase());
row.put(PathwayProductAnalysis.PATHWAY_PRODUCT_HEADER_FIELDS.WELL_COORDINATES.name(), well.getWell().getCoordinatesString());
row.put(PathwayProductAnalysis.PATHWAY_PRODUCT_HEADER_FIELDS.MSID.name(), well.getWell().getMsid());
row.put(PathwayProductAnalysis.PATHWAY_PRODUCT_HEADER_FIELDS.CONSTRUCT_ID.name(), well.getWell().getComposition());
row.put(PathwayProductAnalysis.PATHWAY_PRODUCT_HEADER_FIELDS.METLIN_ION.name(), pathwayStepIon);
writer.append(row);
writer.flush();
}
}
}