Update major raw data flat files to latest versions pulled from EcoCyc

models/ecoli/analysis/causality_network/build_network.py

@@ -104,15 +104,15 @@
 
 # Equilibrium complexes that are formed from deactivated equilibrium reactions,
 # but are reactants in a complexation reaction
-EQUILIBRIUM_COMPLEXES_IN_COMPLEXATION = ["CPLX0-7620[c]", "CPLX0-7701[c]",
-	"CPLX0-7677[c]", "MONOMER0-1781[c]", "CPLX0-7702[c]"]
+EQUILIBRIUM_COMPLEXES_IN_COMPLEXATION = ["CPLX0-7701[c]", "CPLX0-7677[c]",
+	"MONOMER0-1781[c]", "CPLX0-7702[c]"]
 
 # Metabolites that are used as ligands in equilibrium, but do not participate
 # in any metabolic reactions
-METABOLITES_ONLY_IN_EQUILIBRIUM = ["4FE-4S[c]", "NITRATE[p]"]
+METABOLITES_ONLY_IN_EQUILIBRIUM = ["4FE-4S[c]"]
 
 # Molecules in 2CS (two component system) reactions that are not proteins
-NONPROTEIN_MOLECULES_IN_2CS = ["ATP[c]", "ADP[c]", "WATER[c]", "PI[c]",
+NONPROTEIN_MOLECULES_IN_2CS = ["ATP[c]", "ADP[c]", "WATER[c]", "Pi[c]",
 	"PROTON[c]", "PHOSPHO-PHOB[c]"]
 
 COMPARTMENTS = {
@@ -197,6 +197,7 @@ def _build_network(self):
 		self._add_metabolism_and_metabolites()
 		self._add_equilibrium()
 		self._add_regulation()
+		self._remove_hanging_edges()
 
 		# Check for network sanity (optional)
 		if self.check_sanity:
@@ -990,6 +991,25 @@ def _find_duplicate_nodes(self):
 				% (len(duplicate_ids), duplicate_ids))
 
 
+	def _remove_hanging_edges(self):
+		# type: () -> None
+		"""
+		Remove edges that are not connected to existing nodes.
+		"""
+		# Get set of all node IDs
+		node_ids = {node.node_id for node in self.node_list}
+		disconnected_edge_indexes = []
+
+		# Find edges whose source and destination nodes are not defined
+		for index, edge in enumerate(self.edge_list):
+			if edge.src_id not in node_ids or edge.dst_id not in node_ids:
+				disconnected_edge_indexes.append(index)
+
+		# Remove these edges
+		for index in disconnected_edge_indexes[::-1]:
+			self.edge_list.pop(index)
+
+
 	def _append_edge(self, type_, src, dst, stoichiometry=""):
 		# type: (str, str, str, Union[None, str, int]) -> None
 		"""

models/ecoli/analysis/cohort/expression_dynamics.py

@@ -170,8 +170,8 @@ def do_plot(self, seedOutDir, plotOutDir, plotOutFileName, simDataFile, validati
 			counts = (
 				proteinMonomerCounts[:, protein_idx],
 				proteinMonomerCounts[:, protein_idx_burst],
-				rnaMonomerCounts[:, sim_data.relation.rna_index_to_monomer_mapping][:, protein_idx],
-				rnaMonomerCounts[:, sim_data.relation.rna_index_to_monomer_mapping][:, protein_idx_burst]
+				rnaMonomerCounts[:, sim_data.relation.RNA_to_monomer_mapping][:, protein_idx],
+				rnaMonomerCounts[:, sim_data.relation.RNA_to_monomer_mapping][:, protein_idx_burst]
 				)
 			line_color = (
 				EXP_COLOR,

models/ecoli/analysis/cohort/proteinFoldChangeVsTranscriptionFrequency.py

@@ -122,7 +122,7 @@ def do_plot(self, seedOutDir, plotOutDir, plotOutFileName, simDataFile, validati
 			initiationEventsPerRna = rnapData.readColumn("rnaInitEvent").sum(axis = 0)
 
 			# Map transcription initiation events to monomers
-			initiationEventsPerMonomer = initiationEventsPerRna[sim_data.relation.rna_index_to_monomer_mapping]
+			initiationEventsPerMonomer = initiationEventsPerRna[sim_data.relation.RNA_to_monomer_mapping]
 
 			# Load cell mass
 			cellMassInitial = TableReader(os.path.join(simOutDir, "Mass")).readColumn("cellMass")[0]

models/ecoli/analysis/multigen/functionalUnits.py

@@ -63,7 +63,7 @@ def do_plot(self, seedOutDir, plotOutDir, plotOutFileName, simDataFile, validati
 
 		# Load sim data
 		sim_data = cPickle.load(open(simDataFile, "rb"))
-		rnaIds = sim_data.process.transcription.rna_data["id"][sim_data.relation.rna_index_to_monomer_mapping] # orders rna IDs to match monomer IDs
+		rnaIds = sim_data.process.transcription.rna_data["id"][sim_data.relation.RNA_to_monomer_mapping] # orders rna IDs to match monomer IDs
 
 		# Make views for monomers
 		ids_complexation = sim_data.process.complexation.molecule_names

models/ecoli/analysis/multigen/functionalUnitsFC.py

@@ -51,7 +51,7 @@ def do_plot(self, seedOutDir, plotOutDir, plotOutFileName, simDataFile, validati
 
 		# Load sim data
 		sim_data = cPickle.load(open(simDataFile, "rb"))
-		rnaIds = sim_data.process.transcription.rna_data["id"][sim_data.relation.rna_index_to_monomer_mapping] # orders rna IDs to match monomer IDs
+		rnaIds = sim_data.process.transcription.rna_data["id"][sim_data.relation.RNA_to_monomer_mapping] # orders rna IDs to match monomer IDs
 		rnaIds = rnaIds.tolist()
 		ids_complexation = sim_data.process.complexation.molecule_names # Complexe of proteins, and protein monomers
 		ids_complexation_complexes = sim_data.process.complexation.ids_complexes # Only complexes

models/ecoli/analysis/multigen/probProteinExistAndDouble.py

@@ -54,7 +54,7 @@ def do_plot(self, seedOutDir, plotOutDir, plotOutFileName, simDataFile, validati
 			initiationEventsPerRna = rnapData.readColumn("rnaInitEvent").sum(axis = 0)
 
 			# Map transcription initiation events to monomers
-			initiationEventsPerMonomer = initiationEventsPerRna[sim_data.relation.rna_index_to_monomer_mapping]
+			initiationEventsPerMonomer = initiationEventsPerRna[sim_data.relation.RNA_to_monomer_mapping]
 
 			# Log data
 			monomerExistMultigen[gen_idx,:] = monomerExist

models/ecoli/analysis/multigen/proteinAvgCountVsBurstSize.py

@@ -53,7 +53,7 @@ def do_plot(self, seedOutDir, plotOutDir, plotOutFileName, simDataFile, validati
 			initiationEventsPerRna = rnapData.readColumn("rnaInitEvent").sum(axis = 0)
 
 			# Map transcription initiation events to monomers
-			initiationEventsPerMonomer = initiationEventsPerRna[sim_data.relation.rna_index_to_monomer_mapping]
+			initiationEventsPerMonomer = initiationEventsPerRna[sim_data.relation.RNA_to_monomer_mapping]
 
 			# Log data
 			monomerExistMultigen[gen_idx,:] = monomerExist

models/ecoli/analysis/multigen/proteinCountVsFoldChange.py

@@ -53,7 +53,7 @@ def do_plot(self, seedOutDir, plotOutDir, plotOutFileName, simDataFile, validati
 			initiationEventsPerRna = rnapData.readColumn("rnaInitEvent").sum(axis = 0)
 
 			# Map transcription initiation events to monomers
-			initiationEventsPerMonomer = initiationEventsPerRna[sim_data.relation.rna_index_to_monomer_mapping]
+			initiationEventsPerMonomer = initiationEventsPerRna[sim_data.relation.RNA_to_monomer_mapping]
 
 			# Log data
 			monomerExistMultigen[gen_idx,:] = monomerExist

models/ecoli/analysis/multigen/proteinCountsValidation.py

@@ -11,6 +11,7 @@
 from six.moves import cPickle
 
 from wholecell.io.tablereader import TableReader
+from wholecell.utils.protein_counts import get_simulated_validation_counts
 
 from models.ecoli.analysis.AnalysisPaths import AnalysisPaths
 from wholecell.analysis.analysis_tools import exportFigure
@@ -22,8 +23,8 @@ def do_plot(self, seedOutDir, plotOutDir, plotOutFileName, simDataFile, validati
 		sim_data = cPickle.load(open(simDataFile, "rb"))
 		validation_data = cPickle.load(open(validationDataFile, "rb"))
 
-		ids_translation = sim_data.process.translation.monomer_data["id"].tolist()
-		schmidt_idx = [ids_translation.index(x) for x in validation_data.protein.schmidt2015Data["monomerId"].tolist()]
+		monomer_ids = sim_data.process.translation.monomer_data["id"]
+		schmidt_ids = validation_data.protein.schmidt2015Data["monomerId"]
 
 		schmidt_counts = validation_data.protein.schmidt2015Data["glucoseCounts"]
 
@@ -41,21 +42,23 @@ def do_plot(self, seedOutDir, plotOutDir, plotOutFileName, simDataFile, validati
 
 			simOutDir = os.path.join(simDir, "simOut")
 
-			monomerCounts = TableReader(os.path.join(simOutDir, "MonomerCounts"))
-			avgCounts = monomerCounts.readColumn("monomerCounts").mean(axis=0)
-			sim_schmidt_counts = avgCounts[schmidt_idx]
+			monomer_counts_reader = TableReader(os.path.join(simOutDir, "MonomerCounts"))
+			monomer_counts = monomer_counts_reader.readColumn("monomerCounts")
+
+			sim_schmidt_counts, val_schmidt_counts = get_simulated_validation_counts(
+				schmidt_counts, monomer_counts, schmidt_ids, monomer_ids)
 
 			sim_schmidt_counts_multigen.append(sim_schmidt_counts)
 
 		sim_schmidt_counts_multigen = (np.array(sim_schmidt_counts_multigen)).mean(axis = 0)
 
 		axis = plt.subplot(1,1,1)
 
-		axis.plot(np.log10(schmidt_counts + 1), np.log10(sim_schmidt_counts_multigen + 1), 'o', color = "black", markersize = 6, alpha = 0.1, zorder = 1, markeredgewidth = 0.0)
+		axis.plot(np.log10(val_schmidt_counts + 1), np.log10(sim_schmidt_counts_multigen + 1), 'o', color = "black", markersize = 6, alpha = 0.1, zorder = 1, markeredgewidth = 0.0)
 		# print(pearsonr( np.log10(sim_schmidt_counts_mulitgen + 1), np.log10(schmidtCounts + 1) )[0])
 
 		maxLine = np.ceil(
-						max((np.log10(schmidt_counts + 1)).max(),
+						max((np.log10(val_schmidt_counts + 1)).max(),
 						(np.log10(sim_schmidt_counts_multigen + 1)).max())
 					)
 		plt.plot([0, maxLine], [0, maxLine], '-k')

models/ecoli/analysis/multigen/proteinExistVsBurstSize.py

@@ -51,7 +51,7 @@ def do_plot(self, seedOutDir, plotOutDir, plotOutFileName, simDataFile, validati
 			initiationEventsPerRna = rnapData.readColumn("rnaInitEvent").sum(axis = 0)
 
 			# Map transcription initiation events to monomers
-			initiationEventsPerMonomer = initiationEventsPerRna[sim_data.relation.rna_index_to_monomer_mapping]
+			initiationEventsPerMonomer = initiationEventsPerRna[sim_data.relation.RNA_to_monomer_mapping]
 
 			# Log data
 			monomerExistMultigen[gen_idx,:] = monomerExist

models/ecoli/analysis/multigen/proteinFoldChangeVsRnaDeg.py

@@ -60,7 +60,7 @@ def do_plot(self, seedOutDir, plotOutDir, plotOutFileName, simDataFile, validati
 
 		# uniqueBurstSizes = np.unique(initiationEventsPerMonomerMultigen)
 		degradationRates = sim_data.process.transcription.rna_data['deg_rate'].asNumber(1 / units.s)
-		degradationRatesByMonomer = degradationRates[sim_data.relation.rna_index_to_monomer_mapping]
+		degradationRatesByMonomer = degradationRates[sim_data.relation.RNA_to_monomer_mapping]
 		uniqueDegRate = np.unique(degradationRatesByMonomer)
 
 		# burstSizeToPlot = np.zeros(0)

models/ecoli/analysis/multigen/proteinFoldChangeVsTranscriptionEvents.py

@@ -47,7 +47,7 @@ def do_plot(self, seedOutDir, plotOutDir, plotOutFileName, simDataFile, validati
 			initiationEventsPerRna = rnapData.readColumn("rnaInitEvent").sum(axis = 0)
 
 			# Map transcription initiation events to monomers
-			initiationEventsPerMonomer = initiationEventsPerRna[sim_data.relation.rna_index_to_monomer_mapping]
+			initiationEventsPerMonomer = initiationEventsPerRna[sim_data.relation.RNA_to_monomer_mapping]
 
 			# Log data
 			ratioFinalToInitialCount[np.logical_and(ratioFinalToInitialCount == 1., proteinMonomerCounts[0,:] == 0)] = np.nan

models/ecoli/analysis/multigen/proteinFoldChangeVsTranslationEff.py

@@ -52,7 +52,7 @@ def do_plot(self, seedOutDir, plotOutDir, plotOutFileName, simDataFile, validati
 			initiationEventsPerRna = rnapData.readColumn("rnaInitEvent").sum(axis = 0)
 
 			# Map transcription initiation events to monomers
-			initiationEventsPerMonomer = initiationEventsPerRna[sim_data.relation.rna_index_to_monomer_mapping]
+			initiationEventsPerMonomer = initiationEventsPerRna[sim_data.relation.RNA_to_monomer_mapping]
 
 			# Log data
 			monomerExistMultigen[gen_idx,:] = monomerExist

models/ecoli/analysis/multigen/rnaVsProteinPerCell.py

@@ -54,7 +54,7 @@ def do_plot(self, seedOutDir, plotOutDir, plotOutFileName, simDataFile, validati
 
 		# Load sim data
 		sim_data = cPickle.load(open(simDataFile, "rb"))
-		rnaIds = sim_data.process.transcription.rna_data["id"][sim_data.relation.rna_index_to_monomer_mapping] # orders rna IDs to match monomer IDs
+		rnaIds = sim_data.process.transcription.rna_data["id"][sim_data.relation.RNA_to_monomer_mapping] # orders rna IDs to match monomer IDs
 
 		# Make views for monomers
 		ids_complexation = sim_data.process.complexation.molecule_names

models/ecoli/analysis/multigen/subgenerationalTranscription.py

@@ -241,8 +241,12 @@ def do_plot(self, seedOutDir, plotOutDir, plotOutFileName, simDataFile, validati
 		plotGreen = 0
 		plotBlue = 0
 
+		mRNA_id_set = set(mRnaIdsOrdered)
+
 		essential_genes_rna = validation_data.essential_genes.essential_RNAs
 		for g in essential_genes_rna:
+			if g not in mRNA_id_set:
+				continue
 			i = np.where(mRnaIdsOrdered == str(g))[0][0]
 			f = transcribedBoolOrdered[i]
 
@@ -294,8 +298,12 @@ def do_plot(self, seedOutDir, plotOutDir, plotOutFileName, simDataFile, validati
 		unknown = {"r": 0, "g": 0, "b": 0}
 		resistance = {"r": 0, "g": 0, "b": 0}
 		for frameID, function_ in six.viewitems(geneFunctions):
-			if function_ in ["Unknown function", "Unclear/under-characterized"]:
+			try:
 				i = np.where([frameID in x for x in mRnaIdsOrdered])[0][0]
+			# Skip over genes that aren't used in the simulation
+			except IndexError:
+				continue
+			if function_ in ["Unknown function", "Unclear/under-characterized"]:
 				f = transcribedBoolOrdered[i]
 				if f == 0.0:
 					unknown["r"] += 1
@@ -305,7 +313,6 @@ def do_plot(self, seedOutDir, plotOutDir, plotOutFileName, simDataFile, validati
 					unknown["g"] += 1
 
 			elif function_ in ["Antibiotic resistance", "Toxin/antitoxin"]:
-				i = np.where([frameID in x for x in mRnaIdsOrdered])[0][0]
 				f = transcribedBoolOrdered[i]
 				if f == 0.0:
 					resistance["r"] += 1

models/ecoli/analysis/single/cell_wall_expression.py

@@ -3,8 +3,7 @@
 Note: This is not an exhaustive list of all components of cell wall.
 """
 
-from __future__ import absolute_import, division, print_function
-
+import pickle
 import os
 
 from matplotlib import pyplot as plt
@@ -20,49 +19,58 @@
 CMAP_COLORS = [[shade/255. for shade in color] for color in CMAP_COLORS_255]
 CMAP_OVER = [0, 1, 0.75]
 
-cell_wall_proteins = {
+CELL_WALL_PROTEINS = {
     # murein
-    'CPD-12261[p]': 'murein',
-    'CPD-12231[p]':'peptidoglycan dimer',
-    'C6[p]':'lipid II',
-    'C5[c]': 'lipid I',
-    'CPD-9646[c]': 'undecaprenyl phosphate',
-    'UNDECAPRENYL-DIPHOSPHATE[c]': 'undecaprenyl diphosphate',
+    'CPD-12261': 'murein',
+    'CPD-12231':'peptidoglycan dimer',
+    'C6':'lipid II',
+    'C5': 'lipid I',
+    'CPD-9646': 'undecaprenyl phosphate',
+    'UNDECAPRENYL-DIPHOSPHATE': 'undecaprenyl diphosphate',
 
     # catalyzes lipid II transfer to outer membrane
-    'G6561-MONOMER[i]': 'lipid II flippase (MurJ)',
-    'EG10344-MONOMER[i]': 'essential cell division protein (FtsW)',
+    'G6561-MONOMER': 'lipid II flippase (MurJ)',
+    'EG10344-MONOMER': 'essential cell division protein (FtsW)',
 
     # catalyzes lipid II formation
-    'NACGLCTRANS-MONOMER[c]': 'MurG',
+    'NACGLCTRANS-MONOMER': 'MurG',
 
     # catalyzes lipid I formation
-    'PHOSNACMURPENTATRANS-MONOMER[i]': 'MraY',
+    'PHOSNACMURPENTATRANS-MONOMER': 'MraY',
 
     # PBP
-    'CPLX0-7717[i]': 'PBP1A', # transglycosylase-transpeptidase ~100
-    'CPLX0-3951[i]': 'PBP1B', # transglycosylase-transpeptidase ~100
-    'G7322-MONOMER[i]': 'PBP1C', # transglycosylase
-    'EG10606-MONOMER[i]': 'PBP2', # transpeptidase ~20
-    'EG10341-MONOMER[i]': 'PBP3', # transglycosylase-transpeptidase ~50
-    'EG10202-MONOMER[p]': 'PBP4', # DD-endopeptidase, DD-carboxypeptidase ~110
-    'EG10201-MONOMER[i]': 'PBP5', # DD-caroxypeptidase ~1,800
-    'EG10203-MONOMER[i]': 'PBP6', # DD-carbocypeptidase ~600
-    'EG12015-MONOMER[p]': 'PBP7', # DD-endopeptidase
+    'CPLX0-7717': 'PBP1A', # transglycosylase-transpeptidase ~100
+    'CPLX0-3951': 'PBP1B', # transglycosylase-transpeptidase ~100
+    'G7322-MONOMER': 'PBP1C', # transglycosylase
+    'EG10606-MONOMER': 'PBP2', # transpeptidase ~20
+    'EG10341-MONOMER': 'PBP3', # transglycosylase-transpeptidase ~50
+    'EG10202-MONOMER': 'PBP4', # DD-endopeptidase, DD-carboxypeptidase ~110
+    'EG10201-MONOMER': 'PBP5', # DD-caroxypeptidase ~1,800
+    'EG10203-MONOMER': 'PBP6', # DD-carbocypeptidase ~600
+    'EG12015-MONOMER': 'PBP7', # DD-endopeptidase
 
 }
 
 class Plot(singleAnalysisPlot.SingleAnalysisPlot):
 	def do_plot(self, simOutDir, plotOutDir, plotOutFileName,
 				simDataFile, validationDataFile, metadata):
+		with open(simDataFile, 'rb') as f:
+			sim_data = pickle.load(f)
+
 		# Listeners used
 		main_reader = TableReader(os.path.join(simOutDir, 'Main'))
 
 		# Load data
 		initial_time = main_reader.readAttribute('initialTime')
 		time = main_reader.readColumn('time') - initial_time
-		(counts,) = read_bulk_molecule_counts(simOutDir,
-											  (cell_wall_proteins.keys(),))
+
+		existing_molecule_ids = [
+			mol_id + sim_data.getter.get_compartment_tag(mol_id)
+			for mol_id in CELL_WALL_PROTEINS.keys()
+			if sim_data.getter.is_valid_molecule(mol_id)
+		]
+		(counts,) = read_bulk_molecule_counts(
+			simOutDir, (existing_molecule_ids,))
 		counts = counts.astype(float).T
 
 		row_max = counts.max(axis=1)
@@ -73,10 +81,10 @@ def do_plot(self, simOutDir, plotOutDir, plotOutFileName,
 			cmap='RdBu',
 			interpolation='nearest',
 			aspect='auto',
-			extent=[time[0],time[-1],len(cell_wall_proteins)-0.5,-0.5])
-		ax.set_yticks(np.arange(0, len(cell_wall_proteins), 1))
-		ax.set_yticklabels([cell_wall_proteins[mol_id] for
-							mol_id in cell_wall_proteins.keys()], fontsize=8)
+			extent=[time[0],time[-1],len(existing_molecule_ids)-0.5,-0.5])
+		ax.set_yticks(np.arange(0, len(existing_molecule_ids), 1))
+		ax.set_yticklabels([CELL_WALL_PROTEINS[mol_id[:-3]] for
+							mol_id in existing_molecule_ids], fontsize=8)
 		plt.xlabel('time (s)')
 		plt.title('Cell wall major components expression')
 
@@ -85,8 +93,8 @@ def do_plot(self, simOutDir, plotOutDir, plotOutFileName,
 		final_counts.reverse()
 
 		ax2 = fig.add_subplot(111, sharex=ax, frameon=False)
-		ax2.set_ylim([-0.5, len(cell_wall_proteins.keys())-0.5])
-		ax2.set_yticks(np.arange(0, len(cell_wall_proteins), 1))
+		ax2.set_ylim([-0.5, len(existing_molecule_ids)-0.5])
+		ax2.set_yticks(np.arange(0, len(existing_molecule_ids), 1))
 		ax2.set_yticklabels(final_counts, fontsize=8)
 		ax2.tick_params(length=0)
 		ax2.yaxis.tick_right()

models/ecoli/analysis/single/equilibriumComparison.py

@@ -48,8 +48,8 @@ def do_plot(self, simOutDir, plotOutDir, plotOutFileName, simDataFile, validatio
 		cellVolume = cellMass / cellDensity
 
 		fig = plt.figure(figsize = (20, 20))
-		rows = 6
-		cols = 6
+		rows = 8
+		cols = 8
 		num_subentries = 3
 
 		for idx in range(stoichMatrix.shape[1]):