com.act.lcms.plotter.WriteAndPlotMS1Results.java Source code

Java tutorial

Introduction

Here is the source code for com.act.lcms.plotter.WriteAndPlotMS1Results.java

Source

/*************************************************************************
*                                                                        *
*  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();
        }
    }
}