Java tutorial
/************************************************************************* * * * 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.biointerpretation.Utils; import chemaxon.marvin.io.MolExportException; import chemaxon.reaction.ConcurrentReactorProcessor; import chemaxon.reaction.ReactionException; import chemaxon.reaction.Reactor; import chemaxon.sss.SearchConstants; import chemaxon.sss.search.MolSearch; import chemaxon.sss.search.MolSearchOptions; import chemaxon.sss.search.SearchException; import chemaxon.struc.Molecule; import chemaxon.struc.RxnMolecule; import chemaxon.util.iterator.MoleculeIterator; import chemaxon.util.iterator.MoleculeIteratorFactory; import com.act.analysis.chemicals.molecules.MoleculeExporter; import com.act.analysis.chemicals.molecules.MoleculeFormat; import com.act.analysis.chemicals.molecules.MoleculeFormat$; import org.apache.commons.collections4.Bag; import org.apache.commons.collections4.bag.HashBag; import org.apache.commons.collections4.iterators.PermutationIterator; import org.apache.logging.log4j.LogManager; import org.apache.logging.log4j.Logger; import java.io.IOException; import java.util.ArrayList; import java.util.Arrays; import java.util.HashMap; import java.util.HashSet; import java.util.Iterator; import java.util.List; import java.util.Map; import java.util.Optional; import java.util.Set; import java.util.stream.Collectors; public class ReactionProjector { private static final Logger LOGGER = LogManager.getFormatterLogger(ReactionProjector.class); private static final String MOL_NOT_FOUND = "NOT_FOUND"; private static final MolSearchOptions LAX_SEARCH_OPTIONS = new MolSearchOptions(SearchConstants.SUBSTRUCTURE); private static final MolSearch DEFAULT_SEARCHER = new MolSearch(); private static final String DEFAULT_MOLECULE_FORMAT = MoleculeFormat.noAuxInchi$.MODULE$.toString(); private final String moleculeFormat; /** * Set searcher to ignore stereo and bond type. This will allow us to optimistically match products so that we don't * end up throwing out a reactor that produces the right compound in a slightly different form. */ static { LAX_SEARCH_OPTIONS.setStereoSearchType(SearchConstants.STEREO_IGNORE); LAX_SEARCH_OPTIONS.setVagueBondLevel(SearchConstants.VAGUE_BOND_LEVEL4); DEFAULT_SEARCHER.setSearchOptions(LAX_SEARCH_OPTIONS); } private transient ThreadLocal<MolSearch> substrateSearcher = new ThreadLocal<MolSearch>() { @Override protected MolSearch initialValue() { MolSearch search = new MolSearch(); /* These parameters were identified by trial and error as reasonable defaults for doing substructure searches for * reactor substrate matching. */ MolSearchOptions options = new MolSearchOptions(SearchConstants.SUBSTRUCTURE); // This allows H's in RO strings to match implicit hydrogens in our target molecules. options.setImplicitHMatching(SearchConstants.IMPLICIT_H_MATCHING_ENABLED); /* This allows for vague bond matching in ring structures. From the Chemaxon Docs: * In the query all single ring bonds are replaced by "single or aromatic" and all double ring bonds are * replaced by "double or aromatic" prior to search. * (https://www.chemaxon.com/jchem/doc/dev/java/api/chemaxon/sss/SearchConstants.html) * * This should allow us to handle aromatized molecules gracefully without handling non-ring single and double * bonds ambiguously. */ options.setVagueBondLevel(SearchConstants.VAGUE_BOND_LEVEL2); // Few if any of our ROs concern stereo chemistry, so we can just ignore it. options.setStereoSearchType(SearchConstants.STEREO_IGNORE); /* Chemaxon's tautomer handling is weird, as sometimes it picks a non-representative tautomer as its default. * As such, we'll allow tautomer matches to avoid excluding viable candidates. */ options.setTautomerSearch(SearchConstants.TAUTOMER_SEARCH_ON); search.setSearchOptions(options); return search; } }; private boolean useSubstructureFiltering = false; private MolSearch searcher; public ReactionProjector() { this(DEFAULT_SEARCHER, DEFAULT_MOLECULE_FORMAT); } public ReactionProjector(MolSearch searcher) { this(searcher, DEFAULT_MOLECULE_FORMAT); } public ReactionProjector(String moleculeFormat) { this(DEFAULT_SEARCHER, moleculeFormat); } public ReactionProjector(MolSearch searcher, String moleculeFormat) { this.searcher = searcher; this.moleculeFormat = moleculeFormat; } public ReactionProjector(boolean useSubstructureFiltering) { // I don't anticipate either of these ever being used if this is the constructor employed. this.searcher = DEFAULT_SEARCHER; this.moleculeFormat = DEFAULT_MOLECULE_FORMAT; this.useSubstructureFiltering = useSubstructureFiltering; } /** * This method should be called if projecting on more chemicals than should be stored in memory * simultaneously. Until this is called, all chemicals acted on by the projector will be cached along * with their inchis. */ public void clearInchiCache() { MoleculeExporter.clearCache(); } /** * Run the given reactor until it produces the expected product. Reactor must produce one product at a time. * * @param reactor The reactor to run. * @param expectedProduct The product we expect to see. * @return The produced product; this is necessary because the reactor will produce the product with atom maps * corresponding to the substrate, whereas the expectedProduct Molecule will not have such atom maps. * @throws ReactionException If the expected product is not produced at all by the Reactor. * @throws IOException */ public Molecule reactUntilProducesProduct(Reactor reactor, Molecule expectedProduct) throws ReactionException { if (reactor.getProductCount() != 1) { throw new IllegalArgumentException("Reactor must produce exactly one product."); } Molecule[] products; while ((products = reactor.react()) != null) { if (testEquality(products[0], expectedProduct)) { return products[0]; } } throw new ReactionException("Expected product not among Reactor's predictions."); } /** * Tests equality of molecules based on two substructure queries. This is desirable because it makes the * equality comparison very configurable as compared to comparing inchis. * * @param A One molecule. * @param B Another molecule. * @return True if they are equivalent. */ private boolean testEquality(Molecule A, Molecule B) { return substructureTest(A, B) && substructureTest(B, A); } private boolean substructureTest(Molecule substructure, Molecule superstructure) { searcher.setQuery(substructure); searcher.setTarget(superstructure); try { return searcher.findFirst() != null; } catch (SearchException e) { // Log error but don't propagate upward. Have never seen this before. LOGGER.error("Error on ReactionProjector.substructureTest(), %s", e.getMessage()); return false; } } /** * Get the results of a reaction in list form, rather than as a map from substrates to products. * * @param mols The substrates. * @param reactor The reactor. * @param maxProducts The maximum number of products to compute before exiting. No limit is imposed for max values * of null or 0. * @return A list of product sets produced by this reaction. * @throws IOException * @throws ReactionException */ public List<Molecule[]> getAllProjectedProductSets(Molecule[] mols, Reactor reactor, Integer maxProducts) throws IOException, ReactionException { Map<Molecule[], List<Molecule[]>> map = getRoProjectionMap(mols, reactor, maxProducts); List<Molecule[]> allProductSets = new ArrayList<>(); for (Map.Entry<Molecule[], List<Molecule[]>> entry : map.entrySet()) { allProductSets.addAll(entry.getValue()); } return allProductSets; } public List<Molecule[]> getAllProjectedProductSets(Molecule[] mols, Reactor reactor) throws IOException, ReactionException { return getAllProjectedProductSets(mols, reactor, null); } /** * This function takes as input an array of molecules and a Reactor and outputs the products of the transformation. * The results are returned as a map from orderings of the substrates to the products produced by those orderings. * In most cases the map will have only one entry, but in some cases different orderings of substrates can lead to * different valid predictions. * <p> * Note that, for efficient two substrate expansion, the specialized method * fastProjectionOfTwoSubstrateRoOntoTwoMolecules should be used instead of this one. * * @param mols An array of molecules representing the chemical reactants. * @param reactor A Reactor representing the reaction to apply. * @param maxProducts The maximum number of products to compute before exiting. No limit is imposed for max values * of null or 0. * @return The substrate -> product map. */ public Map<Molecule[], List<Molecule[]>> getRoProjectionMap(Molecule[] mols, Reactor reactor, Integer maxProducts) throws ReactionException, IOException { Map<Molecule[], List<Molecule[]>> resultsMap = new HashMap<>(); if (mols.length != reactor.getReactantCount()) { LOGGER.debug("Tried to project RO with %d substrates on set of %d substrates", reactor.getReactantCount(), mols.length); return resultsMap; } // If there is only one reactant, we can do just a simple reaction computation. However, if we have multiple reactants, // we have to use the ConcurrentReactorProcessor API since it gives us the ability to combinatorially explore all // possible matching combinations of reactants on the substrates of the RO. if (mols.length == 1) { List<Molecule[]> productSets = getProductsFixedOrder(reactor, mols, maxProducts); if (!productSets.isEmpty()) { resultsMap.put(mols, productSets); } } else if (useSubstructureFiltering) { Optional<List<Set<Integer>>> viableSubstrateIndexes; try { viableSubstrateIndexes = matchCandidatesToSubstrateStructures(reactor, mols); } catch (SearchException e) { throw new ReactionException("Caught exception when performing pre-reaction substructure search", e); } if (viableSubstrateIndexes.isPresent()) { List<Integer> allIndexes = new ArrayList<Integer>(mols.length) { { for (int i = 0; i < mols.length; i++) { add(i); } } }; int permutationIndex = 0; PermutationIterator<Integer> iter = new PermutationIterator<>(allIndexes); while (iter.hasNext()) { permutationIndex++; List<Integer> permutation = iter.next(); if (permutationFitsSubstructureMatches(permutation, viableSubstrateIndexes.get())) { Molecule[] substrates = indexPermutationToMolecules(mols, permutation); reactor.setReactants(substrates); Molecule[] products; List<Molecule[]> results = new ArrayList<>(); while ((products = reactor.react()) != null) { results.add(products); if (maxProducts != null && maxProducts > 0 && results.size() >= maxProducts) { break; } } if (results.size() > 0) { resultsMap.put(substrates, results); } } } } // Otherwise just return the empty map. } else { // TODO: why not make one of these per ReactionProjector object? // TODO: replace this with Apache commons PermutationIterator for clean iteration over distinct permutations. ConcurrentReactorProcessor reactorProcessor = new ConcurrentReactorProcessor(); reactorProcessor.setReactor(reactor); // This step is needed by ConcurrentReactorProcessor for the combinatorial exploration. // The iterator takes in the same substrates in each iteration step to build a matrix of combinations. // For example, if we have A + B -> C, the iterator creates this array: [[A,B], [A,B]]. Based on this, // ConcurrentReactorProcessor will cross the elements in the first index with the second, creating the combinations // A+A, A+B, B+A, B+B and operates all of those on the RO, which is what is desired. MoleculeIterator[] iterator = new MoleculeIterator[mols.length]; for (int i = 0; i < mols.length; i++) { iterator[i] = MoleculeIteratorFactory.createMoleculeIterator(mols); } reactorProcessor.setReactantIterators(iterator, ConcurrentReactorProcessor.MODE_COMBINATORIAL); // Bag is a multi-set class that ships with Apache commons collection, and we already use many commons libs--easy! Bag<Molecule> originalReactantsSet = new HashBag<>(Arrays.asList(mols)); // This set keeps track of substrate combinations we've used, and avoids repeats. Repeats can occur // when several substrates are identical, and can be put in "different" but symmetric orderings. Set<String> substrateHashes = new HashSet<>(); List<Molecule[]> results = null; int reactantCombination = 0; while ((results = reactorProcessor.reactNext()) != null) { reactantCombination++; Molecule[] reactants = reactorProcessor.getReactants(); if (results.isEmpty()) { LOGGER.debug("No results found for reactants combination %d, skipping", reactantCombination); continue; } Bag<Molecule> thisReactantSet = new HashBag<>(Arrays.asList(reactants)); if (!originalReactantsSet.equals(thisReactantSet)) { LOGGER.debug("Reactant set %d does not represent original, complete reactant sets, skipping", reactantCombination); continue; } String thisHash = getStringHash(reactants); if (substrateHashes.contains(thisHash)) { continue; } substrateHashes.add(thisHash); resultsMap.put(reactants, results); if (maxProducts != null && maxProducts > 0 && resultsMap.size() >= maxProducts) { break; } } } return resultsMap; } private Optional<List<Set<Integer>>> matchCandidatesToSubstrateStructures(Reactor reactor, Molecule[] candidates) throws SearchException, MolExportException { RxnMolecule rxnMol = reactor.getReaction(); Molecule[] substrateStructures = rxnMol.getReactants(); MolSearch search = substrateSearcher.get(); /* Each list in `results` is a set of viable indexes into candidates for a given substrate position in the reactor. * So if candidates contains three molecules and results looks like: * 0: 0, 1 * 1: 1 * 2: 0, 2 * * Then candidates[0] can appear as substrate 0 or 2, candidates[1] as substrate 0 or 1, and candidates[2] only as * substrate 2. */ List<Set<Integer>> results = new ArrayList<>(substrateStructures.length); for (int i = 0; i < substrateStructures.length; i++) { search.setQuery(substrateStructures[i]); results.add(new HashSet<>()); /* Test each candidate molecule against each substrate structure of the reactor. If we can't find a match, * there's no way that the molecule could appear in a particular position in the reaction and we can later * eliminate input molecule permutations that don't conform to this structure match. */ for (int j = 0; j < candidates.length; j++) { search.setTarget(candidates[j]); // Any match will do, so just call findFirst and ensure we get some sort of results. int[] res = search.findFirst(); if (res != null) { results.get(i).add(j); } } /* If one of the substrate patterns can't be matched to any input structure, then there are no viable permutations * to consider. We tell the call as much by returning `empty`. */ if (results.get(i).size() == 0) { return Optional.empty(); } } return Optional.of(results); } private Boolean permutationFitsSubstructureMatches(List<Integer> permutation, List<Set<Integer>> viableIndexes) { /* For each permutation, match its value at each location against the list of viable indexes for that location. * If a given permutation is 0, 2, 1 and viableIndexes looks like * 0: 0, 1 * 1: 1 * 2: 0, 2 * as above, then then the second structure index doesn't pass a substructure search and we can skip it. */ for (int i = 0; i < permutation.size(); i++) { if (!viableIndexes.get(i).contains(permutation.get(i))) { return false; } } return true; } private Molecule[] indexPermutationToMolecules(Molecule[] molecules, List<Integer> permutation) { Molecule[] results = new Molecule[molecules.length]; for (int i = 0; i < permutation.size(); i++) { results[i] = molecules[permutation.get(i)]; } return results; } public Map<Molecule[], List<Molecule[]>> getRoProjectionMap(Molecule[] mols, Reactor reactor) throws ReactionException, IOException { return getRoProjectionMap(mols, reactor, null); } /** * This function projects all possible combinations of two input substrates onto a 2 substrate RO, then * cleans and returns the results of that projection. * TODO: Expand this class to handle 3 or 4 substrate reactions. * * @param mols Substrate molecules * @param reactor The two substrate reactor * @return A list of product arrays, where each array represents the products of a given reaction combination. * @throws ReactionException * @throws IOException */ public Map<Molecule[], List<Molecule[]>> fastProjectionOfTwoSubstrateRoOntoTwoMolecules(Molecule[] mols, Reactor reactor) throws ReactionException, IOException { Map<Molecule[], List<Molecule[]>> results = new HashMap<>(); Molecule[] firstCombinationOfSubstrates = new Molecule[] { mols[0], mols[1] }; List<Molecule[]> productSets = getProductsFixedOrder(reactor, firstCombinationOfSubstrates); if (!productSets.isEmpty()) { results.put(firstCombinationOfSubstrates, productSets); } // Second ordering is same if two molecules are equal. if (getMoleculeString(mols[0]).equals(getMoleculeString(mols[1]))) { return results; } Molecule[] secondCombinationOfSubstrates = new Molecule[] { mols[1], mols[0] }; productSets = getProductsFixedOrder(reactor, firstCombinationOfSubstrates); if (!productSets.isEmpty()) { results.put(secondCombinationOfSubstrates, productSets); } return results; } /** * Gets all product sets that a Reactor produces when applies to a given array of substrates, in only the order * that they are already in. * * @param reactor The Reactor. * @param substrates The substrates. * @param maxProducts The maximum number of products to compute before exiting. No limit is imposed for max values * of null or 0. * @return A list of product arrays returned by the Reactor. * @throws ReactionException */ public List<Molecule[]> getProductsFixedOrder(Reactor reactor, Molecule[] substrates, Integer maxProducts) throws ReactionException { reactor.setReactants(substrates); List<Molecule[]> results = new ArrayList<>(); Molecule[] products; while ((products = reactor.react()) != null) { results.add(products); if (maxProducts != null && maxProducts > 0 && results.size() >= maxProducts) { break; } } return results; } public List<Molecule[]> getProductsFixedOrder(Reactor reactor, Molecule[] substrates) throws ReactionException { return getProductsFixedOrder(reactor, substrates, null); } /** * Builds a string which will be identical for two molecule arrays which represent the same molecules * in the same order, and different otherwise. * * @param mols The array of molecules. * @return The string representation. * @throws IOException */ private String getStringHash(Molecule[] mols) throws IOException { StringBuilder builder = new StringBuilder(); for (Molecule molecule : mols) { builder.append(getMoleculeString(molecule)); } return builder.toString(); } /** * Gets a a string of the molecule. The exporter takes care of all caching so we don't need to worry about it * * @param molecule The molecule. * @return The molecule's string presentation in the format that this class declares.. * @throws IOException */ private String getMoleculeString(Molecule molecule) throws IOException { return MoleculeExporter.exportMolecule(molecule, MoleculeFormat$.MODULE$.getName(this.moleculeFormat)); } }