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