Modelling capture efficiency of single-cell RNA-sequencing data improves inference of transcriptome-wide burst kinetics.

MOTIVATION: Gene expression is characterized by stochastic bursts of transcription that occur at brief and random periods of promoter activity. The kinetics of gene expression burstiness differs across the genome and is dependent on the promoter sequence, among other factors. Single-cell RNA sequencing (scRNA-seq) has made it possible to quantify the cell-to-cell variability in transcription at a global genome-wide level. However, scRNA-seq data are prone to technical variability, including low and variable capture efficiency of transcripts from individual cells. RESULTS: Here, we propose a novel mathematical theory for the observed variability in scRNA-seq data. Our method captures burst kinetics and variability in both the cell size and capture efficiency, which allows us to propose several likelihood-based and simulation-based methods for the inference of burst kinetics from scRNA-seq data. Using both synthetic and real data, we show that the simulation-based methods provide an accurate, robust and flexible tool for inferring burst kinetics from scRNA-seq data. In particular, in a supervised manner, a simulation-based inference method based on neural networks proves to be accurate and useful when applied to both allele and nonallele-specific scRNA-seq data. AVAILABILITY AND IMPLEMENTATION: The code for Neural Network and Approximate Bayesian Computation inference is available at https://github.com/WT215/nnRNA and https://github.com/WT215/Julia_ABC, respectively.

A systematic dissection of determinants and consequences of snoRNA-guided pseudouridylation of human mRNA.

RNA can be extensively modified post-transcriptionally with >170 covalent modifications, expanding its functional and structural repertoire. Pseudouridine (Ψ), the most abundant modified nucleoside in rRNA and tRNA, has recently been found within mRNA molecules. It remains unclear whether pseudouridylation of mRNA can be snoRNA-guided, bearing important implications for understanding the physiological target spectrum of snoRNAs and for their potential therapeutic exploitation in genetic diseases. Here, using a massively parallel reporter based strategy we simultaneously interrogate Ψ levels across hundreds of synthetic constructs with predesigned complementarity against endogenous snoRNAs. Our results demonstrate that snoRNA-mediated pseudouridylation can occur on mRNA targets. However, this is typically achieved at relatively low efficiencies, and is constrained by mRNA localization, snoRNA expression levels and the length of the snoRNA:mRNA complementarity stretches. We exploited these insights for the design of snoRNAs targeting pseudouridylation at premature termination codons, which was previously shown to suppress translational termination. However, in this and follow-up experiments in human cells we observe no evidence for significant levels of readthrough of pseudouridylated stop codons. Our study enhances our understanding of the scope, 'design rules', constraints and consequences of snoRNA-mediated pseudouridylation.

bayNorm: Bayesian gene expression recovery, imputation and normalization for single-cell RNA-sequencing data.

MOTIVATION: Normalization of single-cell RNA-sequencing (scRNA-seq) data is a prerequisite to their interpretation. The marked technical variability, high amounts of missing observations and batch effect typical of scRNA-seq datasets make this task particularly challenging. There is a need for an efficient and unified approach for normalization, imputation and batch effect correction. RESULTS: Here, we introduce bayNorm, a novel Bayesian approach for scaling and inference of scRNA-seq counts. The method's likelihood function follows a binomial model of mRNA capture, while priors are estimated from expression values across cells using an empirical Bayes approach. We first validate our assumptions by showing this model can reproduce different statistics observed in real scRNA-seq data. We demonstrate using publicly available scRNA-seq datasets and simulated expression data that bayNorm allows robust imputation of missing values generating realistic transcript distributions that match single molecule fluorescence in situ hybridization measurements. Moreover, by using priors informed by dataset structures, bayNorm improves accuracy and sensitivity of differential expression analysis and reduces batch effect compared with other existing methods. Altogether, bayNorm provides an efficient, integrated solution for global scaling normalization, imputation and true count recovery of gene expression measurements from scRNA-seq data. AVAILABILITY AND IMPLEMENTATION: The R package 'bayNorm' is publishd on bioconductor at https://bioconductor.org/packages/release/bioc/html/bayNorm.html. The code for analyzing data in this article is available at https://github.com/WT215/bayNorm_papercode. SUPPLEMENTARY INFORMATION: Supplementary data are available at Bioinformatics online.

Branch site bulge conformations in domain 6 determine functional sugar puckers in group II intron splicing.

Although group II intron ribozymes are intensively studied the question how structural dynamics affects splicing catalysis has remained elusive. We report for the first time that the group II intron domain 6 exists in a secondary structure equilibrium between a single- and a two-nucleotide bulge conformation, which is directly linked to a switch between sugar puckers of the branch site adenosine. Our study determined a functional sugar pucker equilibrium between the transesterification active C2'-endo conformation of the branch site adenosine in the 1nt bulge and an inactive C3'-endo state in the 2nt bulge fold, allowing the group II intron to switch its activity from the branching to the exon ligation step. Our detailed NMR spectroscopic investigation identified magnesium (II) ions and the branching reaction as regulators of the equilibrium populations. The tuneable secondary structure/sugar pucker equilibrium supports a conformational selection mechanism to up- and downregulate catalytically active and inactive states of the branch site adenosine to orchestrate the multi-step splicing process. The conformational dynamics of group II intron domain 6 is also proposed to be a key aspect for the directionality selection in reversible splicing.

Cell size control driven by the circadian clock and environment in cyanobacteria.

How cells maintain their size has been extensively studied under constant conditions. In the wild, however, cells rarely experience constant environments. Here, we examine how the 24-h circadian clock and environmental cycles modulate cell size control and division timings in the cyanobacterium Synechococcus elongatus using single-cell time-lapse microscopy. Under constant light, wild-type cells follow an apparent sizer-like principle. Closer inspection reveals that the clock generates two subpopulations, with cells born in the subjective day following different division rules from cells born in subjective night. A stochastic model explains how this behavior emerges from the interaction of cell size control with the clock. We demonstrate that the clock continuously modulates the probability of cell division throughout day and night, rather than solely applying an on-off gate to division, as previously proposed. Iterating between modeling and experiments, we go on to identify an effective coupling of the division rate to time of day through the combined effects of the environment and the clock on cell division. Under naturally graded light-dark cycles, this coupling narrows the time window of cell divisions and shifts divisions away from when light levels are low and cell growth is reduced. Our analysis allows us to disentangle, and predict the effects of, the complex interactions between the environment, clock, and cell size control.

Atomic mutagenesis of stop codon nucleotides reveals the chemical prerequisites for release factor-mediated peptide release.

Termination of protein synthesis is triggered by the recognition of a stop codon at the ribosomal A site and is mediated by class I release factors (RFs). Whereas in bacteria, RF1 and RF2 promote termination at UAA/UAG and UAA/UGA stop codons, respectively, eukaryotes only depend on one RF (eRF1) to initiate peptide release at all three stop codons. Based on several structural as well as biochemical studies, interactions between mRNA, tRNA, and rRNA have been proposed to be required for stop codon recognition. In this study, the influence of these interactions was investigated by using chemically modified stop codons. Single functional groups within stop codon nucleotides were substituted to weaken or completely eliminate specific interactions between the respective mRNA and RFs. Our findings provide detailed insight into the recognition mode of bacterial and eukaryotic RFs, thereby revealing the chemical groups of nucleotides that define the identity of stop codons and provide the means to discriminate against noncognate stop codons or UGG sense codons.

Catalytic H2O2 Activation by a Diiron Complex for Methanol Oxidation.

In nature, C-H bond oxidation of CH4 involves a peroxo intermediate that decays to the high-valent active species of either a "closed" {FeIV(µ-O)2FeIV} core or an "open" {FeIV(O)(µ-O)FeIV(O)} core. To mimic and to obtain more mechanistic insight in this reaction mode, we have investigated the reactivity of the bioinspired diiron complex [(susan){Fe(OH)(µ-O)Fe(OH)}]2+ [susan = 4,7-dimethyl-1,1,10,10-tetrakis(2-pyridylmethyl)-1,4,7,10-tetraazadecane], which catalyzes CH3OH oxidation with H2O2 to HCHO and HCO2H. The kinetics is faster in the presence of a proton. 18O-labeling experiments show that the active species, generated by a decay of the initially formed peroxo intermediate [(susan){FeIII(µ-O)(µ-O2)FeIII}]2+, contains one reactive oxygen atom from the µ-oxo and another from the µ-peroxo bridge of its peroxo precursor. Considering an FeIVFeIV active species, a "closed" {FeIV(µ-O)2FeIV} core explains the observed labeling results, while a scrambling of the terminal and bridging oxo ligands is required to account for an "open" {FeIV(O)(µ-O)FeIV(O)} core.

Bounding the stationary distributions of the chemical master equation via mathematical programming.

The stochastic dynamics of biochemical networks are usually modeled with the chemical master equation (CME). The stationary distributions of CMEs are seldom solvable analytically, and numerical methods typically produce estimates with uncontrolled errors. Here, we introduce mathematical programming approaches that yield approximations of these distributions with computable error bounds which enable the verification of their accuracy. First, we use semidefinite programming to compute increasingly tighter upper and lower bounds on the moments of the stationary distributions for networks with rational propensities. Second, we use these moment bounds to formulate linear programs that yield convergent upper and lower bounds on the stationary distributions themselves, their marginals, and stationary averages. The bounds obtained also provide a computational test for the uniqueness of the distribution. In the unique case, the bounds form an approximation of the stationary distribution with a computable bound on its error. In the nonunique case, our approach yields converging approximations of the ergodic distributions. We illustrate our methodology through several biochemical examples taken from the literature: Schlögl's model for a chemical bifurcation, a two-dimensional toggle switch, a model for bursty gene expression, and a dimerization model with multiple stationary distributions.

Reversible Carboxylate Shift in a µ-Oxo Diferric Complex in Solution by Acid-/Base-Addition.

A reversible carboxylate shift has been observed in a µ-oxo diferric complex in solution by UV-vis-NIR and FTIR spectroscopy triggered by the addition of a base or an acid. A terminal acetate decoordinates upon the addition of a proton, resulting in a shift of the remaining terminal acetato to a µ-η1:η1 bridge. The addition of a base restores the original structure containing only terminal acetates. The implications for metalloenzymes with carboxylate-bridged nonheme diiron active sites are discussed.

Two Unsupported Terminal Hydroxido Ligands in a µ-Oxo-Bridged Ferric Dimer: Protonation and Kinetic Lability Studies.

The dinuclear complex [(susan){FeIII(OH)(µ-O)FeIII(OH)}](ClO4)2 (Fe2(OH)2(ClO4)2; susan = 4,7-dimethyl-1,1,10,10-tetra(2-pyridylmethyl)-1,4,7,10-tetraazadecane) with two unsupported terminal hydroxido ligands and for comparison the fluorido-substituted complex [(susan){FeIIIF(µ-O)FeIIIF}](ClO4)2 (Fe2F2(ClO4)2) have been synthesized and characterized in the solid state as well in acetonitrile (CH3CN) and water (H2O) solutions. The Fe-OH bonds are strongly modulated by intermolecular hydrogen bonds (1.85 and 1.90 Å). UV-vis-near-IR (NIR) and Mössbauer spectroscopies prove that Fe2F22+ and Fe2(OH)22+ retain their structural integrity in a CH3CN solution. The OH- ligand induces a weaker ligand field than the F- ligand because of stronger π donation. This increased electron donation shifts the potential for the irreversible oxidation by 610 mV cathodically from 1.40 V in Fe2F22+ to 0.79 V versus Fc+/Fc in Fe2(OH)22+. Protonation/deprotonation studies in CH3CN and aqueous solutions of Fe2(OH)22+ provide two reversible acid-base equilibria. UV-vis-NIR, Mössbauer, and cryo electrospray ionization mass spectrometry experiments show conservation of the mono(µ-oxo) bridging motif, while the terminal OH- ligands are protonated to H2O. Titration experiments in aqueous solution at room temperature provide the p Ka values as p K1 = 4.9 and p K2 = 6.8. Kinetic studies by temperature- and pressure-dependent 17O NMR spectrometry revealed for the first time the water-exchange parameters [ kex298 = (3.9 ± 0.2) × 105 s-1, Δ H⧧ = 39.6 ± 0.2 kJ mol-1, Δ S⧧ = -5.1 ± 1 J mol-1 K-1, and Δ V⧧ = +3.0 ± 0.2 cm3 mol-1] and the underlying Id mechanism for a {FeIII(OH2)(µ-O)FeIII(OH2)} core. The same studies suggest that in solution the monoprotonated {FeIII(OH)(µ-O)FeIII(OH2)} complex has µ-O and µ-O2H3 bridges between the two Fe centers.

Nucleotide modifications within bacterial messenger RNAs regulate their translation and are able to rewire the genetic code.

Nucleotide modifications within RNA transcripts are found in every organism in all three domains of life. 6-methyladeonsine (m(6)A), 5-methylcytosine (m(5)C) and pseudouridine (Ψ) are highly abundant nucleotide modifications in coding sequences of eukaryal mRNAs, while m(5)C and m(6)A modifications have also been discovered in archaeal and bacterial mRNAs. Employing in vitro translation assays, we systematically investigated the influence of nucleotide modifications on translation. We introduced m(5)C, m(6)A, Ψ or 2'-O-methylated nucleotides at each of the three positions within a codon of the bacterial ErmCL mRNA and analyzed their influence on translation. Depending on the respective nucleotide modification, as well as its position within a codon, protein synthesis remained either unaffected or was prematurely terminated at the modification site, resulting in reduced amounts of the full-length peptide. In the latter case, toeprint analysis of ribosomal complexes was consistent with stalling of translation at the modified codon. When multiple nucleotide modifications were introduced within one codon, an additive inhibitory effect on translation was observed. We also identified the m(5)C modification to alter the amino acid identity of the corresponding codon, when positioned at the second codon position. Our results suggest a novel mode of gene regulation by nucleotide modifications in bacterial mRNAs.

Stochastic Simulation of Biomolecular Networks in Dynamic Environments.

Simulation of biomolecular networks is now indispensable for studying biological systems, from small reaction networks to large ensembles of cells. Here we present a novel approach for stochastic simulation of networks embedded in the dynamic environment of the cell and its surroundings. We thus sample trajectories of the stochastic process described by the chemical master equation with time-varying propensities. A comparative analysis shows that existing approaches can either fail dramatically, or else can impose impractical computational burdens due to numerical integration of reaction propensities, especially when cell ensembles are studied. Here we introduce the Extrande method which, given a simulated time course of dynamic network inputs, provides a conditionally exact and several orders-of-magnitude faster simulation solution. The new approach makes it feasible to demonstrate-using decision-making by a large population of quorum sensing bacteria-that robustness to fluctuations from upstream signaling places strong constraints on the design of networks determining cell fate. Our approach has the potential to significantly advance both understanding of molecular systems biology and design of synthetic circuits.

Inference for Stochastic Chemical Kinetics Using Moment Equations and System Size Expansion.

Quantitative mechanistic models are valuable tools for disentangling biochemical pathways and for achieving a comprehensive understanding of biological systems. However, to be quantitative the parameters of these models have to be estimated from experimental data. In the presence of significant stochastic fluctuations this is a challenging task as stochastic simulations are usually too time-consuming and a macroscopic description using reaction rate equations (RREs) is no longer accurate. In this manuscript, we therefore consider moment-closure approximation (MA) and the system size expansion (SSE), which approximate the statistical moments of stochastic processes and tend to be more precise than macroscopic descriptions. We introduce gradient-based parameter optimization methods and uncertainty analysis methods for MA and SSE. Efficiency and reliability of the methods are assessed using simulation examples as well as by an application to data for Epo-induced JAK/STAT signaling. The application revealed that even if merely population-average data are available, MA and SSE improve parameter identifiability in comparison to RRE. Furthermore, the simulation examples revealed that the resulting estimates are more reliable for an intermediate volume regime. In this regime the estimation error is reduced and we propose methods to determine the regime boundaries. These results illustrate that inference using MA and SSE is feasible and possesses a high sensitivity.

A Mixed-Valence Fluorido-Bridged FeIIFeIII Complex.

The reaction of the new dinucleating ligand susan6-Me with Fe(BF4)2·6H2O results in formation of the homovalent FeIIFeII complex [(susan6-Me){FeII(µ-F)2FeII}]2+ and the mixed-valence FeIIFeIII complex [(susan6-Me){FeIIF(µ-F)FeIIIF}]2+ depending on the absence or presence of dioxygen, respectively. Complex [(susan6-Me){FeIIF(µ-F)FeIIIF}]2+ is the first molecular mixed-valence complex with a fluorido bridge. The short FeIII-µ-F bond of 1.87 Å causes a large reorganization energy, resulting in a localized class II system with an intervalence charge-transfer band of high energy at 10000 cm-1.

Phenotypic switching in gene regulatory networks.

Noise in gene expression can lead to reversible phenotypic switching. Several experimental studies have shown that the abundance distributions of proteins in a population of isogenic cells may display multiple distinct maxima. Each of these maxima may be associated with a subpopulation of a particular phenotype, the quantification of which is important for understanding cellular decision-making. Here, we devise a methodology which allows us to quantify multimodal gene expression distributions and single-cell power spectra in gene regulatory networks. Extending the commonly used linear noise approximation, we rigorously show that, in the limit of slow promoter dynamics, these distributions can be systematically approximated as a mixture of Gaussian components in a wide class of networks. The resulting closed-form approximation provides a practical tool for studying complex nonlinear gene regulatory networks that have thus far been amenable only to stochastic simulation. We demonstrate the applicability of our approach in a number of genetic networks, uncovering previously unidentified dynamical characteristics associated with phenotypic switching. Specifically, we elucidate how the interplay of transcriptional and translational regulation can be exploited to control the multimodality of gene expression distributions in two-promoter networks. We demonstrate how phenotypic switching leads to birhythmical expression in a genetic oscillator, and to hysteresis in phenotypic induction, thus highlighting the ability of regulatory networks to retain memory.

mRNA modifications: Dynamic regulators of gene expression?

The expression of a gene is a tightly regulated process and is exerted by a myriad of different mechanisms. Recently, RNA modifications located in coding sequences of mRNAs, have been identified as potential regulators of gene expression. N(6)-methyladenosine (m(6)A), 5-methylcytosine (m(5)C), pseudouridine (Ψ) and N(1)-methyladenosine (m(1)A) have been found within open reading frames of mRNAs. The presence of these mRNA modifications has been implicated to modulate the fate of an mRNA, ranging from maturation to its translation and even degradation. However, many aspects concerning the biological functions of mRNA modifications remain elusive. Recently, systematic in vitro studies allowed a first glimpse of the direct interplay of mRNA modifications and the efficiency and fidelity of ribosomal translation. It thereby became evident that the effects of mRNA modifications were, astonishingly versatile, depending on the type, position or sequence context. The incorporation of a single modification could either prematurely terminate protein synthesis, reduce the peptide yield or alter the amino acid sequence identity. These results implicate that mRNA modifications are a powerful mechanism to post-transcriptionally regulate gene expression.

Feedback between stochastic gene networks and population dynamics enables cellular decision-making.

Phenotypic selection occurs when genetically identical cells are subject to different reproductive abilities due to cellular noise. Such noise arises from fluctuations in reactions synthesizing proteins and plays a crucial role in how cells make decisions and respond to stress or drugs. We propose a general stochastic agent-based model for growing populations capturing the feedback between gene expression and cell division dynamics. We devise a finite state projection approach to analyze gene expression and division distributions and infer selection from single-cell data in mother machines and lineage trees. We use the theory to quantify selection in multi-stable gene expression networks and elucidate that the trade-off between phenotypic switching and selection enables robust decision-making essential for synthetic circuits and developmental lineage decisions. Using live-cell data, we demonstrate that combining theory and inference provides quantitative insights into bet-hedging-like response to DNA damage and adaptation during antibiotic exposure in Escherichia coli.

Global transcription regulation revealed from dynamical correlations in time-resolved single-cell RNA sequencing.

Single-cell transcriptomics reveals significant variations in transcriptional activity across cells. Yet, it remains challenging to identify mechanisms of transcription dynamics from static snapshots. It is thus still unknown what drives global transcription dynamics in single cells. We present a stochastic model of gene expression with cell size- and cell cycle-dependent rates in growing and dividing cells that harnesses temporal dimensions of single-cell RNA sequencing through metabolic labeling protocols and cel lcycle reporters. We develop a parallel and highly scalable approximate Bayesian computation method that corrects for technical variation and accurately quantifies absolute burst frequency, burst size, and degradation rate along the cell cycle at a transcriptome-wide scale. Using Bayesian model selection, we reveal scaling between transcription rates and cell size and unveil waves of gene regulation across the cell cycle-dependent transcriptome. Our study shows that stochastic modeling of dynamical correlations identifies global mechanisms of transcription regulation. A record of this paper's transparent peer review process is included in the supplemental information.

How reliable is the linear noise approximation of gene regulatory networks?

BACKGROUND: The linear noise approximation (LNA) is commonly used to predict how noise is regulated and exploited at the cellular level. These predictions are exact for reaction networks composed exclusively of first order reactions or for networks involving bimolecular reactions and large numbers of molecules. It is however well known that gene regulation involves bimolecular interactions with molecule numbers as small as a single copy of a particular gene. It is therefore questionable how reliable are the LNA predictions for these systems. RESULTS: We implement in the software package intrinsic Noise Analyzer (iNA), a system size expansion based method which calculates the mean concentrations and the variances of the fluctuations to an order of accuracy higher than the LNA. We then use iNA to explore the parametric dependence of the Fano factors and of the coefficients of variation of the mRNA and protein fluctuations in models of genetic networks involving nonlinear protein degradation, post-transcriptional, post-translational and negative feedback regulation. We find that the LNA can significantly underestimate the amplitude and period of noise-induced oscillations in genetic oscillators. We also identify cases where the LNA predicts that noise levels can be optimized by tuning a bimolecular rate constant whereas our method shows that no such regulation is possible. All our results are confirmed by stochastic simulations. CONCLUSION: The software iNA allows the investigation of parameter regimes where the LNA fares well and where it does not. We have shown that the parametric dependence of the coefficients of variation and Fano factors for common gene regulatory networks is better described by including terms of higher order than LNA in the system size expansion. This analysis is considerably faster than stochastic simulations due to the extensive ensemble averaging needed to obtain statistically meaningful results. Hence iNA is well suited for performing computationally efficient and quantitative studies of intrinsic noise in gene regulatory networks.

Signatures of nonlinearity in single cell noise-induced oscillations.

A class of theoretical models seeks to explain rhythmic single cell data by postulating that they are generated by intrinsic noise in biochemical systems whose deterministic models exhibit only damped oscillations. The main features of such noise-induced oscillations are quantified by the power spectrum which measures the dependence of the oscillatory signal's power with frequency. In this paper we derive an approximate closed-form expression for the power spectrum of any monostable biochemical system close to a Hopf bifurcation, where noise-induced oscillations are most pronounced. Unlike the commonly used linear noise approximation which is valid in the macroscopic limit of large volumes, our theory is valid over a wide range of volumes and hence affords a more suitable description of single cell noise-induced oscillations. Our theory predicts that the spectra have three universal features: (i) a dominant peak at some frequency, (ii) a smaller peak at twice the frequency of the dominant peak and (iii) a peak at zero frequency. Of these, the linear noise approximation predicts only the first feature while the remaining two stem from the combination of intrinsic noise and nonlinearity in the law of mass action. The theoretical expressions are shown to accurately match the power spectra determined from stochastic simulations of mitotic and circadian oscillators. Furthermore it is shown how recently acquired single cell rhythmic fibroblast data displays all the features predicted by our theory and that the experimental spectrum is well described by our theory but not by the conventional linear noise approximation.