New metabolic pathway #1361
New metabolic pathway #1361
Conversation
There was a problem hiding this comment.
I think "external pathway" would be a more appropriate name for this file/class to make it clear that we are treating this pathway separately from the existing metabolic network.
| """ | ||
| Metabolism for new(external) pathway | ||
| Metabolism for new pathway sub-model | ||
|
|
||
| """ |
There was a problem hiding this comment.
I think adding good explanation on why this process needs to exist (why we are treating added pathways to be separate from the existing metabolic network) would be good to have here.
| Metabolism for new pathway sub-model | ||
|
|
||
| """ | ||
| from __future__ import absolute_import, division, print_function |
There was a problem hiding this comment.
This line is no longer needed after completely transitioning to Python 3.
|
|
||
| import wholecell.processes.process | ||
| from wholecell.utils import units | ||
| from wholecell.utils.constants import REQUEST_PRIORITY_METABOLISM_NEW |
There was a problem hiding this comment.
I think it would be okay to just use the existing REQUEST_PRIORITY_METABOLISM value here, since they have the same value anyway. Do you think there will be cases when this value will be different?
| from wholecell.utils.constants import REQUEST_PRIORITY_METABOLISM_NEW | ||
|
|
||
|
|
||
| class MetabolismNewPathway(wholecell.processes.process.Process): |
There was a problem hiding this comment.
Same comment as the filename. "New" pathway seems a little bit too vague.
| "NG002" "vioB" ["vioB"] 1257 4253 "+" ["NG002_RNA"] | ||
| "NG003" "vioC" ["vioC"] 4254 5543 "+" ["NG003_RNA"] | ||
| "NG004" "vioD" ["vioD"] 5544 6665 "+" ["NG004_RNA"] | ||
| "NG005" "vioE" ["vioE"] 6666 7241 "+" ["NG005_RNA"] |
There was a problem hiding this comment.
Same comment about the newline character here, and many of the new vioAE files.
| @@ -0,0 +1,2 @@ | |||
| "id" "stoichiometry" "direction" "catalyzed_by" "forward_rate" "reverse_rate" "kcat" "km" | |||
There was a problem hiding this comment.
What are the purposes of the forward and reverse rates columns?
| # Join datasets | ||
| for row in added_data: | ||
| data.append(row) | ||
| if added_data: |
There was a problem hiding this comment.
I don't think this change would be necessary if we remove the empty file requirements for the metabolism files?
runscripts/fireworks/fw_queue.py
Outdated
There was a problem hiding this comment.
This change should not be included in this PR.
wholecell/fireworks/nohup.out
Outdated
There was a problem hiding this comment.
All fireworks files should not be part of this PR.
|
|
||
| class MetabolismNewPathway(object): | ||
| def __init__(self, raw_data, sim_data): | ||
| # Build the abstractions needed for the metabolism of violacein |
There was a problem hiding this comment.
This comment needs to be more general.
| ratesFwd = [] # rate of reaction | ||
| ratesRev = [] |
There was a problem hiding this comment.
Are the forward and reverse rates used anywhere? If these are just remnants of the 2CS process these should be removed.
| stoichMatrixV = [] # Stoichometric coefficients | ||
|
|
||
| stoichMatrixMass = [] # molecular mass of molecules in stoichMatrixI | ||
| miscrnas_with_singleton_tus = sim_data.getter.get_miscrnas_with_singleton_tus() |
There was a problem hiding this comment.
This probably does not happen for metabolic reactions.
| compartment_ids_to_abbreviations = { | ||
| comp['id']: comp['abbrev'] for comp in raw_data.compartments | ||
| } | ||
| try: |
There was a problem hiding this comment.
I'm not sure how I feel about this big try-except block. As suggested earlier, I think it would be better to set an explicit boolean flag depending on whether external metabolic pathways exist that can be referenced from the simulation itself.
| if new_genes_folder.metabolic_reactions_new: | ||
| for reactionIndex, reaction in enumerate(new_genes_folder.metabolic_reactions_new): | ||
| reactionName = reaction["id"] | ||
| print(reaction) |
There was a problem hiding this comment.
All print statements should be removed from the final version.
| enzymeC = [] | ||
| for e in enzymeNames: | ||
| if e is None: | ||
| enzymeC.append(float(1)) |
There was a problem hiding this comment.
Why 1? Shouldn't this be set to infinity?
| if e is None: | ||
| enzymeC.append(float(1)) | ||
| else: | ||
| enzymeC.append(float(enzymeDict[e+'[c]'])) |
There was a problem hiding this comment.
Again, I think having an enzyme-reaction mapping matrix will be helpful here as well. You can simply get the vector of enzyme counts and multiply it with this matrix to the left to get the number of enzymes for each reaction.
| # Molecules needed are the change but only for the molecules corresponding to the negative terms in the stoich matrix | ||
| moleculesNeeded = abs(np.sum(self.stoich_matrix().clip(max=0), axis=1) * allMoleculesChanges) |
There was a problem hiding this comment.
I'm not sure if I understand the calculation here.. is there a reason we need to deviate from the 2CS method of calculating the molecules needed by simply doing this?
np.negative(allMoleculesChanges).clip(min=0)
| rates = [] | ||
| for colIdx in range(S.shape[1]): | ||
| #here we assume that the substrate is the first element in the reaction | ||
| flux = self.kcats[colIdx] * y[0]/ (self.kms + y[0]) |
There was a problem hiding this comment.
Should self.kms be self.kms[colIdx]?
There was a problem hiding this comment.
Also, is y[0] really the first element in the given reaction? Isn't this just the first metabolite in the stoichiometric matrix?
|
|
||
| rates = [] | ||
| for colIdx in range(S.shape[1]): | ||
| #here we assume that the substrate is the first element in the reaction |
There was a problem hiding this comment.
This is a pretty strong assumption that I think needs to be clearly highlighted elsewhere (maybe in the new reactions flat file as well). Ideally we should be able to clearly specify which metabolite is the substrate that goes into the MM equation.
ioanagherman
force-pushed
the
new-metabolic-pathway
branch
from
9a2e7b0 to
728d183
Compare
There was a problem hiding this comment.
Hi Ioana,
Thank you for finishing up the PR! I reviewed this together with @rjuenemann, and we had a few comments below, mostly related to code style. We also had a few bigger suggestions:
- It would be nice to have an analysis script that goes along with this new function to test and verify that the new process works as it should. The analysis script should test that the new genes are being expressed, the new metabolic reactions have flux, and the fluxes are responding to enzyme counts. I believe you are probably already testing and confirming this using your own pipeline, but it would be better to have something in our codebase to test this new process in future updates. We can also discuss setting up an automated daily tests for this new optional features, but I think an analysis script would be sufficient for now.
- We have a markdown file at
docs/misc/add_new_gene.mdwhere we document how we can add new genes. It would be great if you could add a section here describing the process of adding new metabolic genes and reactions. - We use PR descriptions as a way to provide and record detailed information about the changes to the codebase introduced by the PR. It would be great if you could update the descriptions to provide more information about this particular update, as well as including the output of the new analysis script.
.gitignore
Outdated
| wholecell/fireworks/my_launchpad.yaml | ||
| wholecell/fireworks/my_qadapter.yaml |
There was a problem hiding this comment.
These two files should be ignored if you place them at the top level of the wcEcoli directory.
.gitignore
Outdated
| wholecell/fireworks/nohup.out | ||
| runscripts/fireworks/vioAE.sh |
There was a problem hiding this comment.
We try not to have custom files that are specific to a particular branch on this list. If you want to put in the vioAE.sh file somewhere, it would be better to put it in a new directory under the runscripts directory, similar to how we have done for runscripts/paper.
There was a problem hiding this comment.
I think this file should be renamed to metabolism_external_pathway.py such that the class name matches the file name, as is the case for all other files in the same directory.
| @@ -241,6 +241,11 @@ def __init__(self, operons_on: bool, remove_rrna_operons: bool, remove_rrff: boo | |||
| new_gene_shared_files = ['genes', 'rnas', 'proteins', | |||
| 'rna_half_lives', | |||
| 'protein_half_lives_measured'] | |||
|
|
|||
| if os.path.isfile(os.path.join(FLAT_DIR+'/'+new_gene_path, 'metabolic_reactions_new.tsv')): | |||
There was a problem hiding this comment.
You can do os.path.join(FLAT_DIR, new_gene_path, 'metabolic_reactions_new.tsv') to remove the OS-dependent slash character.
| @@ -241,6 +241,11 @@ def __init__(self, operons_on: bool, remove_rrna_operons: bool, remove_rrff: boo | |||
| new_gene_shared_files = ['genes', 'rnas', 'proteins', | |||
| 'rna_half_lives', | |||
| 'protein_half_lives_measured'] | |||
|
|
|||
| if os.path.isfile(os.path.join(FLAT_DIR+'/'+new_gene_path, 'metabolic_reactions_new.tsv')): | |||
| self.list_of_dict_filenames.append(os.path.join(new_gene_path, 'metabolic_reactions_new.tsv')) | |||
There was a problem hiding this comment.
I think it would be better to name this file metabolic_reactions_external.tsv to better match the naming for the other components.
| self.molecule_names = sim_data.process.metabolism_external_pathway.molecule_names | ||
| self.molecules = self.bulkMoleculesView(self.molecule_names) | ||
|
|
||
| self.enzyme_names = [item for sublist in sim_data.process.metabolism_external_pathway.enzymes.values() for |
There was a problem hiding this comment.
This variable is not used anywhere outside of this method, so it would make sense to turn this into a local variable (enzyme_names, not self.enzyme_names). Also list(chain.from_iterable(enzymes.values())) might be a more readable way of doing this concatenation.
There was a problem hiding this comment.
You will need to import chain through from itertools import chain.
There was a problem hiding this comment.
A simpler trick to concatenate all the lists in a list, e.g. enzymes[0] + enzymes[1] + enzymes[2]..., is sum(enzymes, []). This begins with [] then "adds" all the elements.
|
|
||
|
|
||
| def molecules_to_next_time_step(self, moleculeDict, enzymeDict, cellVolume, | ||
| nAvogadro, timeStepSec, method="LSODA", jit=True): |
There was a problem hiding this comment.
jit is not being used anywhere inside this function, so this could be removed.
| for sub in self.substrates: | ||
| subCounts.append(moleculeDict[sub[0]]) | ||
| subNames.append(sub[0]) |
There was a problem hiding this comment.
Are we only looking at the counts of the first substrate for each reaction? If so, we should make it clear in a comment that this is an assumption we are currently making for enzyme kinetics.
| if self.enzymes[reaction]: | ||
| enzymeCounts.append(enzymeDict[self.enzymes[reaction][0]]) | ||
| enzymeNames.append(self.enzymes[reaction][0]) |
There was a problem hiding this comment.
Same comment as above -- there should be a comment on why we're only picking out the first enzyme from a list.
| #this might need to be modelled differently | ||
| else: | ||
| enzymeCounts.append(1) | ||
| enzymeNames.append('None') |
There was a problem hiding this comment.
It's unclear what's being assumed here. If a reaction does not have an enzyme, are we assuming that the reaction is spontaneous? Would setting the enzyme count to 1 ensure that this reaction does not bottleneck the pathway? I guess for vioAE we have enzymes for each reaction so this does not matter, but it would be great if the comment was a bit more specific on what will need to be done for this block if we have reactions without enzymes.
There was a problem hiding this comment.
Also, using the string 'None' as a filler value for the list of enzyme names is a little bit risky. It looks like this list is not being used later in the function -- was this made just for debugging purposes?
There was a problem hiding this comment.
So actually the last reaction in the violacein pathway is spontaneous according to Metacyc https://biocyc.org/META/NEW-IMAGE?type=PATHWAY&object=PWY-7040&detail-level=2
This is something that I wanted to discuss with you actually because I am not sure how to model this situation. Maybe we can have a quick call or discuss via email.
There was a problem hiding this comment.
+1 to what @ggsun said - let us know if you have any questions!
|
Hi everyone, I updated the code and addressed your comments. Could you please have another look @ggsun @rjuenemann? Thank you! |
There was a problem hiding this comment.
Hi Ioana,
Thank you for addressing all of our comments and adding those analysis scripts! We have left specific comments on the code in addition to these more general comments:
- It would be great if you could go through our code style guide and make some code style changes before merging this branch. Some things that stood out to us during review were some camelCase variable names that should be changed to underscore, some lines that were longer than the line length limits, and typos in variable names and comments.
- It would be great if you could provide in this PR some sample outputs from the analysis plots that you added, even if it's for a much simpler set of simulations. This will help us see what the plot should look like, and also help us review the code.
- In the rest of our model, we use the units of mmol/g DCW/hr to plot flux values. It might be best to make these plots consistent with this system of units as well, for easier comparison with existing flux values for internal metabolic reactions. One example is in the file
models/ecoli/analysis/comparison/proteomics_fluxomics_validation.py.
All of these comments were coauthored with @rjuenemann.
| cell_paths = self.ap.get_cells() | ||
| sim_dir = cell_paths[0] | ||
| plt.figure() | ||
| #for sim_dir in cell_paths: |
| for m in range(fluxes_external_pathway.shape[1]): | ||
| plt.plot(time / 60., fluxes_external_pathway[:, m], | ||
| label='Flux ' + str(m)) |
There was a problem hiding this comment.
If there's an ID associated with each of these reactions, it might be more useful to label each of the lines with their IDs than their indexes.
| sim_dir = cell_paths[0] | ||
| plt.figure() | ||
| #for sim_dir in cell_paths: | ||
| simOutDir = os.path.join(sim_dir, 'simOut') |
There was a problem hiding this comment.
simOutDir is not being used anywhere, so this line should be removed. This would also let you remove line 29.
| bulk_molecules = read_stacked_columns(cell_paths, 'BulkMolecules', 'counts') | ||
|
|
||
| with open(os.path.join(simOutDir, 'BulkMolecules', 'attributes.json')) as f: | ||
| data = json.load(f, object_pairs_hook=OrderedDict) | ||
| data_b = data['objectNames'] | ||
| # Molecule counts | ||
| plt.figure() | ||
| plt.xlabel("Time (min)") | ||
| plt.title("Molecule counts") | ||
|
|
||
| for m in range(len(molecule_names)): | ||
| plt.subplot(int(math.sqrt(len(molecule_names)))+1, int(math.sqrt(len(molecule_names)))+1, m+1) | ||
| molecule_couts = bulk_molecules[:, data_b.index(molecule_names[m])] | ||
| plt.plot(time / 60., molecule_couts) | ||
| plt.ylabel('\n'.join(molecule_names[m].split('-'))) |
There was a problem hiding this comment.
We have a helper function called read_stacked_bulk_molecules (this should actually be imported in the template file so you don't even need to add an import statement) that can be used to specify which molecules counts to read in from the bulk molecules table. Using this helper function will both simplify the code and save memory.
| bulk_molecules = read_stacked_columns(cell_paths, 'BulkMolecules', 'counts') | ||
|
|
||
| with open(os.path.join(simOutDir, 'BulkMolecules', 'attributes.json')) as f: | ||
| data = json.load(f, object_pairs_hook=OrderedDict) | ||
| data_b = data['objectNames'] | ||
| # Molecule counts | ||
| plt.figure() | ||
| plt.xlabel("Time (min)") | ||
| plt.title("Molecule counts") | ||
|
|
||
| for m in range(len(molecule_names)): | ||
| plt.subplot(int(math.sqrt(len(molecule_names)))+1, int(math.sqrt(len(molecule_names)))+1, m+1) | ||
| molecule_couts = bulk_molecules[:, data_b.index(molecule_names[m])] | ||
| plt.plot(time / 60., molecule_couts) | ||
| plt.ylabel('\n'.join(molecule_names[m].split('-'))) |
There was a problem hiding this comment.
I'm not entirely sure what's going on with the string manipulation here. Is there a specific use for the dash character in the IDs of these molecules?
There was a problem hiding this comment.
The names were too long (CPD-, OXYGEN-MOLECULE etc) so the labels were overlapping.
|
|
||
| * id (e.g. ‘rxnpath1’), | ||
|
|
||
| * Stoichiometry represented as a dictionary where the keys are the reactants and the products, and the values represent the coefficients for each reactant (positive coefficient) or product (negative coefficient) |
There was a problem hiding this comment.
Should reactants have negative coefficients, and products positive?
|
|
||
| * metabolites.tsv - list of the metabolites produced by the reaction. Specify for each of them: | ||
|
|
||
| * Their ID as shown on Ecocyc (e.g. CPD-11890) |
There was a problem hiding this comment.
The second 'c' should be capitalized in EcoCyc.
|
|
||
| * Molecular charge. | ||
|
|
||
| * Smiles. |
There was a problem hiding this comment.
SMILES should be in all caps.
|
|
||
| A new process was created to represent the metabolism corresponding to these reactions. This new process is based on Michaelis Menten kinetics and it is run at each time step before the Metabolism process. To achieve this, two Python files were added: | ||
|
|
||
| * reconstruction/ecoli/dataclasses/process/metabolism_external_pathway.py an – This script loads the raw data from the files populated above and stores it in the sim_data object. It also makes it accessible as an instance of the new process class. In this file we also define the methods that represent the Michaelis Menten model. This script is called when runParca is run. |
There was a problem hiding this comment.
Filenames should be wrapped inside backticks (`) to have them formatted in a fixed-width font.
|
|
||
| * If these files are not empty, the WCM will run the metabolic_reactions_external process. | ||
|
|
||
| A new process was created to represent the metabolism corresponding to these reactions. This new process is based on Michaelis Menten kinetics and it is run at each time step before the Metabolism process. To achieve this, two Python files were added: |
There was a problem hiding this comment.
It might be worth highlighting here the exact formulation of the Michaelis-Menten kinetics as part of the documentation.
…genes (#1442) * Add analysis script to generate tables of subgenerationally expressed genes * Add IGNORE_FIRST_N_GENS option to subgenerational expression table plot * Skip analysis if not enough generations are run
… structures (#1443) * Skip rRNA operon structures analysis if using nonstandard rRNA operon structures * Fix numpy deprecation warning in 2CS
* Add docs/maintaining_the_sherlock_environment.md. * Cross-links and minor updates in other docs.
* update to ecocyc 28.1 * remove new hybrid TUs * handle new compartments and too short transcription units * remove debugging statement * updated ecocyc files * update ecocyc files * unremove corrected metabolic reactions * remove corrected compartment
…1446) * initial implementation of activating fewer ribosomes (if needed) * initial implementation for inter ribosome distance analysis * change estimation of number of RNA init events per copy * change estimation of number of ribosomes per cistron * add more print statements for analyzing impact of activating fewer ribos * move setting protein init prob for overcrowded proteins inside rescaling check * add ppgpp_concentration single analysis plot * add comments to single ppgpp conc plot * remove debugging statements from process that activates fewer ribosomes * add listener for whether n active ribos had to be reduced to prevent crowding * change mRNA_is_overcrowded listener to bool data type * change tu_is_overcrowded listener initialization to bool * add marker for monomers with close ribosomes but low cistron counts
* add listeners for max_p in transcript and polypeptide initiation * correct typo in ribosomal subunit label in a comment
* Switch from the deprecated Google Container Registry to Artifact Registry "wcm". * Update the GitHub xianyi repo URL to OpenMathLib. * Update to OpenBLAS v0.3.27 . It has a bunch of bug fixes and support for Apple Silicon. TODO items for the WCM runtime Docker container: * Compile OpenBLAS by default? * Upgrade from Python 3.11.3 to 3.11.8? * Use a Python virtual env in the container? `WARNING: Running pip as the 'root' user can result in broken permissions and conflicting behaviour with the system package manager. It is recommended to use a virtual environment instead: https://pip.pypa.io/warnings/venv` * Fix runtime deprecation warnings about numpy.distutils, msvccompiler.py, and Blas libraries.
* add in placeholders to handle special complexes separately * include active RNAP and ribosomes in reported complex counts for EcoCyc * full ribosome complex id does not have a compartment * import ecocyc 28.1 preview files * correct peroxisome compartment
…justments based on N-end rule for ADENYLATECYC-MONOMER and SPOT-MONOMER, which now have updated sequences. (#1455)
* Manually update pip, then let requirements.txt install the other installers so it can require `setuptools==73.0.1` which is the last version that's compatible with Aesara. It's OK to use some older versions (the pyenv on Sherlock uses 69.0.2) but not `> 73.0.1`. * macOS: Install `cmake` via brew. * Sherlock: `module load git/2.45.1` (git/2.27.0 is no longer available). * Add comments about the different ways to run `test_blas.py` and more consistent results from openblas/0.3.28 and macOS Accelerate. Newer numpy and scipy releases can use Accelerate by default and maybe openblas/0.3.28, but they're not compatible with Aesara. **Note:** Aesara is no longer actively maintained. It's not compatible with newer versions of numpy, scipy, and numba. At some point we'll need to replace it. Its only use in wcEcoli is its Automatic Differentiation feature for `rna_decay.py`. It's likely that PyTorch or TensorFlow can do that, and will remain well maintained. **Note:** A pytest warning says `It is recommended to use setuptools < 60.0 for those Python versions` (< 3.12) but cvxpy 1.3.2 requires `setuptools > 65.5.1`.
conflicts with remote branch resolved
Reverting accidental push to upstream.
Revert accidental push to upstream.
Reverting accidental push to upstream.
The changes made on this branch correspond to the addition of an external metabolic pathway, where reactions are catalyzed by exogenous genes. The main changes include:
This process is modelled using Michaelis Menten kinetics and it uses the reactions and parameters specified in the metabolic_reactions_external.tsv file.
a.models/ecoli/analysis/single/flux_external_metabolic_pathway.py - creates a plot of the fluxes corresponding
to the external reactions added. This will be one plot for each cell.
b. models/ecoli/analysis/multigen/flux_external_metabolic_pathway.py - creates a plot of the fluxes corresponding
to the external reactions added. This will show all these fluxes for all generations simulated in one variant.
c. models/ecoli/analysis/single/molecules_external_pathway.py - creates plots of molecule counts for each
molecule involved in the external reactions, for each cell separately, from every generation.
d. models/ecoli/analysis/multigen/molecules_external_pathway.py - creates plots of molecule counts for each
molecule involved in the external reactions, for all cells from all generations corresponding to a variant.