com.act.lcms.db.analysis.PathwayProductAnalysis.java Source code

Java tutorial

Introduction

Here is the source code for com.act.lcms.db.analysis.PathwayProductAnalysis.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.db.analysis;

import com.act.lcms.db.model.CuratedStandardMetlinIon;
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.io.DB;
import com.act.lcms.db.io.LoadPlateCompositionIntoDB;
import com.act.lcms.db.model.ChemicalAssociatedWithPathway;
import com.act.lcms.db.model.LCMSWell;
import com.act.lcms.db.model.MS1ScanForWellAndMassCharge;
import com.act.lcms.db.model.Plate;
import com.act.lcms.db.model.ScanFile;
import com.act.lcms.db.model.StandardIonResult;
import com.act.lcms.db.model.StandardWell;
import com.act.lcms.plotter.WriteAndPlotMS1Results;
import org.apache.commons.cli.CommandLine;
import org.apache.commons.cli.CommandLineParser;
import org.apache.commons.cli.DefaultParser;
import org.apache.commons.cli.HelpFormatter;
import org.apache.commons.cli.Option;
import org.apache.commons.cli.Options;
import org.apache.commons.cli.ParseException;
import org.apache.commons.lang3.StringUtils;
import org.apache.commons.lang3.tuple.Pair;

import java.io.File;
import java.io.FileOutputStream;
import java.io.IOException;
import java.sql.SQLException;
import java.util.ArrayList;
import java.util.Collections;
import java.util.Comparator;
import java.util.HashMap;
import java.util.HashSet;
import java.util.LinkedHashMap;
import java.util.List;
import java.util.Map;
import java.util.Set;

public class PathwayProductAnalysis {
    public static final String DEFAULT_SEARCH_ION = "M+H";

    public static final String OPTION_DIRECTORY = "d";
    public static final String OPTION_STRAINS = "s";
    public static final String OPTION_CONSTRUCT = "c";
    public static final String OPTION_NEGATIVE_STRAINS = "S";
    public static final String OPTION_NEGATIVE_CONSTRUCTS = "C";
    public static final String OPTION_OUTPUT_PREFIX = "o";
    public static final String OPTION_STANDARD_PLATE_BARCODE = "sp";
    public static final String OPTION_STANDARD_WELLS = "sw";
    public static final String OPTION_FILTER_BY_PLATE_BARCODE = "p";
    public static final String OPTION_USE_HEATMAP = "e";
    public static final String OPTION_SEARCH_ION = "i";
    public static final String OPTION_PATHWAY_SEARCH_IONS = "I";
    public static final String OPTION_ALLOW_MISSING_STANDARDS = "M";
    public static final String OPTION_USE_SNR = "r";
    public static final String OPTION_PLOTTING_DIR = "D";

    public static final String HELP_MESSAGE = StringUtils.join(new String[] {
            "This class applies the MS1 LCMS analysis to a combination of ",
            "standards and samples for the specified ion of all intermediate, side-reaction, ",
            "and final products of a given construct (optionally filtering samples by strains).  ",
            "An appropriate standard and any specified negative controls will be plotted alongside ",
            "a sample analysis for each product.",
            "Example cl options:  -d MNT_DATA_LEVEL1/lcms-ms1 -c ca1 -C ta1,on1,cr1 -o ca1-p7447-pathway -sp 7291 "
                    + "--intermediate-ions M+H,M+H,M+H,M+H -p 7447 --heat-map", },
            "");
    public static final HelpFormatter HELP_FORMATTER = new HelpFormatter();
    static {
        HELP_FORMATTER.setWidth(100);
    }

    public static final List<Option.Builder> OPTION_BUILDERS = new ArrayList<Option.Builder>() {
        {
            add(Option.builder(OPTION_DIRECTORY).argName("directory")
                    .desc("The directory where LCMS analysis results live").hasArg().required()
                    .longOpt("data-dir"));
            add(Option.builder(OPTION_OUTPUT_PREFIX).argName("output prefix")
                    .desc("A prefix for the output data/pdf files").hasArg().required().longOpt("output-prefix"));
            add(Option.builder(OPTION_CONSTRUCT).argName("construct id").desc(
                    "A construct whose intermediate/side-reaction products should be searched for in the traces")
                    .hasArg().longOpt("construct-id"));
            add(Option.builder(OPTION_STRAINS).argName("strains")
                    .desc("Filter analyzed LCMS samples to only these strains").hasArgs().valueSeparator(',')
                    .longOpt("msids"));
            add(Option.builder(OPTION_NEGATIVE_STRAINS).argName("negative-strains")
                    .desc("A strains to use as a negative control (the first novel LCMS sample will be used)")
                    .hasArgs().valueSeparator(',').longOpt("negative-msids"));
            add(Option.builder(OPTION_NEGATIVE_CONSTRUCTS).argName("constructs")
                    .desc("A constructs to use as a negative control (the first novel LCMS sample will be used)")
                    .hasArgs().valueSeparator(',').longOpt("negative-construct-ids"));
            add(Option.builder(OPTION_STANDARD_PLATE_BARCODE).argName("standard plate barcode")
                    .desc("The plate barcode to use when searching for a compatible standard").hasArg()
                    .longOpt("standard-plate"));
            add(Option.builder(OPTION_STANDARD_WELLS).argName("standard wells")
                    .desc("A list of well coordinates for standards, either offset in the pathway or a mapping of "
                            + "intermediate product to well (like paracetamol=A1,chorismate=C2)")
                    .hasArgs().valueSeparator(',').longOpt("standard-wells"));
            add(Option.builder(OPTION_SEARCH_ION).argName("search ion")
                    .desc("The ion for which to search (default is " + DEFAULT_SEARCH_ION + "); if used with - "
                            + OPTION_PATHWAY_SEARCH_IONS + " this will be the default for unspecified steps")
                    .hasArg().longOpt("search-ion"));
            add(Option.builder(OPTION_PATHWAY_SEARCH_IONS).desc(
                    "A list of ions per step, either by offset in the pathway (ultimate target first), or a mapping of "
                            + "intermediate product to ion (like paracetamol=M+H,chorismate=M+K)")
                    .hasArgs().valueSeparator(',').longOpt("intermediate-ions"));
            add(Option.builder(OPTION_FILTER_BY_PLATE_BARCODE).argName("plate barcode list")
                    .desc("A list of plate barcodes to consider, all other plates will be ignored").hasArgs()
                    .valueSeparator(',').longOpt("include-plates"));
            add(Option.builder(OPTION_USE_HEATMAP).desc("Produce a heat map rather than a 2d line plot")
                    .longOpt("heat-map"));
            add(Option.builder(OPTION_ALLOW_MISSING_STANDARDS)
                    .desc("Don't error when the standard for a pathway step can't be found")
                    .longOpt("allow-missing-standards"));
            add(Option.builder(OPTION_USE_SNR)
                    .desc("Use signal-to-noise ratio instead of max intensity for peak identification")
                    .longOpt("use-snr"));

            add(Option.builder().argName("font scale").desc(
                    "A Gnuplot fontscale value, should be between 0.1 and 0.5 (0.4 works if the graph text is large")
                    .hasArg().longOpt("font-scale"));
            add(Option.builder()
                    .desc(String.format(
                            "Use fine-grained M/Z tolerance (%.3f) when conducting the MS1 analysis "
                                    + "instead of default M/Z tolerance %.3f",
                            MS1.MS1_MZ_TOLERANCE_FINE, MS1.MS1_MZ_TOLERANCE_DEFAULT))
                    .longOpt("fine-grained-mz"));

            // Everybody needs a little help from their friends.
            add(Option.builder("h").argName("help").desc("Prints this help message").longOpt("help"));

            add(Option.builder(OPTION_PLOTTING_DIR).argName("plotting directory")
                    .desc("The absolute path of the plotting directory").hasArg().required()
                    .longOpt("plotting-dir"));
        }
    };
    static {
        // Add DB connection options.
        OPTION_BUILDERS.addAll(DB.DB_OPTION_BUILDERS);
    }

    public enum PATHWAY_PRODUCT_HEADER_FIELDS {
        TARGET_CHEMICAL, TYPE, FED_CHEMICAL, DETECTED, INTENSITY, TIME, PLATE_BARCODE, MODE, WELL_COORDINATES, MSID, CONSTRUCT_ID, METLIN_ION
    };

    private static final Set<String> EMPTY_SET = Collections.unmodifiableSet(new HashSet<>(0));

    public static void main(String[] args) throws Exception {
        Options opts = new Options();
        for (Option.Builder b : OPTION_BUILDERS) {
            opts.addOption(b.build());
        }

        CommandLine cl = null;
        try {
            CommandLineParser parser = new DefaultParser();
            cl = parser.parse(opts, args);
        } catch (ParseException e) {
            System.err.format("Argument parsing failed: %s\n", e.getMessage());
            HELP_FORMATTER.printHelp(LoadPlateCompositionIntoDB.class.getCanonicalName(), HELP_MESSAGE, opts, null,
                    true);
            System.exit(1);
        }

        if (cl.hasOption("help")) {
            HELP_FORMATTER.printHelp(LoadPlateCompositionIntoDB.class.getCanonicalName(), HELP_MESSAGE, opts, null,
                    true);
            return;
        }

        File lcmsDir = new File(cl.getOptionValue(OPTION_DIRECTORY));
        if (!lcmsDir.isDirectory()) {
            System.err.format("File at %s is not a directory\n", lcmsDir.getAbsolutePath());
            HELP_FORMATTER.printHelp(LoadPlateCompositionIntoDB.class.getCanonicalName(), HELP_MESSAGE, opts, null,
                    true);
            System.exit(1);
        }

        Double fontScale = null;
        if (cl.hasOption("font-scale")) {
            try {
                fontScale = Double.parseDouble(cl.getOptionValue("font-scale"));
            } catch (IllegalArgumentException e) {
                System.err.format("Argument for font-scale must be a floating point number.\n");
                System.exit(1);
            }
        }

        try (DB db = DB.openDBFromCLI(cl)) {
            Set<Integer> takeSamplesFromPlateIds = null;
            if (cl.hasOption(OPTION_FILTER_BY_PLATE_BARCODE)) {
                String[] plateBarcodes = cl.getOptionValues(OPTION_FILTER_BY_PLATE_BARCODE);
                System.out.format("Considering only sample wells in plates: %s\n",
                        StringUtils.join(plateBarcodes, ", "));
                takeSamplesFromPlateIds = new HashSet<>(plateBarcodes.length);
                for (String plateBarcode : plateBarcodes) {
                    Plate p = Plate.getPlateByBarcode(db, plateBarcode);
                    if (p == null) {
                        System.err.format("WARNING: unable to find plate in DB with barcode %s\n", plateBarcode);
                    } else {
                        takeSamplesFromPlateIds.add(p.getId());
                    }
                }
                // Allow filtering on barcode even if we couldn't find any in the DB.
            }

            System.out.format("Loading/updating LCMS scan files into DB\n");
            ScanFile.insertOrUpdateScanFilesInDirectory(db, lcmsDir);

            System.out.format("Processing LCMS scans\n");
            Pair<List<LCMSWell>, Set<Integer>> positiveWellsAndPlateIds = Utils.extractWellsAndPlateIds(db,
                    cl.getOptionValues(OPTION_STRAINS), cl.getOptionValues(OPTION_CONSTRUCT),
                    takeSamplesFromPlateIds, false);
            List<LCMSWell> positiveWells = positiveWellsAndPlateIds.getLeft();
            if (positiveWells.size() == 0) {
                throw new RuntimeException("Found no LCMS wells for specified strains/constructs");
            }
            // Only take negative samples from the plates where we found the positive samples.
            Pair<List<LCMSWell>, Set<Integer>> negativeWellsAndPlateIds = Utils.extractWellsAndPlateIds(db,
                    cl.getOptionValues(OPTION_NEGATIVE_STRAINS), cl.getOptionValues(OPTION_NEGATIVE_CONSTRUCTS),
                    positiveWellsAndPlateIds.getRight(), true);
            List<LCMSWell> negativeWells = negativeWellsAndPlateIds.getLeft();
            if (negativeWells == null || negativeWells.size() == 0) {
                System.err.format("WARNING: no valid negative samples found in same plates as positive samples\n");
            }

            // Extract the chemicals in the pathway and their product masses, then look up info on those chemicals
            List<Pair<ChemicalAssociatedWithPathway, Double>> productMasses = Utils
                    .extractMassesForChemicalsAssociatedWithConstruct(db, cl.getOptionValue(OPTION_CONSTRUCT));
            List<Pair<String, Double>> searchMZs = new ArrayList<>(productMasses.size());
            List<ChemicalAssociatedWithPathway> pathwayChems = new ArrayList<>(productMasses.size());
            for (Pair<ChemicalAssociatedWithPathway, Double> productMass : productMasses) {
                String chemName = productMass.getLeft().getChemical();
                searchMZs.add(Pair.of(chemName, productMass.getRight()));
                pathwayChems.add(productMass.getLeft());
            }
            System.out.format("Searching for intermediate/side-reaction products:\n");
            for (Pair<String, Double> searchMZ : searchMZs) {
                System.out.format("  %s: %.3f\n", searchMZ.getLeft(), searchMZ.getRight());
            }

            // Look up the standard by name.
            List<StandardWell> standardWells = new ArrayList<>();
            if (cl.hasOption(OPTION_STANDARD_WELLS)) {
                Plate standardPlate = Plate.getPlateByBarcode(db, cl.getOptionValue(OPTION_STANDARD_PLATE_BARCODE));
                Map<Integer, StandardWell> pathwayIdToStandardWell = extractStandardWellsFromOptionsList(db,
                        pathwayChems, cl.getOptionValues(OPTION_STANDARD_WELLS), standardPlate);
                for (ChemicalAssociatedWithPathway c : pathwayChems) { // TODO: we can avoid this loop.
                    StandardWell well = pathwayIdToStandardWell.get(c.getId());
                    if (well != null) {
                        standardWells.add(well);
                    }
                }
            } else {
                for (ChemicalAssociatedWithPathway c : pathwayChems) {
                    String standardName = c.getChemical();
                    System.out.format("Searching for well containing standard %s\n", standardName);
                    List<StandardWell> wells = StandardIonAnalysis.getStandardWellsForChemical(db, c.getChemical());
                    if (wells != null) {
                        standardWells.addAll(wells);
                    }
                }
            }

            boolean useFineGrainedMZ = cl.hasOption("fine-grained-mz");
            boolean useSNR = cl.hasOption(OPTION_USE_SNR);

            /* Process the standard, positive, and negative wells, producing ScanData containers that will allow them to be
             * iterated over for graph writing. We do not need to specify granular includeIons and excludeIons since
             * this would not take advantage of our caching strategy which uses a list of metlin ions as an index. */
            HashMap<Integer, Plate> plateCache = new HashMap<>();
            Pair<List<ScanData<StandardWell>>, Double> allStandardScans = AnalysisHelper.processScans(db, lcmsDir,
                    searchMZs, ScanData.KIND.STANDARD, plateCache, standardWells, useFineGrainedMZ, EMPTY_SET,
                    EMPTY_SET, useSNR);
            Pair<List<ScanData<LCMSWell>>, Double> allPositiveScans = AnalysisHelper.processScans(db, lcmsDir,
                    searchMZs, ScanData.KIND.POS_SAMPLE, plateCache, positiveWells, useFineGrainedMZ, EMPTY_SET,
                    EMPTY_SET, useSNR);
            Pair<List<ScanData<LCMSWell>>, Double> allNegativeScans = AnalysisHelper.processScans(db, lcmsDir,
                    searchMZs, ScanData.KIND.NEG_CONTROL, plateCache, negativeWells, useFineGrainedMZ, EMPTY_SET,
                    EMPTY_SET, useSNR);

            String fmt = "pdf";
            String outImg = cl.getOptionValue(OPTION_OUTPUT_PREFIX) + "." + fmt;
            String outData = cl.getOptionValue(OPTION_OUTPUT_PREFIX) + ".data";
            String outAnalysis = cl.getOptionValue(OPTION_OUTPUT_PREFIX) + ".tsv";

            System.err.format("Writing combined scan data to %s and graphs to %s\n", outData, outImg);
            String plottingDirectory = cl.getOptionValue(OPTION_PLOTTING_DIR);

            List<ScanData<LCMSWell>> posNegWells = new ArrayList<>();
            posNegWells.addAll(allPositiveScans.getLeft());
            posNegWells.addAll(allNegativeScans.getLeft());

            Map<Integer, String> searchIons;
            if (cl.hasOption(OPTION_PATHWAY_SEARCH_IONS)) {
                searchIons = extractPathwayStepIons(pathwayChems, cl.getOptionValues(OPTION_PATHWAY_SEARCH_IONS),
                        cl.getOptionValue(OPTION_SEARCH_ION, "M+H"));
                /* This is pretty lazy, but works with the existing API.  Extract all selected ions for all search masses when
                 * performing the scan, then filter down to the desired ions for the plot at the end.
                 * TODO: specify the masses and scans per sample rather than batching everything together.  It might be slower,
                 * but it'll be clearer to read. */
            } else {
                // We need to make sure that the standard metlin ion we choose is consistent with the ion modes of
                // the given positive, negative and standard scan files. For example, we should not pick a negative
                // metlin ion if all our available positive control scan files are in the positive ion mode.
                Map<Integer, Pair<Boolean, Boolean>> ionModes = new HashMap<>();
                for (ChemicalAssociatedWithPathway chemical : pathwayChems) {
                    boolean isPositiveScanPresent = false;
                    boolean isNegativeScanPresent = false;

                    for (ScanData<StandardWell> scan : allStandardScans.getLeft()) {
                        if (chemical.getChemical().equals(scan.getWell().getChemical())
                                && chemical.getChemical().equals(scan.getTargetChemicalName())) {
                            if (MS1.IonMode.valueOf(
                                    scan.getScanFile().getMode().toString().toUpperCase()) == MS1.IonMode.POS) {
                                isPositiveScanPresent = true;
                            }

                            if (MS1.IonMode.valueOf(
                                    scan.getScanFile().getMode().toString().toUpperCase()) == MS1.IonMode.NEG) {
                                isNegativeScanPresent = true;
                            }
                        }
                    }

                    for (ScanData<LCMSWell> scan : posNegWells) {
                        if (chemical.getChemical().equals(scan.getWell().getChemical())
                                && chemical.getChemical().equals(scan.getTargetChemicalName())) {
                            if (MS1.IonMode.valueOf(
                                    scan.getScanFile().getMode().toString().toUpperCase()) == MS1.IonMode.POS) {
                                isPositiveScanPresent = true;
                            }

                            if (MS1.IonMode.valueOf(
                                    scan.getScanFile().getMode().toString().toUpperCase()) == MS1.IonMode.NEG) {
                                isNegativeScanPresent = true;
                            }
                        }
                    }

                    ionModes.put(chemical.getId(), Pair.of(isPositiveScanPresent, isNegativeScanPresent));
                }

                // Sort in descending order of media where MeOH and Water related media are promoted to the top and
                // anything derived from yeast media are demoted. We do this because we want to first process the water
                // and meoh media before processing the yeast media since the yeast media depends on the analysis of the former.
                Collections.sort(standardWells, new Comparator<StandardWell>() {
                    @Override
                    public int compare(StandardWell o1, StandardWell o2) {
                        if (StandardWell.doesMediaContainYeastExtract(o1.getMedia())
                                && !StandardWell.doesMediaContainYeastExtract(o2.getMedia())) {
                            return 1;
                        } else {
                            return 0;
                        }
                    }
                });

                searchIons = extractPathwayStepIonsFromStandardIonAnalysis(pathwayChems, lcmsDir, db, standardWells,
                        plottingDirectory, ionModes);
            }

            produceLCMSPathwayHeatmaps(lcmsDir, outData, outImg, outAnalysis, pathwayChems, allStandardScans,
                    allPositiveScans, allNegativeScans, fontScale, cl.hasOption(OPTION_USE_HEATMAP), searchIons);
        }
    }

    private static Map<Integer, String> extractPathwayStepIonsFromStandardIonAnalysis(
            List<ChemicalAssociatedWithPathway> pathwayChems, File lcmsDir, DB db, List<StandardWell> standardWells,
            String plottingDir, Map<Integer, Pair<Boolean, Boolean>> ionModesAvailable) throws Exception {

        Map<Integer, String> result = new HashMap<>();

        for (ChemicalAssociatedWithPathway pathwayChem : pathwayChems) {

            Map<StandardWell, StandardIonResult> wellToIonRanking = StandardIonAnalysis
                    .getBestMetlinIonsForChemical(pathwayChem.getChemical(), lcmsDir, db, standardWells,
                            plottingDir);

            Pair<Boolean, Boolean> modes = ionModesAvailable.get(pathwayChem.getId());

            Map<StandardIonResult, String> chemicalToCuratedMetlinIon = new HashMap<>();
            List<StandardIonResult> standardIonResults = new ArrayList<>(wellToIonRanking.values());

            for (StandardIonResult standardIonResult : standardIonResults) {
                Integer manualOverrideId = standardIonResult.getManualOverrideId();
                if (manualOverrideId != null) {
                    chemicalToCuratedMetlinIon.put(standardIonResult,
                            CuratedStandardMetlinIon.getBestMetlinIon(db, manualOverrideId).getBestMetlinIon());
                }
            }

            String bestMetlinIon = AnalysisHelper.scoreAndReturnBestMetlinIonFromStandardIonResults(
                    standardIonResults, chemicalToCuratedMetlinIon, modes.getLeft(), modes.getRight());

            if (bestMetlinIon != null) {
                result.put(pathwayChem.getId(), bestMetlinIon);
            } else {
                result.put(pathwayChem.getId(), DEFAULT_SEARCH_ION);
            }
        }
        return result;
    }

    private static Map<Integer, StandardWell> extractStandardWellsFromOptionsList(DB db,
            List<ChemicalAssociatedWithPathway> pathwayChems, String[] optionValues, Plate standardPlate)
            throws SQLException, IOException, ClassNotFoundException {
        Map<String, String> chemToWellByName = new HashMap<>();
        Map<Integer, String> chemToWellByIndex = new HashMap<>();
        if (optionValues != null && optionValues.length > 0) {
            for (int i = 0; i < optionValues.length; i++) {
                String[] fields = StringUtils.split(optionValues[i], "=");
                if (fields != null && fields.length == 2) {
                    chemToWellByName.put(fields[0], fields[1]);
                } else {
                    chemToWellByIndex.put(i, optionValues[i]);
                }
            }
        }

        Map<Integer, StandardWell> results = new HashMap<>();
        for (int i = 0; i < pathwayChems.size(); i++) {
            ChemicalAssociatedWithPathway chem = pathwayChems.get(i);
            String coords = null;
            if (chemToWellByName.containsKey(chem.getChemical())) {
                coords = chemToWellByName.remove(chem.getChemical());
            } else if (chemToWellByIndex.containsKey(i)) {
                coords = chemToWellByIndex.remove(i);
            }

            if (coords == null) {
                System.err.format("No coordinates specified for %s, skipping\n", chem.getChemical());
                continue;
            }
            Pair<Integer, Integer> intCoords = Utils.parsePlateCoordinates(coords);
            StandardWell well = StandardWell.getInstance().getStandardWellsByPlateIdAndCoordinates(db,
                    standardPlate.getId(), intCoords.getLeft(), intCoords.getRight());
            if (well == null) {
                System.err.format("ERROR: Could not find well %s in plate %s\n", coords,
                        standardPlate.getBarcode());
                System.exit(-1);
            } else if (!well.getChemical().equals(chem.getChemical())) {
                System.err.format(
                        "WARNING: pathway chemical %s and chemical in specified standard well %s don't match!\n",
                        chem.getChemical(), well.getChemical());
            }

            System.out.format("Using standard well %s : %s for pathway chemical %s (step %d)\n",
                    standardPlate.getBarcode(), coords, chem.getChemical(), chem.getIndex());

            results.put(chem.getId(), well);
        }

        return results;
    }

    private static Map<Integer, String> extractPathwayStepIons(List<ChemicalAssociatedWithPathway> pathwayChems,
            String[] optionValues, String defaultIon) {
        Map<String, String> pathwayToIonByName = new HashMap<>();
        Map<Integer, String> pathwayToIonByIndex = new HashMap<>();
        if (optionValues != null && optionValues.length > 0) {
            for (int i = 0; i < optionValues.length; i++) {
                String[] fields = StringUtils.split(optionValues[i], "=");
                if (fields != null && fields.length == 2) {
                    if (!MS1.VALID_MS1_IONS.contains(fields[1])) {
                        System.err.format("WARNING: found invalid intermediate/ion pair, skipping: %s\n",
                                optionValues[i]);
                        continue;
                    }
                    pathwayToIonByName.put(fields[0], fields[1]);
                } else {
                    pathwayToIonByIndex.put(i, optionValues[i]);
                }
            }
        }

        Map<Integer, String> results = new HashMap<>();
        for (int i = 0; i < pathwayChems.size(); i++) {
            ChemicalAssociatedWithPathway chem = pathwayChems.get(i);
            String ion = defaultIon;
            if (pathwayToIonByName.containsKey(chem.getChemical())) {
                ion = pathwayToIonByName.remove(chem.getChemical());
            } else if (pathwayToIonByIndex.containsKey(i)) {
                ion = pathwayToIonByIndex.remove(i);
            }
            System.out.format("Using ion %s for pathway chemical %s (step %d)\n", ion, chem.getChemical(),
                    chem.getIndex());
            results.put(chem.getId(), ion);
        }

        if (!(pathwayToIonByName.isEmpty() && pathwayToIonByIndex.isEmpty())) {
            System.err.format("WARNING: unable to assign some pathway ions by name/index:\n");
            for (Map.Entry<String, String> entry : pathwayToIonByName.entrySet()) {
                System.err.format("  %s => %s\n", entry.getKey(), entry.getValue());
            }
            for (Map.Entry<Integer, String> entry : pathwayToIonByIndex.entrySet()) {
                System.err.format("  %d => %s\n", entry.getKey(), entry.getValue());
            }
        }

        return results;
    }

    private static final Comparator<ScanData<LCMSWell>> LCMS_SCAN_COMPARATOR = new Comparator<ScanData<LCMSWell>>() {
        @Override
        public int compare(ScanData<LCMSWell> o1, ScanData<LCMSWell> o2) {
            int c;
            // TODO: consider feeding conditions in sort to match condition order to steps.
            c = o1.getWell().getMsid().compareTo(o2.getWell().getMsid());
            if (c != 0)
                return c;
            c = o1.getPlate().getBarcode().compareTo(o2.getPlate().getBarcode());
            if (c != 0)
                return c;
            c = o1.getWell().getPlateRow().compareTo(o2.getWell().getPlateRow());
            if (c != 0)
                return c;
            c = o1.getWell().getPlateColumn().compareTo(o2.getWell().getPlateColumn());
            if (c != 0)
                return c;
            c = o1.getScanFile().getFilename().compareTo(o2.getScanFile().getFilename());
            return c;
        }
    };

    private static final ScanData<LCMSWell> BLANK_SCAN = new ScanData<>(ScanData.KIND.BLANK, null, null, null, null,
            null, null);

    public static void produceLCMSPathwayHeatmaps(File lcmsDir, String outData, String outImg, String outAnalysis,
            List<ChemicalAssociatedWithPathway> pathwayChems,
            Pair<List<ScanData<StandardWell>>, Double> allStandardScans,
            Pair<List<ScanData<LCMSWell>>, Double> allPositiveScans,
            Pair<List<ScanData<LCMSWell>>, Double> allNegativeScans, Double fontScale, boolean makeHeatmaps,
            Map<Integer, String> searchIons) throws Exception {
        String fmt = "pdf";
        System.err.format("Writing combined scan data to %s and graphs to %s\n", outData, outImg);

        // Generate the data file and graphs.
        try (FileOutputStream fos = new FileOutputStream(outData)) {
            List<String> graphLabels = new ArrayList<>();
            List<Double> yMaxList = new ArrayList<>();

            List<String> pathwayProductHeaderFields = new ArrayList<>();
            for (PATHWAY_PRODUCT_HEADER_FIELDS field : PATHWAY_PRODUCT_HEADER_FIELDS.values()) {
                pathwayProductHeaderFields.add(field.name());
            }

            TSVWriter<String, String> resultsWriter = new TSVWriter<>(pathwayProductHeaderFields);
            resultsWriter.open(new File(outAnalysis));

            for (ChemicalAssociatedWithPathway chem : pathwayChems) {
                System.out.format("Processing data for pathway chemical %s\n", chem.getChemical());

                Double maxIntensity = 0.0d;

                String pathwayStepIon = null;
                if (searchIons != null && searchIons.containsKey(chem.getId())) {
                    pathwayStepIon = searchIons.get(chem.getId());
                }

                MS1.IonMode mode = MS1.getIonModeOfIon(pathwayStepIon);

                // Extract the first available
                List<ScanData<StandardWell>> stdScan = new ArrayList<>();
                for (ScanData<StandardWell> scan : allStandardScans.getLeft()) {
                    if (chem.getChemical().equals(scan.getWell().getChemical())
                            && chem.getChemical().equals(scan.getTargetChemicalName())) {
                        if (mode.toString().toLowerCase()
                                .equals(scan.getScanFile().getMode().toString().toLowerCase())) {
                            stdScan.add(scan);
                            MS1ScanForWellAndMassCharge scanResults = scan.getMs1ScanResults();
                            Double intensity = pathwayStepIon == null ? scanResults.getMaxYAxis()
                                    : scanResults.getMaxIntensityForIon(pathwayStepIon);
                            if (intensity != null) {
                                maxIntensity = Math.max(maxIntensity, intensity);
                            }
                        }
                    }
                }
                if (stdScan.size() == 0) {
                    System.err.format("WARNING: unable to find standard well scan for chemical %s\n",
                            chem.getChemical());
                }

                List<ScanData<LCMSWell>> matchingPosScans = new ArrayList<>();
                for (ScanData<LCMSWell> scan : allPositiveScans.getLeft()) {
                    if (chem.getChemical().equals(scan.getTargetChemicalName())) {
                        if (mode.toString().toLowerCase()
                                .equals(scan.getScanFile().getMode().toString().toLowerCase())) {
                            matchingPosScans.add(scan);
                            MS1ScanForWellAndMassCharge scanResults = scan.getMs1ScanResults();
                            Double intensity = pathwayStepIon == null ? scanResults.getMaxYAxis()
                                    : scanResults.getMaxIntensityForIon(pathwayStepIon);
                            if (intensity != null) {
                                maxIntensity = Math.max(maxIntensity, intensity);
                            }
                        }
                    }
                }
                matchingPosScans.sort(LCMS_SCAN_COMPARATOR);

                List<ScanData<LCMSWell>> matchingNegScans = new ArrayList<>();
                for (ScanData<LCMSWell> scan : allNegativeScans.getLeft()) {
                    if (chem.getChemical().equals(scan.getTargetChemicalName())) {
                        if (mode.toString().toLowerCase()
                                .equals(scan.getScanFile().getMode().toString().toLowerCase())) {
                            matchingNegScans.add(scan);
                            MS1ScanForWellAndMassCharge scanResults = scan.getMs1ScanResults();
                            Double intensity = pathwayStepIon == null ? scanResults.getMaxYAxis()
                                    : scanResults.getMaxIntensityForIon(pathwayStepIon);
                            if (intensity != null) {
                                maxIntensity = Math.max(maxIntensity, intensity);
                            }
                        }
                    }
                }
                matchingNegScans.sort(LCMS_SCAN_COMPARATOR);

                List<ScanData> allScanData = new ArrayList<>();
                List<ScanData<LCMSWell>> positiveAndNegativeData = new ArrayList<>();

                positiveAndNegativeData.addAll(matchingPosScans);
                positiveAndNegativeData.addAll(matchingNegScans);

                allScanData.addAll(stdScan);
                allScanData.addAll(positiveAndNegativeData);
                allScanData.add(BLANK_SCAN);

                Map<ScanData<LCMSWell>, XZ> result = WaveformAnalysis
                        .pickBestRepresentativeRetentionTimeFromStandardWells(stdScan, pathwayStepIon,
                                positiveAndNegativeData);

                WriteAndPlotMS1Results.writePathwayProductOutput(resultsWriter, chem.getChemical(),
                        positiveAndNegativeData, pathwayStepIon, result);

                Set<String> pathwayStepIons = pathwayStepIon == null ? null : Collections.singleton(pathwayStepIon);
                // Write all the scan data out to a single data file.
                for (ScanData scanData : allScanData) {
                    graphLabels.addAll(AnalysisHelper.writeScanData(fos, lcmsDir, maxIntensity, scanData,
                            makeHeatmaps, false, pathwayStepIons));
                }

                // Save one max intensity per graph so we can plot with them later.
                for (ScanData<LCMSWell> scan : allScanData) {
                    if (!scan.getKind().equals(ScanData.KIND.BLANK)) {
                        Double intensity = pathwayStepIon == null ? scan.getMs1ScanResults().getMaxYAxis()
                                : scan.getMs1ScanResults().getMaxIntensityForIon(pathwayStepIon);
                        yMaxList.add(intensity);
                    } else {
                        // Add a 0 intensity for the blank scan
                        yMaxList.add(0.0d);
                    }
                }
            }

            resultsWriter.close();

            // We need to pass the yMax values as an array to the Gnuplotter.
            Double[] yMaxes = yMaxList.toArray(new Double[yMaxList.size()]);
            Gnuplotter plotter = fontScale == null ? new Gnuplotter() : new Gnuplotter(fontScale);
            if (makeHeatmaps) {
                plotter.plotHeatmap(outData, outImg, graphLabels.toArray(new String[graphLabels.size()]), null, fmt,
                        11.0, 8.5, yMaxes, outImg + ".gnuplot");
            } else {
                plotter.plot2D(outData, outImg, graphLabels.toArray(new String[graphLabels.size()]), "time", null,
                        "intensity", fmt, null, null, yMaxes, outImg + ".gnuplot");
            }
        }
    }
}