Source code for reconstruction.ecoli.dataclasses.process.complexation

"""
SimulationData for the Complexation process
"""

import numpy as np
from wholecell.utils import units
from wholecell.utils.mc_complexation import mccBuildMatrices


[docs] class ComplexationError(Exception): pass
[docs] class MoleculeNotFoundError(ComplexationError): pass
[docs] class Complexation(object): def __init__(self, raw_data, sim_data): # Build the abstractions needed for complexation molecules = [] # List of all molecules involved in complexation subunits = [] # List of all molecules that participate as subunits complexes = [] # List of all molecules that participate as complexes stoichMatrixI = [] # Molecule indices stoichMatrixJ = [] # Reaction indices stoichMatrixV = [] # Stoichiometric coefficients stoichMatrixMass = [] # Molecular masses of molecules in stoichMatrixI self.ids_reactions = [] self.reaction_stoichiometry_unknown = [] reaction_index = 0 miscrnas_with_singleton_tus = sim_data.getter.get_miscrnas_with_singleton_tus() # Build stoichiometric matrix from given complexation reactions for reaction in raw_data.complexation_reactions: self.ids_reactions.append(reaction["id"]) stoichiometry_unknown = False for mol_id, coeff in reaction["stoichiometry"].items(): # Replace miscRNA subunit IDs with TU IDs if mol_id in miscrnas_with_singleton_tus: mol_id = sim_data.getter.get_singleton_tu_id(mol_id) mol_id_with_compartment = "{}[{}]".format( mol_id, sim_data.getter.get_compartment(mol_id)[0] ) if mol_id_with_compartment not in molecules: molecules.append(mol_id_with_compartment) molecule_index = len(molecules) - 1 else: molecule_index = molecules.index(mol_id_with_compartment) # Flag reactions whose stoichioemtric coefficients are given # as null and replace with -1 if coeff is None: stoichiometry_unknown = True coeff = -1 assert (coeff % 1) == 0 stoichMatrixI.append(molecule_index) stoichMatrixJ.append(reaction_index) stoichMatrixV.append(coeff) # Classify molecule into subunit or complex depending on sign # of the stoichiometric coefficient - Note that a molecule can # be both a subunit and a complex if coeff < 0: subunits.append(mol_id_with_compartment) else: complexes.append(mol_id_with_compartment) # Find molecular mass of the molecule and add to mass matrix molecularMass = sim_data.getter.get_mass( mol_id_with_compartment ).asNumber(units.g / units.mol) stoichMatrixMass.append(molecularMass) self.reaction_stoichiometry_unknown.append(stoichiometry_unknown) reaction_index += 1 self.rates = np.full( (reaction_index,), sim_data.constants.complexation_rate.asNumber(1 / units.s), ) self._stoich_matrix_I = np.array(stoichMatrixI) self._stoich_matrix_J = np.array(stoichMatrixJ) self._stoich_matrix_V = np.array(stoichMatrixV) self._stoich_matrix_mass = np.array(stoichMatrixMass) self.molecule_names = molecules self.ids_complexes = [ self.molecule_names[i] for i in np.where(np.any(self.stoich_matrix() > 0, axis=1))[0] ] # Remove duplicate names in subunits and complexes self.subunit_names = set(subunits) self.complex_names = set(complexes) # Create sparse matrix for monomer to complex stoichiometry i, j, v, shape = self._buildStoichMatrixMonomers() self._stoichMatrixMonomersI = i self._stoichMatrixMonomersJ = j self._stoichMatrixMonomersV = v self._stoichMatrixMonomersShape = shape # Mass balance matrix # All reaction mass balances should balance out to numerical zero balanceMatrix = self.stoich_matrix() * self.mass_matrix() massBalanceArray = np.sum(balanceMatrix, axis=0) assert ( np.max(np.absolute(massBalanceArray)) < 1e-8 ) # had to bump this up to 1e-8 because of flagella supercomplex stoichMatrix = self.stoich_matrix().astype(np.int64, order="F") self.prebuilt_matrices = mccBuildMatrices(stoichMatrix) # Add boolean array to mark reactions with unknown stoichiometries self.reaction_stoichiometry_unknown = np.array( self.reaction_stoichiometry_unknown )
[docs] def stoich_matrix(self): """ Builds a stoichiometric matrix based on each given complexation reaction. One reaction corresponds to one column in the stoichiometric matrix. """ shape = (self._stoich_matrix_I.max() + 1, self._stoich_matrix_J.max() + 1) out = np.zeros(shape, np.float64) out[self._stoich_matrix_I, self._stoich_matrix_J] = self._stoich_matrix_V return out
[docs] def mass_matrix(self): """ Builds a matrix with the same shape as the stoichiometric matrix, but with molecular masses as elements instead of stoichiometric constants """ shape = (self._stoich_matrix_I.max() + 1, self._stoich_matrix_J.max() + 1) out = np.zeros(shape, np.float64) out[self._stoich_matrix_I, self._stoich_matrix_J] = self._stoich_matrix_mass return out
[docs] def stoich_matrix_monomers(self): """ Returns the dense stoichiometric matrix for monomers from each complex """ out = np.zeros(self._stoichMatrixMonomersShape, np.float64) out[self._stoichMatrixMonomersI, self._stoichMatrixMonomersJ] = ( self._stoichMatrixMonomersV ) return out
# TODO: redesign this so it doesn't need to create a stoich matrix
[docs] def get_monomers(self, cplxId): """ Returns subunits for a complex (or any ID passed). If the ID passed is already a monomer returns the monomer ID again with a stoichiometric coefficient of one. """ info = self._moleculeRecursiveSearch( cplxId, self.stoich_matrix(), self.molecule_names ) subunits = [] subunit_stoich = [] for subunit, stoich in sorted(info.items()): subunits.append(subunit) subunit_stoich.append(stoich) return { "subunitIds": np.array(subunits), "subunitStoich": np.array(subunit_stoich), }
[docs] def _buildStoichMatrixMonomers(self): """ Builds a stoichiometric matrix where each column is a reaction that forms a complex directly from its constituent monomers. Since some reactions from the raw data are complexation reactions of complexes, this is different from the stoichiometric matrix generated by stoichMatrix(). """ stoichMatrixMonomersI = [] stoichMatrixMonomersJ = [] stoichMatrixMonomersV = [] for colIdx, id_complex in enumerate(self.ids_complexes): D = self.get_monomers(id_complex) rowIdx = self.molecule_names.index(id_complex) stoichMatrixMonomersI.append(rowIdx) stoichMatrixMonomersJ.append(colIdx) stoichMatrixMonomersV.append(1.0) for subunitId, subunitStoich in zip(D["subunitIds"], D["subunitStoich"]): rowIdx = self.molecule_names.index(subunitId) stoichMatrixMonomersI.append(rowIdx) stoichMatrixMonomersJ.append(colIdx) stoichMatrixMonomersV.append(-1.0 * subunitStoich) stoichMatrixMonomersI = np.array(stoichMatrixMonomersI) stoichMatrixMonomersJ = np.array(stoichMatrixMonomersJ) stoichMatrixMonomersV = np.array(stoichMatrixMonomersV) shape = (stoichMatrixMonomersI.max() + 1, stoichMatrixMonomersJ.max() + 1) return ( stoichMatrixMonomersI, stoichMatrixMonomersJ, stoichMatrixMonomersV, shape, )
[docs] def _findRow(self, product, speciesList): try: row = speciesList.index(product) except ValueError: row = -1 # Flag if not found so not a complex return row
[docs] def _findColumn(self, stoichMatrixRow): for i in range(0, len(stoichMatrixRow)): if int(stoichMatrixRow[i]) == 1: return i return -1 # Flag for monomer
[docs] def _moleculeRecursiveSearch(self, product, stoichMatrix, speciesList): row = self._findRow(product, speciesList) if row == -1: return {product: 1.0} col = self._findColumn(stoichMatrix[row, :]) if col == -1: return {product: 1.0} total = {} for i in range(0, len(speciesList)): if i == row: continue val = stoichMatrix[i][col] sp = speciesList[i] if val != 0: x = self._moleculeRecursiveSearch(sp, stoichMatrix, speciesList) for j in x: if j in total: total[j] += x[j] * (np.absolute(val)) else: total[j] = x[j] * (np.absolute(val)) return total