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 org.twentyn.proteintodna; import act.server.MongoDB; import act.shared.Chemical; import act.shared.Reaction; import act.shared.Seq; import com.act.biointerpretation.metadata.ProteinMetadata; import com.act.biointerpretation.metadata.RankPathway; import com.act.reachables.Cascade; import com.act.reachables.ReactionPath; import com.act.utils.CLIUtil; import com.mongodb.BasicDBObject; import com.mongodb.DB; import com.mongodb.DBObject; import com.mongodb.MongoClient; import com.mongodb.MongoInternalException; import com.mongodb.ServerAddress; import org.apache.commons.cli.CommandLine; import org.apache.commons.cli.Option; import org.apache.commons.lang3.tuple.Pair; import org.apache.logging.log4j.LogManager; import org.apache.logging.log4j.Logger; import org.json.JSONArray; import org.json.JSONObject; import org.mongojack.DBCursor; import org.mongojack.DBQuery; import org.mongojack.JacksonDBCollection; import org.mongojack.WriteResult; import java.util.ArrayList; import java.util.Arrays; import java.util.Collections; import java.util.HashMap; import java.util.HashSet; import java.util.List; import java.util.Map; import java.util.Set; import java.util.regex.Pattern; import java.util.stream.Collectors; public class ProteinToDNADriver { private static final Logger LOGGER = LogManager.getFormatterLogger(ProteinToDNADriver.class); private static final String DEFAULT_DB_HOST = "localhost"; private static final String DEFAULT_DB_PORT = "27017"; private static final String DEFAULT_OUTPUT_DB_NAME = "wiki_reachables"; private static final String DEFAULT_INPUT_DB_NAME = "SHOULD_COME_FROM_CMDLINE"; // "jarvis_2016-12-09"; public static final String DEFAULT_INPUT_PATHWAY_COLLECTION_NAME = "SHOULD_COME_FROM_CMDLINE"; // "pathways_jarvis"; public static final String DEFAULT_OUTPUT_PATHWAY_COLLECTION_NAME = "SHOULD_COME_FROM_CMDLINE"; // "pathways_vijay"; public static final String DEFAULT_OUTPUT_DNA_SEQ_COLLECTION_NAME = "SHOULD_COME_FROM_CMDLINE"; // "dna_designs"; private static final String OPTION_DB_HOST = "H"; private static final String OPTION_DB_PORT = "p"; private static final String OPTION_OUTPUT_DB_NAME = "o"; private static final String OPTION_INPUT_DB_NAME = "i"; private static final String OPTION_INPUT_PATHWAY_COLLECTION_NAME = "c"; private static final String OPTION_OUTPUT_PATHWAY_COLLECTION_NAME = "d"; private static final String OPTION_OUTPUT_DNA_SEQ_COLLECTION_NAME = "e"; private static final String OPTION_DESIGN_SOME = "m"; private static final Integer HIGHEST_SCORING_INFERRED_SEQ_INDEX = 0; private static final Set<String> BLACKLISTED_WORDS_IN_INFERRED_SEQ = new HashSet<>(Arrays.asList("Fragment")); private static final Pattern REGEX_ID = Pattern.compile("^\\d+$"); private static final Pattern REGEX_INCHI = Pattern.compile("^InChI=1S?/"); // Based on https://en.wikipedia.org/wiki/International_Chemical_Identifier#InChIKey private static final Pattern REGEX_INCHI_KEY = Pattern.compile("^[A-Z]{14}-[A-Z]{10}-[A-Z]$"); public static final List<Option.Builder> OPTION_BUILDERS = new ArrayList<Option.Builder>() { { add(Option.builder(OPTION_DB_HOST).argName("hostname") .desc(String.format("The DB host to which to connect (default: %s)", DEFAULT_DB_HOST)).hasArg() .longOpt("db-host")); add(Option.builder(OPTION_DB_PORT).argName("port") .desc(String.format("The DB port to which to connect (default: %s)", DEFAULT_DB_PORT)).hasArg() .longOpt("db-port")); add(Option.builder(OPTION_OUTPUT_DB_NAME).argName("output-db-name") .desc(String.format("The name of the database to read pathways and write seqs to (default: %s)", DEFAULT_OUTPUT_DB_NAME)) .hasArg().longOpt("output-db-name")); add(Option.builder(OPTION_INPUT_DB_NAME).argName("input-db-name").desc(String .format("The name of the database to read reactions from (default: %s)", DEFAULT_INPUT_DB_NAME)) .hasArg().required().longOpt("source-db-name")); add(Option.builder(OPTION_INPUT_PATHWAY_COLLECTION_NAME).argName("collection-name") .desc(String.format("The name of the input pathway collection to read from (default: %s)", DEFAULT_INPUT_PATHWAY_COLLECTION_NAME)) .hasArg().required().longOpt("input-pathway-collection")); add(Option.builder(OPTION_OUTPUT_PATHWAY_COLLECTION_NAME).argName("collection-name") .desc(String.format("The name of the output pathway collection to write to (default: %s)", DEFAULT_OUTPUT_PATHWAY_COLLECTION_NAME)) .hasArg().required().longOpt("output-pathway-collection")); add(Option.builder(OPTION_OUTPUT_DNA_SEQ_COLLECTION_NAME).argName("collection-name") .desc(String.format("The name of the output dna seq collection to write to (default: %s)", DEFAULT_OUTPUT_DNA_SEQ_COLLECTION_NAME)) .hasArg().required().longOpt("output-dna-seq-collection")); add(Option.builder(OPTION_DESIGN_SOME).argName("molecule").desc( "Generate designs for specified molecules; can be a numeric ids, InChIs, or InChI Keys, *separated by '|'* " + "for InChI compatibility. Molecules must exist in chemical source DB if InChI or InChI Key is used.") .hasArgs().valueSeparator('|').longOpt("design-this")); } }; public static final String HELP_MESSAGE = "This class is the driver to extract protein sequences from pathways and construct DNA designs from these proteins."; private static final CLIUtil CLI_UTIL = new CLIUtil(ProteinToDNADriver.class, HELP_MESSAGE, OPTION_BUILDERS); /** * This function gets all protein combinations of a pathway from candidate protein sequences from each reaction on * the pathway. * @param listOfSetOfProteinSequences A list of sets of candidate protein sequences in the pathway * @return A set of all possible combinations of proteins from all the reactions in the pathway. */ public static Set<List<String>> makePermutations(List<Set<String>> listOfSetOfProteinSequences) { Set<List<String>> accum = new HashSet<>(); makePermutationsHelper(accum, listOfSetOfProteinSequences, new ArrayList<>()); return accum; } private static void makePermutationsHelper(Set<List<String>> accum, List<Set<String>> input, List<String> prefix) { // Base case: no more sequences to add. Accumulate and return. if (input.isEmpty()) { accum.add(prefix); return; } // Recursive case: iterate through next level of input sequences, appending each to a prefix and recurring. Set<String> head = input.get(0); // Avoid index out of bounds exception. List<Set<String>> rest = input.size() > 1 ? input.subList(1, input.size()) : Collections.emptyList(); for (String next : head) { List<String> newPrefix = new ArrayList<>(prefix); newPrefix.add(next); makePermutationsHelper(accum, rest, newPrefix); } } public static void main(String[] args) throws Exception { CommandLine cl = CLI_UTIL.parseCommandLine(args); String reactionDbName = cl.getOptionValue(OPTION_INPUT_DB_NAME, DEFAULT_INPUT_DB_NAME); String dbHost = cl.getOptionValue(OPTION_DB_HOST, DEFAULT_DB_HOST); Integer dbPort = Integer.valueOf(cl.getOptionValue(OPTION_DB_PORT, DEFAULT_DB_PORT)); MongoDB reactionDB = new MongoDB(dbHost, dbPort, reactionDbName); MongoClient inputClient = new MongoClient(new ServerAddress(dbHost, dbPort)); String outDB = cl.getOptionValue(OPTION_OUTPUT_DB_NAME, DEFAULT_OUTPUT_DB_NAME); DB db = inputClient.getDB(outDB); String inputPathwaysCollectionName = cl.getOptionValue(OPTION_INPUT_PATHWAY_COLLECTION_NAME, DEFAULT_INPUT_PATHWAY_COLLECTION_NAME); String outputPathwaysCollectionName = cl.getOptionValue(OPTION_OUTPUT_PATHWAY_COLLECTION_NAME, DEFAULT_OUTPUT_PATHWAY_COLLECTION_NAME); String outputDnaDeqCollectionName = cl.getOptionValue(OPTION_OUTPUT_DNA_SEQ_COLLECTION_NAME, DEFAULT_OUTPUT_DNA_SEQ_COLLECTION_NAME); JacksonDBCollection<ReactionPath, String> inputPathwayCollection = JacksonDBCollection .wrap(db.getCollection(inputPathwaysCollectionName), ReactionPath.class, String.class); JacksonDBCollection<DNADesign, String> dnaDesignCollection = JacksonDBCollection .wrap(db.getCollection(outputDnaDeqCollectionName), DNADesign.class, String.class); JacksonDBCollection<ReactionPath, String> outputPathwayCollection = JacksonDBCollection .wrap(db.getCollection(outputPathwaysCollectionName), ReactionPath.class, String.class); Map<String, Set<ProteinInformation>> proteinSeqToOrgInfo = new HashMap<>(); ProteinsToDNA2 p2d = ProteinsToDNA2.initiate(); List<Long> reachableIds = cl.hasOption(OPTION_DESIGN_SOME) ? Arrays.stream(cl.getOptionValues(OPTION_DESIGN_SOME)).map(m -> lookupMolecule(reactionDB, m)) .collect(Collectors.toList()) : Collections.emptyList(); // If there is a list of molecules for which to create designs, only consider their pathways. Otherwise, find all. DBObject query = reachableIds.isEmpty() ? new BasicDBObject() : DBQuery.in("target", reachableIds); /* Extract all pathway ids, then read pathways one id at a time. This will reduce the likelihood of the cursor * timing out, which has a tendency to happen when doing expensive operations before advancing (as done here). */ DBCursor<ReactionPath> cursor = inputPathwayCollection.find(query, new BasicDBObject("_id", true)); List<String> pathwayIds = new ArrayList<>(); while (cursor.hasNext()) { pathwayIds.add(cursor.next().getId()); } for (String pathwayId : pathwayIds) { ReactionPath reactionPath = inputPathwayCollection.findOne(DBQuery.is("_id", pathwayId)); List<List<Pair<ProteinMetadata, Integer>>> processedP = RankPathway .processSinglePathAsJava(reactionPath, reactionDbName, outDB); if (processedP == null) { LOGGER.info(String.format( "Process pathway was filtered out possibly because there were more than %s seqs for a given pathway", RankPathway.MAX_PROTEINS_PER_PATH())); outputPathwayCollection.insert(reactionPath); continue; } Boolean atleastOneSeqMissingInPathway = false; List<Set<String>> proteinPaths = new ArrayList<>(); for (Cascade.NodeInformation nodeInformation : reactionPath.getPath().stream() .filter(nodeInfo -> nodeInfo.getIsReaction()).collect(Collectors.toList())) { Set<String> proteinSeqs = new HashSet<>(); for (Long id : nodeInformation.getReactionIds()) { // If the id is negative, it is a reaction in the reverse direction. Moreover, the enzyme for this reverse // reaction is the same, so can use the actual positive reaction id's protein seq reference. // TODO: Add a preference for the positive forward direction compared to the negative backward direction seq. if (id < 0) { LOGGER.info("Found a negative reaction id", id); id = Reaction.reverseID(id); } Reaction reaction = reactionDB.getReactionFromUUID(id); for (JSONObject data : reaction.getProteinData()) { // Get the sequences if (data.has("sequences")) { JSONArray seqs = data.getJSONArray("sequences"); for (int i = 0; i < seqs.length(); i++) { Long s = seqs.getLong(i); if (s != null) { Seq sequenceInfo = reactionDB.getSeqFromID(s); String dnaSeq = sequenceInfo.getSequence(); if (dnaSeq == null) { LOGGER.info(String.format( "Sequence string for seq id %d, reaction id %d and reaction path %s is null", s, id, reactionPath.getId())); continue; } // odd sequence if (dnaSeq.length() <= 80 || dnaSeq.charAt(0) != 'M') { JSONObject metadata = sequenceInfo.getMetadata(); if (!metadata.has("inferred_sequences") || metadata.getJSONArray("inferred_sequences").length() == 0) { continue; } JSONArray inferredSequences = metadata.getJSONArray("inferred_sequences"); // get the first inferred sequence since it has the highest hmmer score JSONObject object = inferredSequences .getJSONObject(HIGHEST_SCORING_INFERRED_SEQ_INDEX); for (String blacklistWord : BLACKLISTED_WORDS_IN_INFERRED_SEQ) { if (object.getString("fasta_header").contains(blacklistWord)) { continue; } } dnaSeq = object.getString("sequence"); } proteinSeqs.add(dnaSeq); ProteinInformation proteinInformation = new ProteinInformation( sequenceInfo.getOrgName(), sequenceInfo.getEc(), sequenceInfo.getSequence(), reaction.getReactionName()); if (!proteinSeqToOrgInfo.containsKey(dnaSeq)) { proteinSeqToOrgInfo.put(dnaSeq, new HashSet<>()); } proteinSeqToOrgInfo.get(dnaSeq).add(proteinInformation); } } } } } if (proteinSeqs.size() == 0) { LOGGER.info("The reaction does not have any viable protein sequences"); atleastOneSeqMissingInPathway = true; break; } // Now we select two representative protein seqs from the reaction. In order to do this deterministically, // we sort and pick the first and middle index protein seqs. List<String> proteinSeqArray = new ArrayList<>(proteinSeqs); Collections.sort(proteinSeqArray); int firstIndex = 0; int middleIndex = proteinSeqs.size() / 2; // get first seq Set<String> combination = new HashSet<>(); combination.add(proteinSeqArray.get(firstIndex)); // get middle index of the protein seq array if (proteinSeqs.size() > 1) { combination.add(proteinSeqArray.get(middleIndex)); } proteinPaths.add(combination); } if (atleastOneSeqMissingInPathway) { LOGGER.info(String.format("There is at least one reaction with no sequence in reaction path id: %s", reactionPath.getId())); } else { LOGGER.info(String.format("All reactions in reaction path have at least one viable seq: %s", reactionPath.getId())); // We only compute the dna design if we can find at least one sequence for each reaction in the pathway. Set<List<String>> pathwayProteinCombinations = makePermutations(proteinPaths); Set<DNAOrgECNum> dnaDesigns = new HashSet<>(); for (List<String> proteinsInPathway : pathwayProteinCombinations) { try { Construct dna = p2d.computeDNA(proteinsInPathway, Host.Ecoli); List<Set<ProteinInformation>> seqMetadata = new ArrayList<>(); for (String protein : proteinsInPathway) { seqMetadata.add(proteinSeqToOrgInfo.get(protein)); } DNAOrgECNum instance = new DNAOrgECNum(dna.toSeq(), seqMetadata, proteinsInPathway.size()); dnaDesigns.add(instance); } catch (Exception ex) { LOGGER.error("The error thrown while trying to call computeDNA", ex.getMessage()); } } DNADesign dnaDesignSeq = new DNADesign(dnaDesigns); try { WriteResult<DNADesign, String> result = dnaDesignCollection.insert(dnaDesignSeq); String id = result.getSavedId(); reactionPath.setDnaDesignRef(id); } catch (MongoInternalException e) { // This condition happens whent the protein designs are too big and cannot be inserted in to the JSON object. // TODO: Handle this case without dropping the record LOGGER.error(String.format("Mongo internal exception caught while inserting dna design: %s", e.getMessage())); } } outputPathwayCollection.insert(reactionPath); } } private static Long lookupMolecule(MongoDB db, String someKey) { if (REGEX_ID.matcher(someKey).find()) { // Note: this doesn't verify that the chemical id is valid. Maybe we should do that? return Long.valueOf(someKey); } else if (REGEX_INCHI.matcher(someKey).find()) { Chemical c = db.getChemicalFromInChI(someKey); if (c != null) { return c.getUuid(); } } else if (REGEX_INCHI_KEY.matcher(someKey).find()) { Chemical c = db.getChemicalFromInChIKey(someKey); if (c != null) { return c.getUuid(); } } else { String msg = String.format("Unable to find key type for query '%s'", someKey); LOGGER.error(msg); throw new IllegalArgumentException(msg); } String msg = String.format("Unable to find matching chemical for query %s", someKey); LOGGER.error(msg); throw new IllegalArgumentException(msg); } }