use of cbit.vcell.bionetgen.BNGComplexSpecies in project vcell by virtualcell.
the class BNGExecutorServiceMultipass method doWork.
private CorrectedSR doWork(String oldSpeciesString, String newNetFile) throws ParseException {
// TODO: make a map to preserve the ancestry of each generated species (rule and iteration that generated it)
// each species may be generated multiple times, by different rules, at different iterations
// the key should be the normalized (sorted lexicographically) expression of each species
// parse the .net file with BNGOutputFileParser
// seed species at the beginning of the current iteration
List<BNGSpecies> oldSpeciesList = BNGOutputFileParser.createBngSpeciesOutputSpec(oldSpeciesString);
// same as above, used internally, as we add new species from the current iteration we put them here;
// used to speed up the verifications and reduce the isomorphism calls together with the isomorphismSignaturesMap
Map<String, BNGSpecies> seedSpeciesMap = BNGOutputFileParser.createBngSpeciesSignatureMap(oldSpeciesString);
// .net file content generated during current iteration
BNGOutputSpec workSpec = BNGOutputFileParser.createBngOutputSpec(newNetFile);
// we build here the list of valid (perhaps even corrected) NEW species
List<BNGSpecies> newSpeciesList = new ArrayList<>();
// we build here the list of flattened reactions (corrected at need, if the species were corrected)
ArrayList<BNGReaction> newReactionsList = new ArrayList<BNGReaction>(Arrays.asList(workSpec.getBNGReactions()));
// ArrayList<ObservableGroup> newObservablesList = new ArrayList<ObservableGroup>(); // here we put all the observables after correction
// ArrayList<ObservableGroup> newObservablesList = new ArrayList<ObservableGroup>(Arrays.asList(workSpec.getObservableGroups()));
// identify new products, loop thorough each of them
Pair<List<BNGSpecies>, List<BNGSpecies>> p = BNGSpecies.diff(oldSpeciesList, workSpec.getBNGSpecies());
List<BNGSpecies> removed = p.one;
// after possible needed corrections we may end up with more or less new species for this iteration
List<BNGSpecies> added = p.two;
if (!removed.isEmpty()) {
// this may actually be fired if the BNGSpecies.equals() fails for some reason
throw new RuntimeException("No existing BNGSpecies may be removed: \n" + removed.get(0));
}
if (added.isEmpty()) {
System.out.println("No new species added this iteration. Finished!");
}
// in this map we store which of the indexes of species generated during the current iteration
// were converted into what after correction
// an index may exist already in which case we won't save a duplicate or the same new species (that needs repair)
// may result in more new species after being fixed
// param1: original index of the newly created species before repairing it
// param2: list of the indexes of the species resulting from the original after repairing
Map<String, List<String>> indexesMap = new LinkedHashMap<>();
// some reactions may end up as identity reactions after corrections caused by anchoring
// for instance a transport xA -> yA may become xA -> xA after correction if A is anchored to x
// we make a list of these and get rid of them at the end
Set<BNGReaction> identityReactions = new HashSet<>();
// as we validate and we add new species, we use this index to set their network index
// we can't get the indexes given to the newly generated species for granted, during correction we may
// lose or gain species
int firstAvailableIndex = BNGOutputSpec.getFirstAvailableSpeciesIndex(oldSpeciesList);
int summaryCreated = added.size();
int summaryRepaired = 0;
int summaryExisted = 0;
// species we're adding at the end of this iteration, after repairing and checking for existing
int summarySpecies = 0;
for (BNGSpecies s : added) {
// find the flattened reaction and from the reaction find the rule it's coming from
for (BNGReaction r : newReactionsList) {
// console only message
String message = "";
// array list with position of this species in the product (may be more than 1)
List<Integer> reactionProductPositions = r.findProductPositions(s);
if (reactionProductPositions.isEmpty()) {
// this species is not a product of this reaction, nothing to do
continue;
}
String ruleName = r.getRuleName();
if (r.isRuleReversed()) {
throw new RuntimeException("Reaction " + ruleName + " cannot be 'reverse'.");
}
// System.out.println("Species " + s + " found in rule " + r.getRuleName() + ", at index " + position);
// check the product against the rule to see if it's valid
// sanity check: only "transport" rules can give incorrect products, any rule with all participants in the same
// compartment should only give valid products (is that so?)
String structureName;
if (model.getStructures().length == 1) {
structureName = model.getStructure(0).getName();
} else {
// this is the structure where the product should be, according to the rule
structureName = findRuleProductCompartment(s);
}
ReactionRule rr = model.getRbmModelContainer().getReactionRule(ruleName);
// TODO: the code below may be greatly simplified using the more advanced BNGSpecies classes instead of using the strings
// 'one' is the list of all the compartments mentioned in this product; for single compartment models 'one' will always be 1
// and thus we'll never attempt to correct anything
// 'two' is the BNGSpecies string with the compartment info extracted (that is, the AAA sites extracted)
Pair<List<String>, String> pair = RbmUtils.extractCompartment(s.getName());
boolean needsRepairing = false;
if (pair.one.size() > 1) {
// System.out.println(s.getName() + " multiple compartments, needs repairing.");
message += s.getName() + " needs repairing... ";
needsRepairing = true;
/* At this point we have the 'probable' structure of this product, extracted from the rule (structureName)
However, we need to verify against the lists of anchors to make sure is not excluded from the possible list of anchors
for any of the molecules and choose a compatible structure name, as follows:
compartments x y z
molecules A B C A comes as transport
rule says that this product belongs in x
cases:
1. A~y xB.xC.yA -> yB.yC.yA - can't use x as the rule says because A must be anchored to y and that takes precedence
- if B or C can't be in y throw exception, our anchors are too restrictive
2. A~x~y xB.xC.yA -> xB.xC.xA - correct to what the rule wants since x is acceptable as anchor for A
3. A~z xB.xC.yA -> !!! - impossible, there is no reactant such that A can get in y (throw exception)
4. A~x~z xB.xC.yA -> !!! - same as above
step 1. look for exceptions: make sets with the compartments for each molecule in the species and see if
they all are allowed by the anchor rules (if not, we are in cases 3 or 4 - throw exception
step 2. see if structureName is acceptable for all anchored molecules - if yes, use it
step 3. intersect all sets from step 1
3.1 - if there's 1 compartment that is acceptable for all, use it
3.2 - if there's no compartment then we can't put the species anywhere - throw exception
3.3 - if there's more than 1 compartment, it's ambiguous - throw exception
*/
// ex of one entry: key T, value mem, cyt
Map<String, Set<String>> speciesAnchorMap = extractSpeciesAnchorMap(s);
// sanity check, step 1
for (Map.Entry<String, Set<String>> entry : speciesAnchorMap.entrySet()) {
// molecule in the newly generated species (product)
String molecule = entry.getKey();
if (!anchorsMap.containsKey(molecule)) {
// this molecule can be anywhere, it's not anchored - nothing to verify
continue;
}
// 1 or more location where we try to put it (before correction may be more than 1)
Set<String> structures = entry.getValue();
// these structures MUST all be allowed by the anchor rule, otherwise it would mean they are coming from a reactant that cannot exist (case 3 or 4 above)
for (String wantedStructure : structures) {
Set<String> allowedStructures = anchorsMap.get(molecule);
if (!allowedStructures.contains(wantedStructure)) {
// should never happen: no reactant and no rule should have placed this molecule in this structure
throw new RuntimeException("Error in " + s.getName() + "\nMolecule " + molecule + " cannot possibly be in Structure " + wantedStructure);
}
}
}
// step 2
boolean isCandidateStructureAcceptable = true;
for (Map.Entry<String, Set<String>> entry : speciesAnchorMap.entrySet()) {
// molecule in the newly generated species (product)
String molecule = entry.getKey();
if (!anchorsMap.containsKey(molecule)) {
// this molecule can be anywhere, it's not anchored - nothing to verify
continue;
}
Set<String> allowedStructures = anchorsMap.get(molecule);
if (!allowedStructures.contains(structureName)) {
// the structure 'wanted' by the rule for this product is not acceptable because of some anchor rule
isCandidateStructureAcceptable = false;
break;
}
}
// step 3
if (isCandidateStructureAcceptable == false) {
// we try to find one and only one structure that would satisfy all anchored molecules of this species
TreeSet<String> intersection = new TreeSet<>();
for (Structure str : model.getStructures()) {
// we start by adding all structures in the model
intersection.add(str.getName());
}
for (Map.Entry<String, Set<String>> entry : speciesAnchorMap.entrySet()) {
// molecule in the newly generated species (product)
String molecule = entry.getKey();
if (!anchorsMap.containsKey(molecule)) {
// this molecule can be anywhere, it's not anchored - nothing to verify
continue;
}
// we already know from step 1 that all structures here are acceptable by the anchor
Set<String> structures = anchorsMap.get(molecule);
// intersection retains only the elements in 'structures'
intersection.retainAll(structures);
}
if (intersection.size() == 0) {
// no structure satisfies the anchor rules for the molecules in this species
throw new RuntimeException("Error in " + s.getName() + "\nThe Structures allowed for anchored molecules are mutually exclusive.");
} else if (intersection.size() > 1) {
// ambiguous, don't know which compartment to pick
throw new RuntimeException("Error in " + s.getName() + "\nMultiple Structures allowed by the anchor rules, don't know which to choose.");
} else {
// found one single structure everybody is happy about
structureName = intersection.first();
}
}
} else {
// no transport caused by a 'wild card', just direct rule application (no compartment ambiguity, nothing to correct)
String structure;
if (model.getStructures().length == 1) {
structure = model.getStructure(0).getName();
} else {
structure = pair.one.get(0);
}
if (!structure.equals(structureName)) {
// the rule (no transport caused by '?' is possible)
throw new RuntimeException("Error in " + s.getName() + "\nIf one single structure is present in the species it must match the structure of the rule product");
}
// product structure (from rule) should not conflict with the anchor - we should have error issue in the rule itself reporting the conflict
// TODO: check here too anyway, just in case?
}
// new generated species that we may add the list of species for the next iteration
BNGSpecies candidate;
// deleted or added previously during earlier iterations through newly created species
if (needsRepairing) {
// if not valid, correct the fake "compartment" site to be conform to the compartment of
// the product pattern
String speciesRepairedName = RbmUtils.repairCompartment(s.getName(), structureName);
speciesRepairedName = RbmUtils.resetProductIndex(speciesRepairedName);
candidate = new BNGComplexSpecies(speciesRepairedName, s.getConcentration(), firstAvailableIndex);
// System.out.println(candidate.getName() + " repaired!");
message += "repaired from rule " + rr.getDisplayName() + "... ";
summaryRepaired++;
System.out.println("");
} else {
String speciesName = s.getName();
speciesName = RbmUtils.resetProductIndex(speciesName);
candidate = new BNGComplexSpecies(speciesName, s.getConcentration(), firstAvailableIndex);
}
// At this point we have a valid candidate - but we may not need it if it already exist in the list of old seed species or in the list of
// new seed species we're building now (it may exist in either if correction took place)
//
// We correct the reaction to match the networkFileIndex of this species (useless activity except for the last iteration
// when we want the valid list of flattened reactions
//
// we set this to an existing species if the candidate matches it (directly or through isomorphism)
BNGSpecies existingMatch = null;
long st = System.currentTimeMillis();
String sig = BNGSpecies.getShortSignature(candidate, sigDetailLevel);
if (!shortSignaturesSet.contains(sig)) {
// our candidate has a new short signature, this means it is certainly not among the existing species
// no point in doing any other search to find it
} else {
// short signature match, now we check if the full species exists or is isomorphic with an existing one
String isomorphSpeciesName = isomorphismSignaturesMap.get(candidate.getName());
if (isomorphSpeciesName != null) {
// name of an existing species that is isomorphic with our candidate
// we recover the existing species and will use it instead of the candidate
existingMatch = seedSpeciesMap.get(isomorphSpeciesName);
if (existingMatch == null) {
throw new RuntimeException("Unable to find an isomorph BNGSpecies named " + isomorphSpeciesName + " that should have existed.");
}
} else {
// not present among the existing species or their known isomorphisms
// try to find an isomorphism the hard way
BNGSpecies existingMatchInNew = BNGOutputSpec.findMatch(candidate, newSpeciesList, sigDetailLevel);
if (existingMatchInNew != null) {
isomorphismSignaturesMap.put(candidate.getName(), existingMatchInNew.getName());
existingMatch = existingMatchInNew;
} else {
BNGSpecies existingMatchInOld = BNGOutputSpec.findMatch(candidate, oldSpeciesList, sigDetailLevel);
if (existingMatchInOld != null) {
isomorphismSignaturesMap.put(candidate.getName(), existingMatchInOld.getName());
existingMatch = existingMatchInOld;
}
}
}
}
long et = System.currentTimeMillis();
eltIsomorph += (et - st);
// if it doesn't exist we add it to the list of seed species we are preparing for the next iteration
if (existingMatch == null) {
newSpeciesList.add(candidate);
message += "Candidate " + candidate.getName() + " added to the seed species list.";
summarySpecies++;
firstAvailableIndex++;
for (int reactionProductPosition : reactionProductPositions) {
// correct the reaction
r.getProducts()[reactionProductPosition] = candidate;
}
manageIndexesMap(indexesMap, s, candidate);
seedSpeciesMap.put(candidate.getName(), candidate);
// put self in the list of signature map too
isomorphismSignaturesMap.put(candidate.getName(), candidate.getName());
String signature = BNGSpecies.getShortSignature(candidate, sigDetailLevel);
shortSignaturesSet.add(sig);
} else {
message += "Candidate " + candidate.getName() + " already exists, not added.";
summaryExisted++;
for (int reactionProductPosition : reactionProductPositions) {
// correct the reaction
r.getProducts()[reactionProductPosition] = existingMatch;
}
manageIndexesMap(indexesMap, s, existingMatch);
}
if (needsRepairing && existingMatch != null) {
if (isIdentityReaction(r)) {
// System.out.println("Identity reaction detected after product correction: " + r);
identityReactions.add(r);
}
}
}
// end checking reactions for this species
}
// we'll get rid of these
if (!identityReactions.isEmpty()) {
// we found some
newReactionsList.removeAll(identityReactions);
}
// System.out.println("------------- Finished checking newly generated species for this iteration. Summary:");
System.out.println(" Added " + newSpeciesList.size() + " new species");
System.out.println(" ");
// String summary = " " + summaryCreated + " species were created.\n";
// summary += " " + summaryRepaired + " species needed repairing.\n";
// summary += " " + summaryExisted + " species already present, not added to the seed species list.\n";
// summary += " " + summarySpecies + " new species added.\n";
// consoleNotification(summary);
// total number of species at this moment
currentIterationTotalSpecies = previousIterationTotalSpecies + summarySpecies;
String message = previousIterationTotalSpecies + "," + currentIterationTotalSpecies + "," + needAdjustMaxMolecules;
TaskCallbackMessage newCallbackMessage = new TaskCallbackMessage(TaskCallbackStatus.AdjustAllFlags, message);
broadcastCallbackMessage(newCallbackMessage);
previousIterationTotalSpecies = currentIterationTotalSpecies;
CorrectedSR sr = new CorrectedSR(newSpeciesList, newReactionsList);
return sr;
}
use of cbit.vcell.bionetgen.BNGComplexSpecies in project vcell by virtualcell.
the class BNGExecutorServiceMultipass method findRuleProductCompartment.
private static String findRuleProductCompartment(BNGSpecies s) {
boolean needsCorrection = false;
// the structures we encounter in this species (use tree set to avoid duplicates)
Set<String> structureNameSet = new TreeSet<>();
// if it's a complex species we break it down in simple ones and populate this list
List<BNGSpecies> simpleSpeciesList = new ArrayList<>();
if (s instanceof BNGComplexSpecies) {
simpleSpeciesList.addAll(Arrays.asList(s.parseBNGSpeciesName()));
} else {
// if it's a simple species to begin with we'll only have one element in this list
simpleSpeciesList.add(s);
}
for (BNGSpecies simpleSpecies : simpleSpeciesList) {
// must be multi state (actually multi-site!), we have at least 2 components (sites): RbmUtils.SiteStruct and RbmUtils.SiteProduct
if (!(simpleSpecies instanceof BNGMultiStateSpecies)) {
throw new RuntimeException("Species " + s.getName() + " must be instance of BNGMultiStateSpecies");
}
BNGMultiStateSpecies mss = (BNGMultiStateSpecies) simpleSpecies;
// System.out.println(" " + mss.toString());
List<BNGSpeciesComponent> componentsList = new ArrayList<>(Arrays.asList(mss.getComponents()));
String structName = "";
String prodPosition = "";
for (BNGSpeciesComponent sc : componentsList) {
// System.out.println(" " + sc.getComponentName() + "~" + sc.getCurrentState());
if (sc.getComponentName().equals(RbmUtils.SiteStruct)) {
structName = sc.getCurrentState();
} else if (sc.getComponentName().equals(RbmUtils.SiteProduct)) {
prodPosition = sc.getCurrentState();
}
}
if (structName == null || structName.isEmpty() || structName.equals("null")) {
throw new RuntimeException("Structure name site missing from BNGSpecies " + mss);
}
if (prodPosition == null || prodPosition.isEmpty() || prodPosition.equals("null")) {
throw new RuntimeException("Index of original rule product site missing from BNGSpecies " + mss);
}
switch(prodPosition) {
case // comes from matching a wildcard with a seed species, its compartment may be wrong
"0":
break;
case // comes by directly applying a rule, this is the correct compartment inherited from the rule product
"1":
structureNameSet.add(structName);
break;
default:
throw new RuntimeException("Product position must be 0 or 1 for " + mss.getName());
}
}
if (structureNameSet.size() != 1) {
throw new RuntimeException("The number of structures for " + s.getName() + " must be exactly 1.");
}
String structName = structureNameSet.iterator().next();
return structName;
}
use of cbit.vcell.bionetgen.BNGComplexSpecies in project vcell by virtualcell.
the class BNGExecutorServiceMultipass method extractPolymerObservablesAsString.
private String extractPolymerObservablesAsString(String prefix, String sBngInputString) {
if (polymerEqualObservables.isEmpty() && polymerGreaterObservables.isEmpty()) {
return prefix;
}
String observablesString = "";
List<BNGSpecies> bngSpeciesList = BNGOutputFileParser.createBngSpeciesOutputSpec(sBngInputString);
// ordered by the network file index which is also key
Map<Integer, Map<String, Integer>> masterSignaturesMap = new LinkedHashMap<>();
// key = network file index, value = compartment of species
Map<Integer, String> masterCompartmentMap = new LinkedHashMap<>();
for (BNGSpecies species : bngSpeciesList) {
List<BNGSpecies> ourList = new ArrayList<>();
if (species instanceof BNGComplexSpecies) {
ourList.addAll(Arrays.asList(species.parseBNGSpeciesName()));
} else {
// if it's a simple species to begin with we'll only have one element in this list
ourList.add(species);
}
// key = name of molecules, value = number of occurrences
Map<String, Integer> signatureMap = new HashMap<>();
for (BNGSpecies mtp : ourList) {
if (!(mtp instanceof BNGMultiStateSpecies)) {
throw new RuntimeException("Species " + mtp.getName() + " must be instance of BNGMultiStateSpecies");
}
BNGMultiStateSpecies ss = (BNGMultiStateSpecies) mtp;
String mtpName = ss.extractMolecularTypeName();
if (signatureMap.containsKey(mtpName)) {
int value = signatureMap.get(mtpName);
value += 1;
signatureMap.put(mtpName, value);
} else {
signatureMap.put(mtpName, 1);
}
}
int networkFileIndex = species.getNetworkFileIndex();
masterSignaturesMap.put(networkFileIndex, signatureMap);
// we look in first mtp of this seed species, the compartment is the same in all
BNGMultiStateSpecies ss = (BNGMultiStateSpecies) ourList.get(0);
String compartment = ss.extractCompartment();
if (compartment == null) {
// single compartment, we get it from the model
if (model.getStructures().length > 1) {
throw new RuntimeException("BNGExecutorServiceMultipass: Unable to extract compartment");
}
compartment = model.getStructure(0).getName();
}
masterCompartmentMap.put(networkFileIndex, compartment);
}
// we need to find the next available index for observables
int nextAvailableIndex;
if (prefix != null) {
nextAvailableIndex = extractNextAvailableIndexFromObservables(prefix);
} else {
nextAvailableIndex = 1;
}
// we assume this is an observable made of only 1 sp that contains only one mtp, ex: A()=xx where xx is number of occurrences of A
for (RbmObservable oo : polymerEqualObservables) {
// for each polymer observable, here we build the string with the network indexes of the species that satisfy the criteria
String speciesFoundIndexes = "";
String mtpName = oo.getSpeciesPatternList().get(0).getMolecularTypePatterns().get(0).getMolecularType().getDisplayName();
int sequenceLength = oo.getSequenceLength();
System.out.println(mtpName + "=" + sequenceLength);
// comma counter
int i = 0;
boolean found = false;
// check all seed species for those that satisfy this observable and build comma delimited string of network file indexes
for (Map.Entry<Integer, Map<String, Integer>> entry : masterSignaturesMap.entrySet()) {
Integer networkFileIndex = entry.getKey();
// only interested if compartment is the same for both obs and seed species
String compartment = masterCompartmentMap.get(networkFileIndex);
if (!compartment.equals(oo.getStructure().getName())) {
// seed species in other compartment than our observable
continue;
}
Map<String, Integer> value = entry.getValue();
if (value.containsKey(mtpName)) {
int occurences = value.get(mtpName);
if (sequenceLength == occurences) {
if (i > 0) {
speciesFoundIndexes += ",";
}
speciesFoundIndexes += networkFileIndex;
found = true;
i++;
}
}
}
if (found == false) {
// desired number of occurrences for the polymer species has not been found in any seed species
continue;
}
// finished for this observable
System.out.println("Observable " + oo.getDisplayName() + ": " + mtpName + "()=" + sequenceLength + " found in species " + speciesFoundIndexes + ".");
observablesString += "\t" + nextAvailableIndex + " " + oo.getDisplayName() + "\t" + speciesFoundIndexes + "\n";
nextAvailableIndex++;
}
for (RbmObservable oo : polymerGreaterObservables) {
// same as above for the A()>xx polymer observables
String speciesFoundIndexes = "";
String mtpName = oo.getSpeciesPatternList().get(0).getMolecularTypePatterns().get(0).getMolecularType().getDisplayName();
int sequenceLength = oo.getSequenceLength();
System.out.println(mtpName + "=" + sequenceLength);
int i = 0;
boolean found = false;
// check all seed species for those that satisfy this observable and build comma delimited string of network file indexes
for (Map.Entry<Integer, Map<String, Integer>> entry : masterSignaturesMap.entrySet()) {
Integer networkFileIndex = entry.getKey();
String compartment = masterCompartmentMap.get(networkFileIndex);
if (!compartment.equals(oo.getStructure().getName())) {
// seed species in other compartment than our observable
continue;
}
Map<String, Integer> value = entry.getValue();
if (value.containsKey(mtpName)) {
int occurences = value.get(mtpName);
if (sequenceLength < occurences) {
if (i > 0) {
speciesFoundIndexes += ",";
}
speciesFoundIndexes += networkFileIndex;
found = true;
i++;
}
}
}
if (found == false) {
continue;
}
System.out.println("Observable " + oo.getDisplayName() + ": " + mtpName + "()>" + sequenceLength + " found in species " + speciesFoundIndexes + ".");
observablesString += "\t" + nextAvailableIndex + " " + oo.getDisplayName() + "\t" + speciesFoundIndexes + "\n";
nextAvailableIndex++;
}
if (observablesString.isEmpty()) {
// may be null
return prefix;
}
if (prefix != null) {
observablesString = prefix + observablesString;
}
return observablesString;
}
Aggregations