"""
=====================
Transcript Initiation
=====================
This process models the binding of RNA polymerase to each gene.
The number of RNA polymerases to activate in each time step is determined
such that the average fraction of RNA polymerases that are active throughout
the simulation matches measured fractions, which are dependent on the
cellular growth rate. This is done by assuming a steady state concentration
of active RNA polymerases.
TODO:
- use transcription units instead of single genes
- match sigma factors to promoters
"""
import numpy as np
import scipy.sparse
import matplotlib.pyplot as plt
from typing import cast
from vivarium.core.composition import simulate_process
from ecoli.library.schema import (
create_unique_indexes,
listener_schema,
numpy_schema,
counts,
attrs,
bulk_name_to_idx,
)
from wholecell.utils import units
from wholecell.utils.random import stochasticRound
from wholecell.utils.unit_struct_array import UnitStructArray
from ecoli.library.data_predicates import monotonically_decreasing, all_nonnegative
from scipy.stats import chisquare
from ecoli.processes.registries import topology_registry
from ecoli.processes.partition import PartitionedProcess
# Register default topology for this process, associating it with process name
NAME = "ecoli-transcript-initiation"
TOPOLOGY = {
"environment": ("environment",),
"full_chromosomes": ("unique", "full_chromosome"),
"RNAs": ("unique", "RNA"),
"active_RNAPs": ("unique", "active_RNAP"),
"promoters": ("unique", "promoter"),
"bulk": ("bulk",),
"listeners": ("listeners",),
"timestep": ("timestep",),
}
topology_registry.register(NAME, TOPOLOGY)
[docs]
class TranscriptInitiation(PartitionedProcess):
"""Transcript Initiation PartitionedProcess
**Defaults:**
- **fracActiveRnapDict** (``dict``): Dictionary with keys corresponding to
media, values being the fraction of active RNA Polymerase (RNAP)
for that media.
- **rnaLengths** (``numpy.ndarray[int]``): lengths of RNAs for each transcription
unit (TU), in nucleotides
- **rnaPolymeraseElongationRateDict** (``dict``): Dictionary with keys
corresponding to media, values being RNAP's elongation rate in
that media, in nucleotides/s
- **variable_elongation** (``bool``): Whether to add amplified elongation rates
for rRNAs. False by default.
- **make_elongation_rates** (``func``): Function for making elongation rates
(see :py:meth:`~reconstruction.ecoli.dataclasses.process.transcription.Transcription.make_elongation_rates`).
Takes PRNG, basal elongation rate, timestep, and ``variable_elongation``.
Returns an array of elongation rates for all genes.
- **active_rnap_footprint_size** (``int``): Maximum physical footprint of RNAP
in nucleotides to cap initiation probabilities
- **basal_prob** (``numpy.ndarray[float]``): Baseline probability of synthesis for
every TU.
- **delta_prob** (``dict``): Dictionary with four keys, used to create a matrix
encoding the effect of transcription factors (TFs) on transcription
probabilities::
{'deltaV' (np.ndarray[float]): deltas associated with the effects of
TFs on TUs,
'deltaI' (np.ndarray[int]): index of the affected TU for each delta,
'deltaJ' (np.ndarray[int]): index of the acting TF for each delta,
'shape' (tuple): (m, n) = (# of TUs, # of TFs)}
- **perturbations** (``dict``): Dictionary of genetic perturbations (optional,
can be empty)
- **rna_data** (``numpy.ndarray``): Structured array with an entry for each TU,
where entries look like::
(id, deg_rate, length (nucleotides), counts_AGCU, mw
(molecular weight), is_mRNA, is_miscRNA, is_rRNA, is_tRNA,
is_23S_rRNA, is_16S_rRNA, is_5S_rRNA, is_ribosomal_protein,
is_RNAP, gene_id, Km_endoRNase, replication_coordinate,
direction)
- **idx_rRNA** (``numpy.ndarray[int]``): indexes of TUs corresponding to rRNAs
- **idx_mRNA** (``numpy.ndarray[int]``): indexes of TUs corresponding to mRNAs
- **idx_tRNA** (``numpy.ndarray[int]``): indexes of TUs corresponding to tRNAs
- **idx_rprotein** (``numpy.ndarray[int]``): indexes of TUs corresponding ribosomal
proteins
- **idx_rnap** (``numpy.ndarray[int]``): indexes of TU corresponding to RNAP
- **rnaSynthProbFractions** (``dict``): Dictionary where keys are media types,
values are sub-dictionaries with keys 'mRna', 'tRna', 'rRna', and
values being probabilities of synthesis for each RNA type
- **rnaSynthProbRProtein** (``dict``): Dictionary where keys are media types,
values are arrays storing the (fixed) probability of synthesis for
each rProtein TU, under that media condition.
- **rnaSynthProbRnaPolymerase** (``dict``): Dictionary where keys are media
types, values are arrays storing the (fixed) probability of
synthesis for each RNAP TU, under that media condition.
- **replication_coordinate** (``numpy.ndarray[int]``): Array with chromosome
coordinates for each TU
- **transcription_direction** (``numpy.ndarray[bool]``): Array of transcription
directions for each TU
- **n_avogadro** (``unum.Unum``): Avogadro's number (constant)
- **cell_density** (``unum.Unum``): Density of cell (constant)
- **ppgpp** (``str``): id of ppGpp
- **inactive_RNAP** (``str``): id of inactive RNAP
- **synth_prob** (``Callable[[Unum, int], numpy.ndarrray[float]]``):
Function used in model of ppGpp regulation
(see :py:func:`~reconstruction.ecoli.dataclasses.process.transcription.Transcription.synth_prob_from_ppgpp`).
Takes ppGpp concentration (mol/volume) and copy number, returns
normalized synthesis probability for each gene
- **copy_number** (``Callable[[Unum, int], numpy.ndarrray[float]]``):
see :py:func:`~reconstruction.ecoli.dataclasses.process.replication.Replication.get_average_copy_number`.
Takes expected doubling time in minutes and chromosome coordinates of genes,
returns average copy number of each gene expected at doubling time
- **ppgpp_regulation** (``bool``): Whether to include model of ppGpp regulation
- **get_rnap_active_fraction_from_ppGpp** (``Callable[[Unum], float]``):
Returns elongation rate for a given ppGpp concentration
- **seed** (``int``): random seed to initialize PRNG
"""
name = NAME
topology = TOPOLOGY
defaults = {
"fracActiveRnapDict": {},
"rnaLengths": np.array([]),
"rnaPolymeraseElongationRateDict": {},
"variable_elongation": False,
"make_elongation_rates": (
lambda random, rate, timestep, variable: np.array([])
),
"active_rnap_foorprint_size": 1,
"basal_prob": np.array([]),
"delta_prob": {"deltaI": [], "deltaJ": [], "deltaV": [], "shape": tuple()},
"get_delta_prob_matrix": None,
"perturbations": {},
"rna_data": {},
"active_rnap_footprint_size": 24 * units.nt,
"get_rnap_active_fraction_from_ppGpp": lambda x: 0.1,
"idx_rRNA": np.array([]),
"idx_mRNA": np.array([]),
"idx_tRNA": np.array([]),
"idx_rprotein": np.array([]),
"idx_rnap": np.array([]),
"rnaSynthProbFractions": {},
"rnaSynthProbRProtein": {},
"rnaSynthProbRnaPolymerase": {},
"replication_coordinate": np.array([]),
"transcription_direction": np.array([]),
"n_avogadro": 6.02214076e23 / units.mol,
"cell_density": 1100 * units.g / units.L,
"ppgpp": "ppGpp",
"inactive_RNAP": "APORNAP-CPLX[c]",
"synth_prob": lambda concentration, copy: 0.0,
"copy_number": lambda x: x,
"ppgpp_regulation": False,
# attenuation
"trna_attenuation": False,
"attenuated_rna_indices": np.array([]),
"attenuation_adjustments": np.array([]),
# random seed
"seed": 0,
"emit_unique": False,
}
# Constructor
def __init__(self, parameters=None):
super().__init__(parameters)
# Load parameters
self.fracActiveRnapDict = self.parameters["fracActiveRnapDict"]
self.rnaLengths = self.parameters["rnaLengths"]
self.rnaPolymeraseElongationRateDict = self.parameters[
"rnaPolymeraseElongationRateDict"
]
self.variable_elongation = self.parameters["variable_elongation"]
self.make_elongation_rates = self.parameters["make_elongation_rates"]
self.active_rnap_footprint_size = self.parameters["active_rnap_footprint_size"]
# Initialize matrices used to calculate synthesis probabilities
self.basal_prob = self.parameters["basal_prob"].copy()
self.trna_attenuation = self.parameters["trna_attenuation"]
if self.trna_attenuation:
self.attenuated_rna_indices = self.parameters["attenuated_rna_indices"]
self.attenuation_adjustments = self.parameters["attenuation_adjustments"]
self.basal_prob[self.attenuated_rna_indices] += self.attenuation_adjustments
self.n_TUs = len(self.basal_prob)
self.delta_prob = self.parameters["delta_prob"]
if self.parameters["get_delta_prob_matrix"] is not None:
self.delta_prob_matrix = self.parameters["get_delta_prob_matrix"](
dense=True, ppgpp=self.parameters["ppgpp_regulation"]
)
else:
# make delta_prob_matrix without adjustments
self.delta_prob_matrix = scipy.sparse.csr_matrix(
(
self.delta_prob["deltaV"],
(self.delta_prob["deltaI"], self.delta_prob["deltaJ"]),
),
shape=self.delta_prob["shape"],
).toarray()
# Determine changes from genetic perturbations
self.genetic_perturbations = {}
self.perturbations = self.parameters["perturbations"]
self.rna_data = self.parameters["rna_data"]
if len(self.perturbations) > 0:
probability_indexes = [
(index, self.perturbations[rna_data["id"]])
for index, rna_data in enumerate(self.rna_data)
if rna_data["id"] in self.perturbations
]
self.genetic_perturbations = {
"fixedRnaIdxs": [pair[0] for pair in probability_indexes],
"fixedSynthProbs": [pair[1] for pair in probability_indexes],
}
# ID Groups
self.idx_rRNA = self.parameters["idx_rRNA"]
self.idx_mRNA = self.parameters["idx_mRNA"]
self.idx_tRNA = self.parameters["idx_tRNA"]
self.idx_rprotein = self.parameters["idx_rprotein"]
self.idx_rnap = self.parameters["idx_rnap"]
# Synthesis probabilities for different categories of genes
self.rnaSynthProbFractions = self.parameters["rnaSynthProbFractions"]
self.rnaSynthProbRProtein = self.parameters["rnaSynthProbRProtein"]
self.rnaSynthProbRnaPolymerase = self.parameters["rnaSynthProbRnaPolymerase"]
# Coordinates and transcription directions of transcription units
self.replication_coordinate = self.parameters["replication_coordinate"]
self.transcription_direction = self.parameters["transcription_direction"]
self.inactive_RNAP = self.parameters["inactive_RNAP"]
# ppGpp control related
self.n_avogadro = self.parameters["n_avogadro"]
self.cell_density = self.parameters["cell_density"]
self.ppgpp = self.parameters["ppgpp"]
self.synth_prob = self.parameters["synth_prob"]
self.copy_number = self.parameters["copy_number"]
self.ppgpp_regulation = self.parameters["ppgpp_regulation"]
self.get_rnap_active_fraction_from_ppGpp = self.parameters[
"get_rnap_active_fraction_from_ppGpp"
]
self.seed = self.parameters["seed"]
self.random_state = np.random.RandomState(seed=self.seed)
# Use separate random state instance to create unique indices
# so results are directly comparable with wcEcoli
self.unique_idx_random_state = np.random.RandomState(seed=self.seed)
# Helper indices for Numpy indexing
self.ppgpp_idx = None
[docs]
def ports_schema(self):
return {
"environment": {"media_id": {"_default": "", "_updater": "set"}},
"bulk": numpy_schema("bulk"),
"full_chromosomes": numpy_schema(
"full_chromosomes", emit=self.parameters["emit_unique"]
),
"promoters": numpy_schema("promoters", emit=self.parameters["emit_unique"]),
"RNAs": numpy_schema("RNAs", emit=self.parameters["emit_unique"]),
"active_RNAPs": numpy_schema(
"active_RNAPs", emit=self.parameters["emit_unique"]
),
"listeners": {
"mass": {"cell_mass": {"_default": 0.0}, "dry_mass": {"_default": 0.0}},
"rna_synth_prob": listener_schema(
{
"target_rna_synth_prob": [0.0],
"actual_rna_synth_prob": [0.0],
"tu_is_overcrowded": (
[False] * self.n_TUs,
self.rna_data["id"],
),
"total_rna_init": 0,
"max_p": 0.0,
}
),
"ribosome_data": listener_schema(
{
"rRNA_initiated_TU": [0] * len(self.idx_rRNA),
"rRNA_init_prob_TU": [0.0] * len(self.idx_rRNA),
"total_rna_init": 0,
}
),
"rnap_data": listener_schema(
{"did_initialize": 0, "rna_init_event": (0, self.rna_data["id"])}
),
},
"timestep": {"_default": self.parameters["time_step"]},
}
[docs]
def calculate_request(self, timestep, states):
# At first update, convert all strings to indices
if self.ppgpp_idx is None:
bulk_ids = states["bulk"]["id"]
self.ppgpp_idx = bulk_name_to_idx(self.ppgpp, bulk_ids)
self.inactive_RNAP_idx = bulk_name_to_idx(self.inactive_RNAP, bulk_ids)
# Get all inactive RNA polymerases
requests = {
"bulk": [
(self.inactive_RNAP_idx, counts(states["bulk"], self.inactive_RNAP_idx))
]
}
# Read current environment
current_media_id = states["environment"]["media_id"]
if states["full_chromosomes"]["_entryState"].sum() > 0:
# Get attributes of promoters
TU_index, bound_TF = attrs(states["promoters"], ["TU_index", "bound_TF"])
if self.ppgpp_regulation:
cell_mass = states["listeners"]["mass"]["cell_mass"] * units.fg
cell_volume = cell_mass / self.cell_density
counts_to_molar = 1 / (self.n_avogadro * cell_volume)
ppgpp_conc = counts(states["bulk"], self.ppgpp_idx) * counts_to_molar
basal_prob, _ = self.synth_prob(ppgpp_conc, self.copy_number)
if self.trna_attenuation:
basal_prob[self.attenuated_rna_indices] += (
self.attenuation_adjustments
)
self.fracActiveRnap = self.get_rnap_active_fraction_from_ppGpp(
ppgpp_conc
)
ppgpp_scale = basal_prob[TU_index]
# Use original delta prob if no ppGpp basal
ppgpp_scale[ppgpp_scale == 0] = 1
else:
basal_prob = self.basal_prob
self.fracActiveRnap = self.fracActiveRnapDict[current_media_id]
ppgpp_scale = 1
# Calculate probabilities of the RNAP binding to each promoter
self.promoter_init_probs = basal_prob[TU_index] + ppgpp_scale * np.multiply(
self.delta_prob_matrix[TU_index, :], bound_TF
).sum(axis=1)
if len(self.genetic_perturbations) > 0:
self._rescale_initiation_probs(
self.genetic_perturbations["fixedRnaIdxs"],
self.genetic_perturbations["fixedSynthProbs"],
TU_index,
)
# Adjust probabilities to not be negative
self.promoter_init_probs[self.promoter_init_probs < 0] = 0.0
self.promoter_init_probs /= self.promoter_init_probs.sum()
if not self.ppgpp_regulation:
# Adjust synthesis probabilities depending on environment
synthProbFractions = self.rnaSynthProbFractions[current_media_id]
# Create masks for different types of RNAs
is_mrna = np.isin(TU_index, self.idx_mRNA)
is_trna = np.isin(TU_index, self.idx_tRNA)
is_rrna = np.isin(TU_index, self.idx_rRNA)
is_rprotein = np.isin(TU_index, self.idx_rprotein)
is_rnap = np.isin(TU_index, self.idx_rnap)
is_fixed = is_trna | is_rrna | is_rprotein | is_rnap
# Rescale initiation probabilities based on type of RNA
self.promoter_init_probs[is_mrna] *= (
synthProbFractions["mRna"] / self.promoter_init_probs[is_mrna].sum()
)
self.promoter_init_probs[is_trna] *= (
synthProbFractions["tRna"] / self.promoter_init_probs[is_trna].sum()
)
self.promoter_init_probs[is_rrna] *= (
synthProbFractions["rRna"] / self.promoter_init_probs[is_rrna].sum()
)
# Set fixed synthesis probabilities for RProteins and RNAPs
self._rescale_initiation_probs(
np.concatenate((self.idx_rprotein, self.idx_rnap)),
np.concatenate(
(
self.rnaSynthProbRProtein[current_media_id],
self.rnaSynthProbRnaPolymerase[current_media_id],
)
),
TU_index,
)
assert self.promoter_init_probs[is_fixed].sum() < 1.0
# Scale remaining synthesis probabilities accordingly
scaleTheRestBy = (
1.0 - self.promoter_init_probs[is_fixed].sum()
) / self.promoter_init_probs[~is_fixed].sum()
self.promoter_init_probs[~is_fixed] *= scaleTheRestBy
# If there are no chromosomes in the cell, set all probs to zero
else:
self.promoter_init_probs = np.zeros(
states["promoters"]["_entryState"].sum()
)
self.rnaPolymeraseElongationRate = self.rnaPolymeraseElongationRateDict[
current_media_id
]
self.elongation_rates = self.make_elongation_rates(
self.random_state,
self.rnaPolymeraseElongationRate.asNumber(units.nt / units.s),
1, # want elongation rate, not lengths adjusted for time step
self.variable_elongation,
)
return requests
[docs]
def evolve_state(self, timestep, states):
update = {
"listeners": {
"rna_synth_prob": {
"target_rna_synth_prob": np.zeros(self.n_TUs),
"actual_rna_synth_prob": np.zeros(self.n_TUs),
"tu_is_overcrowded": np.zeros(self.n_TUs, dtype=np.bool_),
"total_rna_init": 0,
"max_p": 0.0,
},
"ribosome_data": {"total_rna_init": 0},
"rnap_data": {
"did_initialize": 0,
"rna_init_event": np.zeros(self.n_TUs, dtype=np.int64),
},
},
"active_RNAPs": {},
"full_chromosomes": {},
"promoters": {},
"RNAs": {},
}
# no synthesis if no chromosome
if len(states["full_chromosomes"]) == 0:
return update
# Get attributes of promoters
TU_index, domain_index_promoters = attrs(
states["promoters"], ["TU_index", "domain_index"]
)
n_promoters = states["promoters"]["_entryState"].sum()
# Construct matrix that maps promoters to transcription units
TU_to_promoter = scipy.sparse.csr_matrix(
(np.ones(n_promoters), (TU_index, np.arange(n_promoters))),
shape=(self.n_TUs, n_promoters),
dtype=np.int8,
)
# Compute target synthesis probabilities of each transcription unit
target_TU_synth_probs = TU_to_promoter.dot(self.promoter_init_probs)
update["listeners"]["rna_synth_prob"]["target_rna_synth_prob"] = (
target_TU_synth_probs
)
# Calculate RNA polymerases to activate based on probabilities
# Note: ideally we should be using the actual TU synthesis probabilities
# here, but the calculation of actual probabilities requires the number
# of RNAPs to activate. The difference should be small.
self.activationProb = self._calculateActivationProb(
states["timestep"],
self.fracActiveRnap,
self.rnaLengths,
(units.nt / units.s) * self.elongation_rates,
target_TU_synth_probs,
)
n_RNAPs_to_activate = np.int64(
self.activationProb * counts(states["bulk"], self.inactive_RNAP_idx)
)
if n_RNAPs_to_activate == 0:
return update
# Cap the initiation probabilities at the maximum level physically
# allowed from the known RNAP footprint sizes
max_p = (
self.rnaPolymeraseElongationRate
/ self.active_rnap_footprint_size
* (units.s)
* states["timestep"]
/ n_RNAPs_to_activate
).asNumber()
update["listeners"]["rna_synth_prob"]["max_p"] = max_p
is_overcrowded = self.promoter_init_probs > max_p
while np.any(self.promoter_init_probs > max_p):
self.promoter_init_probs[is_overcrowded] = max_p
scale_the_rest_by = (
1.0 - self.promoter_init_probs[is_overcrowded].sum()
) / self.promoter_init_probs[~is_overcrowded].sum()
self.promoter_init_probs[~is_overcrowded] *= scale_the_rest_by
is_overcrowded |= self.promoter_init_probs > max_p
# Compute actual synthesis probabilities of each transcription unit
actual_TU_synth_probs = TU_to_promoter.dot(self.promoter_init_probs)
tu_is_overcrowded = TU_to_promoter.dot(is_overcrowded).astype(bool)
update["listeners"]["rna_synth_prob"]["actual_rna_synth_prob"] = (
actual_TU_synth_probs
)
update["listeners"]["rna_synth_prob"]["tu_is_overcrowded"] = tu_is_overcrowded
# Sample a multinomial distribution of initiation probabilities to
# determine what promoters are initialized
n_initiations = self.random_state.multinomial(
n_RNAPs_to_activate, self.promoter_init_probs
)
# Build array of transcription unit indexes for partially transcribed
# RNAs and domain indexes for RNAPs
TU_index_partial_RNAs = np.repeat(TU_index, n_initiations)
domain_index_rnap = np.repeat(domain_index_promoters, n_initiations)
# Build arrays of starting coordinates and transcription directions
coordinates = self.replication_coordinate[TU_index_partial_RNAs]
is_forward = self.transcription_direction[TU_index_partial_RNAs]
# new RNAPs
RNAP_indexes = create_unique_indexes(
n_RNAPs_to_activate, self.unique_idx_random_state
)
update["active_RNAPs"].update(
{
"add": {
"unique_index": RNAP_indexes,
"domain_index": domain_index_rnap,
"coordinates": coordinates,
"is_forward": is_forward,
}
}
)
# Decrement counts of inactive RNAPs
update["bulk"] = [(self.inactive_RNAP_idx, -n_initiations.sum())]
# Add partially transcribed RNAs
is_mRNA = np.isin(TU_index_partial_RNAs, self.idx_mRNA)
rna_indices = create_unique_indexes(
n_RNAPs_to_activate, self.unique_idx_random_state
)
update["RNAs"].update(
{
"add": {
"unique_index": rna_indices,
"TU_index": TU_index_partial_RNAs,
"transcript_length": np.zeros(cast(int, n_RNAPs_to_activate)),
"is_mRNA": is_mRNA,
"is_full_transcript": np.zeros(
cast(int, n_RNAPs_to_activate), dtype=bool
),
"can_translate": is_mRNA,
"RNAP_index": RNAP_indexes,
}
}
)
rna_init_event = TU_to_promoter.dot(n_initiations)
rRNA_initiations = rna_init_event[self.idx_rRNA]
# Write outputs to listeners
update["listeners"]["ribosome_data"] = {
"rRNA_initiated_TU": rRNA_initiations.astype(int),
"rRNA_init_prob_TU": rRNA_initiations / float(n_RNAPs_to_activate),
"total_rna_init": n_RNAPs_to_activate,
}
update["listeners"]["rnap_data"] = {
"did_initialize": n_RNAPs_to_activate,
"rna_init_event": rna_init_event.astype(np.int64),
}
update["listeners"]["rna_synth_prob"]["total_rna_init"] = n_RNAPs_to_activate
return update
[docs]
def _calculateActivationProb(
self,
timestep,
fracActiveRnap,
rnaLengths,
rnaPolymeraseElongationRates,
synthProb,
):
"""
Calculate expected RNAP termination rate based on RNAP elongation rate
- allTranscriptionTimes: Vector of times required to transcribe each
transcript
- allTranscriptionTimestepCounts: Vector of numbers of timesteps
required to transcribe each transcript
- averageTranscriptionTimeStepCounts: Average number of timesteps
required to transcribe a transcript, weighted by synthesis
probabilities of each transcript
- expectedTerminationRate: Average number of terminations in one
timestep for one transcript
"""
allTranscriptionTimes = 1.0 / rnaPolymeraseElongationRates * rnaLengths
timesteps = (1.0 / (timestep * units.s) * allTranscriptionTimes).asNumber()
allTranscriptionTimestepCounts = np.ceil(timesteps)
averageTranscriptionTimestepCounts = np.dot(
synthProb, allTranscriptionTimestepCounts
)
expectedTerminationRate = 1.0 / averageTranscriptionTimestepCounts
"""
Modify given fraction of active RNAPs to take into account early
terminations in between timesteps
- allFractionTimeInactive: Vector of probabilities an "active" RNAP
will in effect be "inactive" because it has terminated during a
timestep
- averageFractionTimeInactive: Average probability of an "active" RNAP
being in effect "inactive", weighted by synthesis probabilities
- effectiveFracActiveRnap: New higher "goal" for fraction of active
RNAP, considering that the "effective" fraction is lower than what the
listener sees
"""
allFractionTimeInactive = (
1
- (1.0 / (timestep * units.s) * allTranscriptionTimes).asNumber()
/ allTranscriptionTimestepCounts
)
averageFractionTimeInactive = np.dot(allFractionTimeInactive, synthProb)
effectiveFracActiveRnap = fracActiveRnap / (1 - averageFractionTimeInactive)
# Return activation probability that will balance out the expected termination rate
activation_prob = (
effectiveFracActiveRnap
* expectedTerminationRate
/ (1 - effectiveFracActiveRnap)
)
if activation_prob > 1:
activation_prob = 1
return activation_prob
[docs]
def _rescale_initiation_probs(self, fixed_indexes, fixed_synth_probs, TU_index):
"""
Rescales the initiation probabilities of each promoter such that the
total synthesis probabilities of certain types of RNAs are fixed to
a predetermined value. For instance, if there are two copies of
promoters for RNA A, whose synthesis probability should be fixed to
0.1, each promoter is given an initiation probability of 0.05.
"""
for idx, synth_prob in zip(fixed_indexes, fixed_synth_probs):
fixed_mask = TU_index == idx
self.promoter_init_probs[fixed_mask] = synth_prob / fixed_mask.sum()
def test_transcript_initiation(return_data=False):
def make_elongation_rates(random, base, time_step, variable_elongation=False):
size = 9 # number of TUs
lengths = time_step * np.full(size, base, dtype=np.int64)
lengths = stochasticRound(random, lengths) if random else np.round(lengths)
return lengths.astype(np.int64)
rna_data = UnitStructArray(
# id, deg_rate, len, counts, _ACGU mw, mRNA?, miscRNA?, rRNA?, tRNA?, 23S?, 16S?, 5S?, rProt?, RNAP?, geneid,
# Km, coord, direction
np.array(
[
(
"16SrRNA",
0.002,
45,
[10, 11, 12, 12],
13500,
False,
False,
True,
False,
False,
True,
False,
False,
False,
"16SrRNA",
2e-4,
0,
True,
),
(
"23SrRNA",
0.002,
450,
[100, 110, 120, 120],
135000,
False,
False,
True,
False,
True,
False,
False,
False,
False,
"23SrRNA",
2e-4,
1000,
True,
),
(
"5SrRNA",
0.002,
600,
[150, 150, 150, 150],
180000,
False,
False,
True,
False,
False,
False,
True,
False,
False,
"5SrRNA",
2e-4,
2000,
True,
),
(
"rProtein",
0.002,
700,
[175, 175, 175, 175],
210000,
True,
False,
False,
False,
False,
False,
False,
True,
False,
"rProtein",
2e-4,
3000,
False,
),
(
"RNAP",
0.002,
800,
[200, 200, 200, 200],
240000,
True,
False,
False,
False,
False,
False,
False,
False,
True,
"RNAP",
2e-4,
4000,
False,
),
(
"miscProt",
0.002,
900,
[225, 225, 225, 225],
270000,
True,
False,
False,
False,
False,
False,
False,
False,
False,
"miscProt",
2e-4,
5000,
True,
),
(
"tRNA1",
0.002,
1200,
[300, 300, 300, 300],
360000,
False,
False,
False,
True,
False,
False,
False,
False,
False,
"tRNA1",
2e-4,
6000,
False,
),
(
"tRNA2",
0.002,
4000,
[1000, 1000, 1000, 1000],
1200000,
False,
False,
False,
True,
False,
False,
False,
False,
False,
"tRNA2",
2e-4,
7000,
False,
),
(
"tRNA3",
0.002,
7000,
[1750, 1750, 1750, 1750],
2100000,
False,
False,
False,
True,
False,
False,
False,
False,
False,
"tRNA3",
2e-4,
8000,
True,
),
],
dtype=[
("id", "<U15"),
("deg_rate", "<f8"),
("length", "<i8"),
("counts_ACGU", "<i8", (4,)),
("mw", "<f8"),
("is_mRNA", "?"),
("is_miscRNA", "?"),
("is_rRNA", "?"),
("is_tRNA", "?"),
("is_23S_rRNA", "?"),
("is_16S_rRNA", "?"),
("is_5S_rRNA", "?"),
("is_ribosomal_protein", "?"),
("is_RNAP", "?"),
("gene_id", "<U8"),
("Km_endoRNase", "<f8"),
("replication_coordinate", "<i8"),
("direction", "?"),
],
),
{
"id": None,
"deg_rate": 1.0 / units.s,
"length": units.nt,
"counts_ACGU": units.nt,
"mw": units.g / units.mol,
"is_mRNA": None,
"is_miscRNA": None,
"is_rRNA": None,
"is_tRNA": None,
"is_23S_rRNA": None,
"is_16S_rRNA": None,
"is_5S_rRNA": None,
"is_ribosomal_protein": None,
"is_RNAP": None,
"gene_id": None,
"Km_endoRNase": units.mol / units.L,
"replication_coordinate": None,
"direction": None,
},
)
test_config = {
"fracActiveRnapDict": {"minimal": 0.2},
"rnaLengths": np.array([x[2] for x in rna_data.fullArray()]),
"rnaPolymeraseElongationRateDict": {"minimal": 50 * units.nt / units.s},
"make_elongation_rates": make_elongation_rates,
"basal_prob": np.array([1e-7, 1e-7, 1e-7, 1e-7, 1e-6, 1e-6, 1e-6, 1e-5, 1e-5]),
"delta_prob": {
"deltaV": [-1e-3, -1e-5, -1e-6, 1e-7, 1e-6, 1e-6, 1e-5],
"deltaI": [0, 1, 2, 3, 4, 5, 6],
"deltaJ": [0, 1, 2, 3, 0, 1, 2],
"shape": (9, 4),
},
"rna_data": rna_data,
"idx_16SrRNA": np.array([0]),
"idx_23SrRNA": np.array([1]),
"idx_5SrRNA": np.array([2]),
"idx_rRNA": np.array([0, 1, 2]),
"idx_mRNA": np.array([3, 4, 5]),
"idx_tRNA": np.array([6, 7, 8]),
"idx_rprotein": np.array([3]),
"idx_rnap": np.array([4]),
"rnaSynthProbFractions": {"minimal": {"mRna": 0.25, "tRna": 0.6, "rRna": 0.15}},
"rnaSynthProbRProtein": {"minimal": np.array([0.06])},
"rnaSynthProbRnaPolymerase": {"minimal": np.array([0.04])},
"replication_coordinate": np.array([x[-2] for x in rna_data.fullArray()]),
"transcription_direction": np.array([x[-1] for x in rna_data.fullArray()]),
"inactive_RNAP": "APORNAP-CPLX[c]",
"seed": 0,
"time_step": 2,
#'_schema' : {'molecules' : {'APORNAP-CPLX[c]' : {'_updater' : 'null'}}}
}
transcript_initiation = TranscriptInitiation(test_config)
submass_dtypes = [
("massDiff_DNA", "<f8"),
("massDiff_mRNA", "<f8"),
("massDiff_metabolite", "<f8"),
("massDiff_miscRNA", "<f8"),
("massDiff_nonspecific_RNA", "<f8"),
("massDiff_protein", "<f8"),
("massDiff_rRNA", "<f8"),
("massDiff_tRNA", "<f8"),
("massDiff_water", "<f8"),
]
chromosome_dtypes = [
("_entryState", "i1"),
("_globalIndex", "<i8"),
("division_time", "<f8"),
("domain_index", "<i4"),
("has_triggered_division", "?"),
("unique_index", "<i8"),
]
rna_dtypes = [
("RNAP_index", "<i8"),
("TU_index", "<i8"),
("_entryState", "i1"),
("_globalIndex", "<i8"),
("can_translate", "?"),
("is_full_transcript", "?"),
("is_mRNA", "?"),
("unique_index", "<i8"),
("transcript_length", "<i8"),
]
active_rnap_dtypes = [
("_entryState", "i1"),
("_globalIndex", "<i8"),
("coordinates", "<i8"),
("is_forward", "?"),
("domain_index", "<i4"),
("unique_index", "<i8"),
]
promoter_dtypes = [
("TU_index", "<i8"),
("_entryState", "i1"),
("_globalIndex", "<i8"),
("bound_TF", "?", (4,)),
("coordinates", "<i8"),
("domain_index", "<i4"),
("unique_index", "<i8"),
]
initial_state = {
"environment": {"media_id": "minimal"},
"bulk": np.array(
[("APORNAP-CPLX[c]", 1000), ("GUANOSINE-5DP-3DP[c]", 0), ("ppGpp", 0)],
dtype=[("id", "U40"), ("count", int)],
),
"listeners": {"mass": {"cell_mass": 1000, "dry_mass": 350}},
}
unique_state = {
"full_chromosome": np.array(
[
(1, 0, 0, 0, False, 0) + (0,) * 9,
],
dtype=chromosome_dtypes + submass_dtypes,
),
"promoter": np.array(
[
(i, 1, i, [False] * 4, subdata["replication_coordinate"], 0, i)
+ (0,) * 9
for i, subdata in enumerate(rna_data)
],
dtype=promoter_dtypes + submass_dtypes,
),
"RNA": np.array([], dtype=rna_dtypes + submass_dtypes),
"active_RNAP": np.array([], dtype=active_rnap_dtypes + submass_dtypes),
}
initial_state["unique"] = unique_state
settings = {"total_time": 100, "initial_state": initial_state, "topology": TOPOLOGY}
data_noTF = simulate_process(transcript_initiation, settings)
# Also gather data where TFs are bound:
# Assertions =========================================================
# TODO(Michael):
# 1) When no initiations occurred in a timestep, the inits_by_TU is a scalar 0
# 2) Weird things happen when RNAP is limiting, including affecting RNA synth probs
# - for toy model, initial RNAPs <= 132 results in no initiation
# - 1000 initial RNAPs, run for t=1000s results in initiations without depletion of inactive RNAP
# 3) rnaps in data['active_RNAPs'] seem to be missing direction, domain index data (likewise with rnas)
# 4) Effect of TF binding is not tested in the toy model
# - simple test is to compare gene under up/down-regulation with data from same gene without regulation
# 5) Test of fixed synthesis probabilties does not pass
# Unpack data
bulk_timeseries = np.array(data_noTF["bulk"])
inactive_RNAP = bulk_timeseries[:, 0]
d_inactive_RNAP = inactive_RNAP[1:] - inactive_RNAP[:-1]
d_active_RNAP = np.array(data_noTF["listeners"]["rnap_data"]["did_initialize"][1:])
inits_by_TU = np.stack(data_noTF["listeners"]["rnap_data"]["rna_init_event"][1:])
rnap_inits = inits_by_TU[:, test_config["idx_rnap"]]
rprotein_inits = inits_by_TU[:, test_config["idx_rprotein"]]
# Sanity checks
assert monotonically_decreasing(
inactive_RNAP
), "Inactive RNAPs are not monotonically decreasing"
assert all_nonnegative(inactive_RNAP), "Inactive RNAPs fall below zero"
assert all_nonnegative(inits_by_TU), "Negative initiations (!?)"
assert monotonically_decreasing(
d_active_RNAP
), "Change in active RNAPs is not monotonically decreasing"
assert all_nonnegative(
d_active_RNAP
), "One or more timesteps has decrease in active RNAPs"
assert np.sum(d_active_RNAP) == np.sum(
inits_by_TU
), "# of active RNAPs does not match number of initiations"
# Inactive RNAPs deplete as they are activated
np.testing.assert_array_equal(
-d_inactive_RNAP,
d_active_RNAP,
"Depletion of inactive RNAPs does not match counts of RNAPs activated.",
)
# Fixed synthesis probability TUs (RNAP, rProtein) and non-fixed TUs synthesized in correct proportion
expected = np.array(
[
test_config["rnaSynthProbRProtein"]["minimal"][0],
test_config["rnaSynthProbRnaPolymerase"]["minimal"][0],
1,
]
)
expected[2] -= expected[0] + expected[1]
expected *= np.sum(inits_by_TU)
actual = np.array(
[
np.sum(rprotein_inits),
np.sum(rnap_inits),
np.sum(inits_by_TU) - np.sum(rprotein_inits) - np.sum(rnap_inits),
]
)
fixed_prob_test = chisquare(actual, f_exp=expected)
assert fixed_prob_test.pvalue > 0.05, (
"Distribution of RNA types synthesized does "
"not (approximately) match set points for fixed synthesis"
f"(p = {fixed_prob_test.pvalue} <= 0.05)"
)
# mRNA, tRNA, rRNA synthesized in correct proportion
RNA_dist = np.array(
[
np.sum(inits_by_TU[:, test_config["idx_mRNA"]]),
np.sum(inits_by_TU[:, test_config["idx_tRNA"]]),
np.sum(inits_by_TU[:, test_config["idx_rRNA"]]),
]
)
RNA_dist = RNA_dist / RNA_dist.sum()
RNA_synth_prob_test = chisquare(
RNA_dist, [v for v in test_config["rnaSynthProbFractions"]["minimal"].values()]
)
assert RNA_synth_prob_test.pvalue > 0.05, (
"Distribution of RNA types synthesized does"
"not (approximately) match values set by media"
)
if return_data:
return test_config, data_noTF
[docs]
def run_plot(config, data):
N = len(data["time"])
timestep = config["time_step"]
inits_by_TU = np.stack(data["listeners"]["rnap_data"]["rna_init_event"][1:])
synth_probs = np.array(
data["listeners"]["rna_synth_prob"]["actual_rna_synth_prob"][1:]
)
# plot synthesis probabilities over time
plt.subplot(2, 2, 1)
prev = np.zeros(N - 1)
for TU in range(synth_probs.shape[1]):
plt.bar(data["time"][1:], synth_probs[:, TU], bottom=prev, width=timestep)
prev += synth_probs[:, TU]
plt.xlabel("Time (s)")
plt.ylabel("Probability of Synthesis")
plt.title("Theoretical Synthesis Probabilities over Time")
# plot actual probability of synthesis for each RNA
plt.subplot(2, 2, 2)
probs = np.sum(inits_by_TU, axis=0) / np.sum(inits_by_TU)
prev = 0
for i in range(len(probs)):
prob = probs[i]
plt.bar([0], [prob], bottom=prev, width=1)
plt.text(i / len(probs) - 0.5, prev + prob / 2, config["rna_data"][i][0])
prev += prob
plt.tick_params(axis="x", which="both", bottom=False, top=False, labelbottom=False)
plt.ylabel("Probability")
plt.title("Actual Probability of Synthesis by TU")
# plot which RNAs are transcribed
plt.subplot(2, 2, 3)
prev = np.zeros(N - 1)
for TU in range(inits_by_TU.shape[1]):
plt.bar(data["time"][1:], inits_by_TU[:, TU], bottom=prev, width=timestep)
prev += inits_by_TU[:, TU]
plt.xlabel("Time (s)")
plt.ylabel("Transcripts")
plt.title("Transcripts over time, for all TUs")
# plot which RNAs are transcribed, grouped by mRNA/tRNA/rRNA
plt.subplot(2, 2, 4)
grouped_inits = np.concatenate(
[
[np.sum(inits_by_TU[:, config["idx_tRNA"]], axis=1)],
[np.sum(inits_by_TU[:, config["idx_mRNA"]], axis=1)],
[np.sum(inits_by_TU[:, config["idx_rRNA"]], axis=1)],
]
).T
prev = np.zeros(N - 1)
for TU in range(grouped_inits.shape[1]):
plt.bar(data["time"][1:], grouped_inits[:, TU], bottom=prev, width=timestep)
prev += grouped_inits[:, TU]
plt.xlabel("Time (s)")
plt.ylabel("Transcripts")
plt.title("Transcripts over time by RNA type")
plt.legend(
["tRNA", "mRNA", "rRNA"],
bbox_to_anchor=(1.05, 1),
loc="upper left",
borderaxespad=0.0,
)
plt.subplots_adjust(hspace=0.5, wspace=0.25)
plt.gcf().set_size_inches(10, 6)
plt.savefig("out/migration/transcript_initiation_toy_model.png")
[docs]
def main():
config, data = test_transcript_initiation(return_data=True)
run_plot(config, data)
if __name__ == "__main__":
main()