Identifying Subpopulations of Cells in Single-Cell Transcriptomic Data: A Bayesian Mixture Modeling Approach to Zero Inflation of Counts.

In the study of single-cell RNA-seq (scRNA-Seq) data, a key component of the analysis is to identify subpopulations of cells in the data. A variety of approaches to this have been considered, and although many machine learning-based methods have been developed, these rarely give an estimate of uncertainty in the cluster assignment. To allow for this, probabilistic models have been developed, but scRNA-Seq data exhibit a phenomenon known as dropout, whereby a large proportion of the observed read counts are zero. This poses challenges in developing probabilistic models that appropriately model the data. We develop a novel Dirichlet process mixture model that employs both a mixture at the cell level to model multiple populations of cells and a zero-inflated negative binomial mixture of counts at the transcript level. By taking a Bayesian approach, we are able to model the expression of genes within clusters, and to quantify uncertainty in cluster assignments. It is shown that this approach outperforms previous approaches that applied multinomial distributions to model scRNA-Seq counts and negative binomial models that do not take into account zero inflation. Applied to a publicly available data set of scRNA-Seq counts of multiple cell types from the mouse cortex and hippocampus, we demonstrate how our approach can be used to distinguish subpopulations of cells as clusters in the data, and to identify gene sets that are indicative of membership of a subpopulation.

Inhibition of the angiotensin II type 2 receptor AT2R is a novel therapeutic strategy for glioblastoma.

Glioblastoma (GBM) is an aggressive malignant primary brain tumor with limited therapeutic options. We show that the angiotensin II (AngII) type 2 receptor (AT2R) is a therapeutic target for GBM and that AngII, endogenously produced in GBM cells, promotes proliferation through AT2R. We repurposed EMA401, an AT2R antagonist originally developed as a peripherally restricted analgesic, for GBM and showed that it inhibits the proliferation of AT2R-expressing GBM spheroids and blocks their invasiveness and angiogenic capacity. The crystal structure of AT2R bound to EMA401 was determined and revealed the receptor to be in an active-like conformation with helix-VIII blocking G-protein or ß-arrestin recruitment. The architecture and interactions of EMA401 in AT2R differ drastically from complexes of AT2R with other relevant compounds. To enhance central nervous system (CNS) penetration of EMA401, we exploited the crystal structure to design an angiopep-2-tethered EMA401 derivative, A3E. A3E exhibited enhanced CNS penetration, leading to reduced tumor volume, inhibition of proliferation, and increased levels of apoptosis in an orthotopic xenograft model of GBM.

Can trainees safely perform pancreatoenteric anastomosis? A systematic review, meta-analysis, and risk-adjusted analysis of postoperative pancreatic fistula.

BACKGROUND: The complexity of pancreaticoduodenectomy and fear of morbidity, particularly postoperative pancreatic fistula, can be a barrier to surgical trainees gaining operative experience. This meta-analysis sought to compare the postoperative pancreatic fistula rate after pancreatoenteric anastomosis by trainees or established surgeons. METHODS: A systematic review of the literature was performed using Preferred Reporting Items for Systematic Reviews and Meta-Analyses guidelines, with differences in postoperative pancreatic fistula rates after pancreatoenteric anastomosis between trainee-led versus consultant/attending surgeons pooled using meta-analysis. Variation in rates of postoperative pancreatic fistula was further explored using risk-adjusted outcomes using published risk scores and cumulative sum control chart analysis in a retrospective cohort. RESULTS: Across 14 cohorts included in the meta-analysis, trainees tended toward a lower but nonsignificant rate of all postoperative pancreatic fistula (odds ratio: 0.77, P = .45) and clinically relevant postoperative pancreatic fistula (odds ratio: 0.69, P = .37). However, there was evidence of case selection, with trainees being less likely to operate on patients with a pancreatic duct width <3 mm (odds ratio: 0.45, P = .05). Similarly, analysis of a retrospective cohort (N = 756 cases) found patients operated by trainees to have significantly lower predicted all postoperative pancreatic fistula (median: 20 vs 26%, P < .001) and clinically relevant postoperative pancreatic fistula (7 vs 9%, P = .020) rates than consultant/attending surgeons, based on preoperative risk scores. After adjusting for this on multivariable analysis, the risks of all postoperative pancreatic fistula (odds ratio: 1.18, P = .604) and clinically relevant postoperative pancreatic fistula (odds ratio: 0.85, P = .693) remained similar after pancreatoenteric anastomosis by trainees or consultant/attending surgeons. CONCLUSION: Pancreatoenteric anastomosis, when performed by trainees, is associated with acceptable outcomes. There is evidence of case selection among patients undergoing surgery by trainees; hence, risk adjustment provides a critical tool for the objective evaluation of performance.

Topological approximate Bayesian computation for parameter inference of an angiogenesis model.

MOTIVATION: Inferring the parameters of models describing biological systems is an important problem in the reverse engineering of the mechanisms underlying these systems. Much work has focused on parameter inference of stochastic and ordinary differential equation models using Approximate Bayesian Computation (ABC). While there is some recent work on inference in spatial models, this remains an open problem. Simultaneously, advances in topological data analysis (TDA), a field of computational mathematics, have enabled spatial patterns in data to be characterized. RESULTS: Here, we focus on recent work using TDA to study different regimes of parameter space for a well-studied model of angiogenesis. We propose a method for combining TDA with ABC to infer parameters in the Anderson-Chaplain model of angiogenesis. We demonstrate that this topological approach outperforms ABC approaches that use simpler statistics based on spatial features of the data. This is a first step toward a general framework of spatial parameter inference for biological systems, for which there may be a variety of filtrations, vectorizations and summary statistics to be considered. AVAILABILITY AND IMPLEMENTATION: All code used to produce our results is available as a Snakemake workflow from github.com/tt104/tabc_angio.

External validation of postoperative pancreatic fistula prediction scores in pancreatoduodenectomy: a systematic review and meta-analysis.

BACKGROUND: Multiple risk scores claim to predict the probability of postoperative pancreatic fistula (POPF) after pancreatoduodenectomy. It is unclear which scores have undergone external validation and are the most accurate. The aim of this study was to identify risk scores for POPF, and assess the clinical validity of these scores. METHODS: Areas under receiving operator characteristic curve (AUROCs) were extracted from studies that performed external validation of POPF risk scores. These were pooled for each risk score, using intercept-only random-effects meta-regression models. RESULTS: Systematic review identified 34 risk scores, of which six had been subjected to external validation, and so included in the meta-analysis, (Tokyo (N=2 validation studies), Birmingham (N=5), FRS (N=19), a-FRS (N=12), m-FRS (N=3) and ua-FRS (N=3) scores). Overall predictive accuracies were similar for all six scores, with pooled AUROCs of 0.61, 0.70, 0.71, 0.70, 0.70 and 0.72, respectively. Considerably heterogeneity was observed, with I2 statistics ranging from 52.1-88.6%. CONCLUSION: Most risk scores lack external validation; where this was performed, risk scores were found to have limited predictive accuracy. . Consensus is needed for which score to use in clinical practice. Due to the limited predictive accuracy, future studies to derive a more accurate risk score are warranted.

Approximate inference of gene regulatory network models from RNA-Seq time series data.

BACKGROUND: Inference of gene regulatory network structures from RNA-Seq data is challenging due to the nature of the data, as measurements take the form of counts of reads mapped to a given gene. Here we present a model for RNA-Seq time series data that applies a negative binomial distribution for the observations, and uses sparse regression with a horseshoe prior to learn a dynamic Bayesian network of interactions between genes. We use a variational inference scheme to learn approximate posterior distributions for the model parameters. RESULTS: The methodology is benchmarked on synthetic data designed to replicate the distribution of real world RNA-Seq data. We compare our method to other sparse regression approaches and find improved performance in learning directed networks. We demonstrate an application of our method to a publicly available human neuronal stem cell differentiation RNA-Seq time series data set to infer the underlying network structure. CONCLUSIONS: Our method is able to improve performance on synthetic data by explicitly modelling the statistical distribution of the data when learning networks from RNA-Seq time series. Applying approximate inference techniques we can learn network structures quickly with only moderate computing resources.

NetDiff - Bayesian model selection for differential gene regulatory network inference.

Differential networks allow us to better understand the changes in cellular processes that are exhibited in conditions of interest, identifying variations in gene regulation or protein interaction between, for example, cases and controls, or in response to external stimuli. Here we present a novel methodology for the inference of differential gene regulatory networks from gene expression microarray data. Specifically we apply a Bayesian model selection approach to compare models of conserved and varying network structure, and use Gaussian graphical models to represent the network structures. We apply a variational inference approach to the learning of Gaussian graphical models of gene regulatory networks, that enables us to perform Bayesian model selection that is significantly more computationally efficient than Markov Chain Monte Carlo approaches. Our method is demonstrated to be more robust than independent analysis of data from multiple conditions when applied to synthetic network data, generating fewer false positive predictions of differential edges. We demonstrate the utility of our approach on real world gene expression microarray data by applying it to existing data from amyotrophic lateral sclerosis cases with and without mutations in C9orf72, and controls, where we are able to identify differential network interactions for further investigation.

Empirical likelihood tests for nonparametric detection of differential expression from RNA-seq data.

The availability of large quantities of transcriptomic data in the form of RNA-seq count data has necessitated the development of methods to identify genes differentially expressed between experimental conditions. Many existing approaches apply a parametric model of gene expression and so place strong assumptions on the distribution of the data. Here we explore an alternate nonparametric approach that applies an empirical likelihood framework, allowing us to define likelihoods without specifying a parametric model of the data. We demonstrate the performance of our method when applied to gold standard datasets, and to existing experimental data. Our approach outperforms or closely matches performance of existing methods in the literature, and requires modest computational resources. An R package, EmpDiff implementing the methods described in the paper is available from: http://homepages.inf.ed.ac.uk/tthorne/software/packages/EmpDiff_0.99.tar.gz.

Derivative processes for modelling metabolic fluxes.

MOTIVATION: One of the challenging questions in modelling biological systems is to characterize the functional forms of the processes that control and orchestrate molecular and cellular phenotypes. Recently proposed methods for the analysis of metabolic pathways, for example, dynamic flux estimation, can only provide estimates of the underlying fluxes at discrete time points but fail to capture the complete temporal behaviour. To describe the dynamic variation of the fluxes, we additionally require the assumption of specific functional forms that can capture the temporal behaviour. However, it also remains unclear how to address the noise which might be present in experimentally measured metabolite concentrations. RESULTS: Here we propose a novel approach to modelling metabolic fluxes: derivative processes that are based on multiple-output Gaussian processes (MGPs), which are a flexible non-parametric Bayesian modelling technique. The main advantages that follow from MGPs approach include the natural non-parametric representation of the fluxes and ability to impute the missing data in between the measurements. Our derivative process approach allows us to model changes in metabolite derivative concentrations and to characterize the temporal behaviour of metabolic fluxes from time course data. Because the derivative of a Gaussian process is itself a Gaussian process, we can readily link metabolite concentrations to metabolic fluxes and vice versa. Here we discuss how this can be implemented in an MGP framework and illustrate its application to simple models, including nitrogen metabolism in Escherichia coli. AVAILABILITY AND IMPLEMENTATION: R code is available from the authors upon request.

Model selection in systems and synthetic biology.

Developing mechanistic models has become an integral aspect of systems biology, as has the need to differentiate between alternative models. Parameterizing mathematical models has been widely perceived as a formidable challenge, which has spurred the development of statistical and optimisation routines for parameter inference. But now focus is increasingly shifting to problems that require us to choose from among a set of different models to determine which one offers the best description of a given biological system. We will here provide an overview of recent developments in the area of model selection. We will focus on approaches that are both practical as well as build on solid statistical principles and outline the conceptual foundations and the scope for application of such methods in systems biology.

Graphical modelling of molecular networks underlying sporadic inclusion body myositis.

Here we present a novel statistical methodology that allows us to analyze gene expression data that have been collected from a number of different cases or conditions in a unified framework. Using a Bayesian nonparametric framework we develop a hierarchical model wherein genes can maintain a shared set of interactions between different cases, whilst also exhibiting behaviour that is unique to specific cases, sets of conditions, or groups of data points. By doing so we are able to not only combine data from different cases but also to discern the unique regulatory interactions that differentiate the cases. We apply our method to clinical data collected from patients suffering from sporadic Inclusion Body Myositis (sIBM), as well as control samples, and demonstrate the ability of our method to infer regulatory interactions that are unique to the disease cases of interest. The method thus balances the statistical need to include as many patients and controls as possible, and the clinical need to maintain potentially cryptic differences among patients and between patients and controls at the regulatory level.

Inference of temporally varying Bayesian networks.

MOTIVATION: When analysing gene expression time series data, an often overlooked but crucial aspect of the model is that the regulatory network structure may change over time. Although some approaches have addressed this problem previously in the literature, many are not well suited to the sequential nature of the data. RESULTS: Here, we present a method that allows us to infer regulatory network structures that may vary between time points, using a set of hidden states that describe the network structure at a given time point. To model the distribution of the hidden states, we have applied the Hierarchical Dirichlet Process Hidden Markov Model, a non-parametric extension of the traditional Hidden Markov Model, which does not require us to fix the number of hidden states in advance. We apply our method to existing microarray expression data as well as demonstrating is efficacy on simulated test data.

Parameter-free model discrimination criterion based on steady-state coplanarity.

We introduce a procedure for deciding when a mass-action model is incompatible with observed steady-state data that does not require any parameter estimation. Thus, we avoid the difficulties of nonlinear optimization typically associated with methods based on parameter fitting. Instead, we borrow ideas from algebraic geometry to construct a transformation of the model variables such that any set of steady states of the model under that transformation lies on a common plane, irrespective of the values of the model parameters. Model rejection can then be performed by assessing the degree to which the transformed data deviate from coplanarity. We demonstrate our method by applying it to models of multisite phosphorylation and cell death signaling. Our framework offers a parameter-free perspective on the statistical model selection problem, which can complement conventional statistical methods in certain classes of problems where inference has to be based on steady-state data and the model structures allow for suitable algebraic relationships among the steady-state solutions.

A systems biology analysis of long and short-term memories of osmotic stress adaptation in fungi.

BACKGROUND: Saccharomyces cerevisiae senses hyperosmotic conditions via the HOG signaling network that activates the stress-activated protein kinase, Hog1, and modulates metabolic fluxes and gene expression to generate appropriate adaptive responses. The integral control mechanism by which Hog1 modulates glycerol production remains uncharacterized. An additional Hog1-independent mechanism retains intracellular glycerol for adaptation. Candida albicans also adapts to hyperosmolarity via a HOG signaling network. However, it remains unknown whether Hog1 exerts integral or proportional control over glycerol production in C. albicans. RESULTS: We combined modeling and experimental approaches to study osmotic stress responses in S. cerevisiae and C. albicans. We propose a simple ordinary differential equation (ODE) model that highlights the integral control that Hog1 exerts over glycerol biosynthesis in these species. If integral control arises from a separation of time scales (i.e. rapid HOG activation of glycerol production capacity which decays slowly under hyperosmotic conditions), then the model predicts that glycerol production rates elevate upon adaptation to a first stress and this makes the cell adapts faster to a second hyperosmotic stress. It appears as if the cell is able to remember the stress history that is longer than the timescale of signal transduction. This is termed the long-term stress memory. Our experimental data verify this. Like S. cerevisiae, C. albicans mimimizes glycerol efflux during adaptation to hyperosmolarity. Also, transient activation of intermediate kinases in the HOG pathway results in a short-term memory in the signaling pathway. This determines the amplitude of Hog1 phosphorylation under a periodic sequence of stress and non-stressed intervals. Our model suggests that the long-term memory also affects the way a cell responds to periodic stress conditions. Hence, during osmohomeostasis, short-term memory is dependent upon long-term memory. This is relevant in the context of fungal responses to dynamic and changing environments. CONCLUSIONS: Our experiments and modeling have provided an example of identifying integral control that arises from time-scale separation in different processes, which is an important functional module in various contexts.

Graph spectral analysis of protein interaction network evolution.

We present an analysis of protein interaction network data via the comparison of models of network evolution to the observed data. We take a bayesian approach and perform posterior density estimation using an approximate bayesian computation with sequential Monte Carlo method. Our approach allows us to perform model selection over a selection of potential network growth models. The methodology we apply uses a distance defined in terms of graph spectra which captures the network data more naturally than previously used summary statistics such as the degree distribution. Furthermore, we include the effects of sampling into the analysis, to properly correct for the incompleteness of existing datasets, and have analysed the performance of our method under various degrees of sampling. We consider a number of models focusing not only on the biologically relevant class of duplication models, but also including models of scale-free network growth that have previously been claimed to describe such data. We find a preference for a duplication-divergence with linear preferential attachment model in the majority of the interaction datasets considered. We also illustrate how our method can be used to perform multi-model inference of network parameters to estimate properties of the full network from sampled data.

Combinatorial stresses kill pathogenic Candida species.

Pathogenic microbes exist in dynamic niches and have evolved robust adaptive responses to promote survival in their hosts. The major fungal pathogens of humans, Candida albicans and Candida glabrata, are exposed to a range of environmental stresses in their hosts including osmotic, oxidative and nitrosative stresses. Significant efforts have been devoted to the characterization of the adaptive responses to each of these stresses. In the wild, cells are frequently exposed simultaneously to combinations of these stresses and yet the effects of such combinatorial stresses have not been explored. We have developed a common experimental platform to facilitate the comparison of combinatorial stress responses in C. glabrata and C. albicans. This platform is based on the growth of cells in buffered rich medium at 30°C, and was used to define relatively low, medium and high doses of osmotic (NaCl), oxidative (H(2)O(2)) and nitrosative stresses (e.g., dipropylenetriamine (DPTA)-NONOate). The effects of combinatorial stresses were compared with the corresponding individual stresses under these growth conditions. We show for the first time that certain combinations of combinatorial stress are especially potent in terms of their ability to kill C. albicans and C. glabrata and/or inhibit their growth. This was the case for combinations of osmotic plus oxidative stress and for oxidative plus nitrosative stress. We predict that combinatorial stresses may be highly significant in host defences against these pathogenic yeasts.

Calibrating spatio-temporal models of leukocyte dynamics against in vivo live-imaging data using approximate Bayesian computation.

In vivo studies allow us to investigate biological processes at the level of the organism. But not all aspects of in vivo systems are amenable to direct experimental measurements. In order to make the most of such data we therefore require statistical tools that allow us to obtain reliable estimates for e.g. kinetic in vivo parameters. Here we show how we can use approximate Bayesian computation approaches in order to analyse leukocyte migration in zebrafish embryos in response to injuries. We track individual leukocytes using live imaging following surgical injury to the embryos' tail-fins. The signalling gradient that leukocytes follow towards the site of the injury cannot be directly measured but we can estimate its shape and how it changes with time from the directly observed patterns of leukocyte migration. By coupling simple models of immune signalling and leukocyte migration with the unknown gradient shape into a single statistical framework we can gain detailed insights into the tissue-wide processes that are involved in the innate immune response to wound injury. In particular we find conclusive evidence for a temporally and spatially changing signalling gradient that modulates the changing activity of the leukocyte population in the embryos. We conclude with a robustness analysis which highlights the most important factors determining the leukocyte dynamics. Our approach relies only on the ability to simulate numerically the process under investigation and is therefore also applicable in other in vivo contexts and studies.

Prediction of putative protein interactions through evolutionary analysis of osmotic stress response in the model yeast Saccharomyces cerevisae.

The osmotic stress response signalling pathway of the model yeast Saccharomyces cerevisae is crucial for the survival of cells under osmotic stress, and is preserved to varying degrees in other related fungal species. We apply a method for inference of ancestral states of characteristics over a phylogeny to 17 fungal species to infer the maximum likelihood estimate of presence or absence in ancestral genomes of genes involved in osmotic stress response. The same method allows us furthermore to perform a statistical test for correlated evolution between genes. Where such correlations exist within the osmotic stress response pathway of S. cerevisae, we have used this in order to predict and subsequently test for the presence of physical protein-protein interactions in an attempt to detect novel interactions. Finally we assess the relevance of observed evolutionary correlations in predicting protein interactions in light of the experimental results. We do find that correlated evolution provides some useful information for the prediction of protein-protein interactions, but that these alone are not sufficient to explain detectable patterns of correlated evolution.

Estimating the size of the human interactome.

After the completion of the human and other genome projects it emerged that the number of genes in organisms as diverse as fruit flies, nematodes, and humans does not reflect our perception of their relative complexity. Here, we provide reliable evidence that the size of protein interaction networks in different organisms appears to correlate much better with their apparent biological complexity. We develop a stable and powerful, yet simple, statistical procedure to estimate the size of the whole network from subnet data. This approach is then applied to a range of eukaryotic organisms for which extensive protein interaction data have been collected and we estimate the number of interactions in humans to be approximately 650,000. We find that the human interaction network is one order of magnitude bigger than the Drosophila melanogaster interactome and approximately 3 times bigger than in Caenorhabditis elegans.

Generating confidence intervals on biological networks.

BACKGROUND: In the analysis of networks we frequently require the statistical significance of some network statistic, such as measures of similarity for the properties of interacting nodes. The structure of the network may introduce dependencies among the nodes and it will in general be necessary to account for these dependencies in the statistical analysis. To this end we require some form of Null model of the network: generally rewired replicates of the network are generated which preserve only the degree (number of interactions) of each node. We show that this can fail to capture important features of network structure, and may result in unrealistic significance levels, when potentially confounding additional information is available. METHODS: We present a new network resampling Null model which takes into account the degree sequence as well as available biological annotations. Using gene ontology information as an illustration we show how this information can be accounted for in the resampling approach, and the impact such information has on the assessment of statistical significance of correlations and motif-abundances in the Saccharomyces cerevisiae protein interaction network. An algorithm, GOcardShuffle, is introduced to allow for the efficient construction of an improved Null model for network data. RESULTS: We use the protein interaction network of S. cerevisiae; correlations between the evolutionary rates and expression levels of interacting proteins and their statistical significance were assessed for Null models which condition on different aspects of the available data. The novel GOcardShuffle approach results in a Null model for annotated network data which appears better to describe the properties of real biological networks. CONCLUSION: An improved statistical approach for the statistical analysis of biological network data, which conditions on the available biological information, leads to qualitatively different results compared to approaches which ignore such annotations. In particular we demonstrate the effects of the biological organization of the network can be sufficient to explain the observed similarity of interacting proteins.