Julia for biologists.

Major computational challenges exist in relation to the collection, curation, processing and analysis of large genomic and imaging datasets, as well as the simulation of larger and more realistic models in systems biology. Here we discuss how a relative newcomer among programming languages-Julia-is poised to meet the current and emerging demands in the computational biosciences and beyond. Speed, flexibility, a thriving package ecosystem and readability are major factors that make high-performance computing and data analysis available to an unprecedented degree. We highlight how Julia's design is already enabling new ways of analyzing biological data and systems, and we provide a list of resources that can facilitate the transition into Julian computing.

A single-cell resolved cell-cell communication model explains lineage commitment in hematopoiesis.

Cells do not make fate decisions independently. Arguably, every cell-fate decision occurs in response to environmental signals. In many cases, cell-cell communication alters the dynamics of the internal gene regulatory network of a cell to initiate cell-fate transitions, yet models rarely take this into account. Here, we have developed a multiscale perspective to study the granulocyte-monocyte versus megakaryocyte-erythrocyte fate decisions. This transition is dictated by the GATA1-PU.1 network: a classical example of a bistable cell-fate system. We show that, for a wide range of cell communication topologies, even subtle changes in signaling can have pronounced effects on cell-fate decisions. We go on to show how cell-cell coupling through signaling can spontaneously break the symmetry of a homogenous cell population. Noise, both intrinsic and extrinsic, shapes the decision landscape profoundly, and affects the transcriptional dynamics underlying this important hematopoietic cell-fate decision-making system. This article has an associated 'The people behind the papers' interview.

Age-dependent effects of Igf2bp2 on gene regulation, function, and aging of hematopoietic stem cells in mice.

Increasing evidence links metabolism, protein synthesis, and growth signaling to impairments in the function of hematopoietic stem and progenitor cells (HSPCs) during aging. The Lin28b/Hmga2 pathway controls tissue development, and the postnatal downregulation of this pathway limits the self-renewal of adult vs fetal hematopoietic stem cells (HSCs). Igf2bp2 is an RNA binding protein downstream of Lin28b/Hmga2, which regulates messenger RNA stability and translation. The role of Igf2bp2 in HSC aging is unknown. In this study, an analysis of wild-type and Igf2bp2 knockout mice showed that Igf2bp2 regulates oxidative metabolism in HSPCs and the expression of metabolism, protein synthesis, and stemness-related genes in HSCs of young mice. Interestingly, Igf2bp2 expression and function strongly declined in aging HSCs. In young mice, Igf2bp2 deletion mimicked aging-related changes in HSCs, including changes in Igf2bp2 target gene expression and impairment of colony formation and repopulation capacity. In aged mice, Igf2bp2 gene status had no effect on these parameters in HSCs. Unexpectedly, Igf2bp2-deficient mice exhibited an amelioration of the aging-associated increase in HSCs and myeloid-skewed differentiation. The results suggest that Igf2bp2 controls mitochondrial metabolism, protein synthesis, growth, and stemness of young HSCs, which is necessary for full HSC function during young adult age. However, Igf2bp2 gene function is lost during aging, and it appears to contribute to HSC aging in 2 ways: the aging-related loss of Igf2bp2 gene function impairs the growth and repopulation capacity of aging HSCs, and the activity of Igf2bp2 at a young age contributes to aging-associated HSC expansion and myeloid skewing.

RVAgene: generative modeling of gene expression time series data.

MOTIVATION: Methods to model dynamic changes in gene expression at a genome-wide level are not currently sufficient for large (temporally rich or single-cell) datasets. Variational autoencoders offer means to characterize large datasets and have been used effectively to characterize features of single-cell datasets. Here, we extend these methods for use with gene expression time series data. RESULTS: We present RVAgene: a recurrent variational autoencoder to model gene expression dynamics. RVAgene learns to accurately and efficiently reconstruct temporal gene profiles. It also learns a low dimensional representation of the data via a recurrent encoder network that can be used for biological feature discovery, and from which we can generate new gene expression data by sampling the latent space. We test RVAgene on simulated and real biological datasets, including embryonic stem cell differentiation and kidney injury response dynamics. In all cases, RVAgene accurately reconstructed complex gene expression temporal profiles. Via cross validation, we show that a low-error latent space representation can be learnt using only a fraction of the data. Through clustering and gene ontology term enrichment analysis on the latent space, we demonstrate the potential of RVAgene for unsupervised discovery. In particular, RVAgene identifies new programs of shared gene regulation of Lox family genes in response to kidney injury. AVAILABILITY AND IMPLEMENTATION: All datasets analyzed in this manuscript are publicly available and have been published previously. RVAgene is available in Python, at GitHub: https://github.com/maclean-lab/RVAgene; Zenodo archive: http://doi.org/10.5281/zenodo.4271097. SUPPLEMENTARY INFORMATION: Supplementary data are available at Bioinformatics online.

An Ovol2-Zeb1 transcriptional circuit regulates epithelial directional migration and proliferation.

Directional migration is inherently important for epithelial tissue regeneration and repair, but how it is precisely controlled and coordinated with cell proliferation is unclear. Here, we report that Ovol2, a transcriptional repressor that inhibits epithelial-to-mesenchymal transition (EMT), plays a crucial role in adult skin epithelial regeneration and repair. Ovol2-deficient mice show compromised wound healing characterized by aberrant epidermal cell migration and proliferation, as well as delayed anagen progression characterized by defects in hair follicle matrix cell proliferation and subsequent differentiation. Epidermal keratinocytes and bulge hair follicle stem cells (Bu-HFSCs) lacking Ovol2 fail to expand in culture and display molecular alterations consistent with enhanced EMT and reduced proliferation. Live imaging of wound explants and Bu-HFSCs reveals increased migration speed but reduced directionality, and post-mitotic cell cycle arrest. Remarkably, simultaneous deletion of Zeb1 encoding an EMT-promoting factor restores directional migration to Ovol2-deficient Bu-HFSCs. Taken together, our findings highlight the important function of an Ovol2-Zeb1 EMT-regulatory circuit in controlling the directional migration of epithelial stem and progenitor cells to facilitate adult skin epithelial regeneration and repair.

Cell lineage and communication network inference via optimization for single-cell transcriptomics.

The use of single-cell transcriptomics has become a major approach to delineate cell subpopulations and the transitions between them. While various computational tools using different mathematical methods have been developed to infer clusters, marker genes, and cell lineage, none yet integrate these within a mathematical framework to perform multiple tasks coherently. Such coherence is critical for the inference of cell-cell communication, a major remaining challenge. Here, we present similarity matrix-based optimization for single-cell data analysis (SoptSC), in which unsupervised clustering, pseudotemporal ordering, lineage inference, and marker gene identification are inferred via a structured cell-to-cell similarity matrix. SoptSC then predicts cell-cell communication networks, enabling reconstruction of complex cell lineages that include feedback or feedforward interactions. Application of SoptSC to early embryonic development, epidermal regeneration, and hematopoiesis demonstrates robust identification of subpopulations, lineage relationships, and pseudotime, and prediction of pathway-specific cell communication patterns regulating processes of development and differentiation.

Transcriptomics reveals complex kinetics of dorsal-ventral patterning gene expression in the mandibular arch.

The mandibular or first pharyngeal arch forms the upper and lower jaws in all gnathostomes. A gene regulatory network that defines ventral, intermediate, and dorsal domains along the dorsal-ventral (D-V) axis of the arch has emerged from studies in zebrafish and mice, but the temporal dynamics of this process remain unclear. To define cell fate trajectories in the arches we have performed quantitative gene expression analyses of D-V patterning genes in pharyngeal arch primordia in zebrafish and mice. Using NanoString technology to measure transcript numbers per cell directly we show that, in many cases, genes expressed in similar D-V domains and induced by similar signals vary dramatically in their temporal profiles. This suggests that cellular responses to D-V patterning signals are likely shaped by the baseline kinetics of target gene expression. Furthermore, similarities in the temporal dynamics of genes that occupy distinct pathways suggest novel shared modes of regulation. Incorporating these gene expression kinetics into our computational models for the mandibular arch improves the accuracy of patterning, and facilitates temporal comparisons between species. These data suggest that the magnitude and timing of target gene expression help diversify responses to patterning signals during craniofacial development.

scEpath: energy landscape-based inference of transition probabilities and cellular trajectories from single-cell transcriptomic data.

Motivation: Single-cell RNA-sequencing (scRNA-seq) offers unprecedented resolution for studying cellular decision-making processes. Robust inference of cell state transition paths and probabilities is an important yet challenging step in the analysis of these data. Results: Here we present scEpath, an algorithm that calculates energy landscapes and probabilistic directed graphs in order to reconstruct developmental trajectories. We quantify the energy landscape using 'single-cell energy' and distance-based measures, and find that the combination of these enables robust inference of the transition probabilities and lineage relationships between cell states. We also identify marker genes and gene expression patterns associated with cell state transitions. Our approach produces pseudotemporal orderings that are-in combination-more robust and accurate than current methods, and offers higher resolution dynamics of the cell state transitions, leading to new insight into key transition events during differentiation and development. Moreover, scEpath is robust to variation in the size of the input gene set, and is broadly unsupervised, requiring few parameters to be set by the user. Applications of scEpath led to the identification of a cell-cell communication network implicated in early human embryo development, and novel transcription factors important for myoblast differentiation. scEpath allows us to identify common and specific temporal dynamics and transcriptional factor programs along branched lineages, as well as the transition probabilities that control cell fates. Availability and implementation: A MATLAB package of scEpath is available at https://github.com/sqjin/scEpath. Supplementary information: Supplementary data are available at Bioinformatics online.

Concise Review: Stem Cell Population Biology: Insights from Hematopoiesis.

Stem cells are fundamental to human life and offer great therapeutic potential, yet their biology remains incompletely-or in cases even poorly-understood. The field of stem cell biology has grown substantially in recent years due to a combination of experimental and theoretical contributions: the experimental branch of this work provides data in an ever-increasing number of dimensions, while the theoretical branch seeks to determine suitable models of the fundamental stem cell processes that these data describe. The application of population dynamics to biology is amongst the oldest applications of mathematics to biology, and the population dynamics perspective continues to offer much today. Here we describe the impact that such a perspective has made in the field of stem cell biology. Using hematopoietic stem cells as our model system, we discuss the approaches that have been used to study their key properties, such as capacity for self-renewal, differentiation, and cell fate lineage choice. We will also discuss the relevance of population dynamics in models of stem cells and cancer, where competition naturally emerges as an influential factor on the temporal evolution of cell populations. Stem Cells 2017;35:80-88.

Single Cell Phenotyping Reveals Heterogeneity Among Hematopoietic Stem Cells Following Infection.

The hematopoietic stem cell (HSC) niche provides essential microenvironmental cues for the production and maintenance of HSCs within the bone marrow. During inflammation, hematopoietic dynamics are perturbed, but it is not known whether changes to the HSC-niche interaction occur as a result. We visualize HSCs directly in vivo, enabling detailed analysis of the 3D niche dynamics and migration patterns in murine bone marrow following Trichinella spiralis infection. Spatial statistical analysis of these HSC trajectories reveals two distinct modes of HSC behavior: (a) a pattern of revisiting previously explored space and (b) a pattern of exploring new space. Whereas HSCs from control donors predominantly follow pattern (a), those from infected mice adopt both strategies. Using detailed computational analyses of cell migration tracks and life-history theory, we show that the increased motility of HSCs following infection can, perhaps counterintuitively, enable mice to cope better in deteriorating HSC-niche microenvironments following infection. Stem Cells 2017;35:2292-2304.

Bayesian inference of agent-based models: a tool for studying kidney branching morphogenesis.

The adult mammalian kidney has a complex, highly-branched collecting duct epithelium that arises as a ureteric bud sidebranch from an epithelial tube known as the nephric duct. Subsequent branching of the ureteric bud to form the collecting duct tree is regulated by subcellular interactions between the epithelium and a population of mesenchymal cells that surround the tips of outgrowing branches. The mesenchymal cells produce glial cell-line derived neurotrophic factor (GDNF), that binds with RET receptors on the surface of the epithelial cells to stimulate several subcellular pathways in the epithelium. Such interactions are known to be a prerequisite for normal branching development, although competing theories exist for their role in morphogenesis. Here we introduce the first agent-based model of ex vivo kidney uretic branching. Through comparison with experimental data, we show that growth factor-regulated growth mechanisms can explain early epithelial cell branching, but only if epithelial cell division depends in a switch-like way on the local growth factor concentration; cell division occurring only if the driving growth factor level exceeds a threshold. We also show how a recently-developed method, "Approximate Approximate Bayesian Computation", can be used to infer key model parameters, and reveal the dependency between the parameters controlling a growth factor-dependent growth switch. These results are consistent with a requirement for signals controlling proliferation and chemotaxis, both of which are previously identified roles for GDNF.

Parameter-free methods distinguish Wnt pathway models and guide design of experiments.

The canonical Wnt signaling pathway, mediated by ß-catenin, is crucially involved in development, adult stem cell tissue maintenance, and a host of diseases including cancer. We analyze existing mathematical models of Wnt and compare them to a new Wnt signaling model that targets spatial localization; our aim is to distinguish between the models and distill biological insight from them. Using Bayesian methods we infer parameters for each model from mammalian Wnt signaling data and find that all models can fit this time course. We appeal to algebraic methods (concepts from chemical reaction network theory and matroid theory) to analyze the models without recourse to specific parameter values. These approaches provide insight into aspects of Wnt regulation: the new model, via control of shuttling and degradation parameters, permits multiple stable steady states corresponding to stem-like vs. committed cell states in the differentiation hierarchy. Our analysis also identifies groups of variables that should be measured to fully characterize and discriminate between competing models, and thus serves as a guide for performing minimal experiments for model comparison.

The ecology in the hematopoietic stem cell niche determines the clinical outcome in chronic myeloid leukemia.

Chronic myeloid leukemia (CML) is a blood disease that disrupts normal function of the hematopoietic system. Despite the great progress made in terms of molecular therapies for CML, there remain large gaps in our understanding. By comparing mathematical models that describe CML progression and etiology we sought to identify those models that provide the best description of disease dynamics and their underlying mechanisms. Data for two clinical outcomes--disease remission or relapse--are considered, and we investigate these using Bayesian inference techniques throughout. We find that it is not possible to choose between the models based on fits to the data alone; however, by studying model predictions we can discard models that fail to take niche effects into account. More detailed analysis of the remaining models reveals mechanistic differences: for one model, leukemia stem cell dynamics determine the disease outcome; and for the other model disease progression is determined at the stage of progenitor cells, in particular by differences in progenitor death rates. This analysis also reveals distinct transient dynamics that will be experimentally accessible, but are currently at the limits of what is possible to measure. To resolve these differences we need to be able to probe the hematopoietic stem cell niche directly. Our analysis highlights the importance of further mapping of the bone marrow hematopoietic niche microenvironment as the "ecological" interactions between cells in this niche appear to be intricately linked to disease outcome.

Feedback mechanisms control coexistence in a stem cell model of acute myeloid leukaemia.

Haematopoietic stem cell dynamics regulate healthy blood cell production and are disrupted during leukaemia. Competition models of cellular species help to elucidate stem cell dynamics in the bone marrow microenvironment (or niche), and to determine how these dynamics impact leukaemia progression. Here we develop two models that target acute myeloid leukaemia with particular focus on the mechanisms that control proliferation via feedback signalling. It is within regions of parameter space permissive of coexistence that the effects of competition are most subtle and the clinical outcome least certain. Steady state and linear stability analyses identify parameter regions that allow for coexistence to occur, and allow us to characterise behaviour near critical points. Where analytical expressions are no longer informative, we proceed statistically and sample parameter space over a coexistence region. We find that the rates of proliferation and differentiation of healthy progenitors exert key control over coexistence. We also show that inclusion of a regulatory feedback onto progenitor cells promotes healthy haematopoiesis at the expense of leukaemia, and that - somewhat paradoxically - within the coexistence region feedback increases the sensitivity of the system to dominance by one lineage over another.

Developmental hematopoietic stem cell variation explains clonal hematopoiesis later in life.

Clonal hematopoiesis becomes increasingly common with age, but its cause is enigmatic because driver mutations are often absent. Serial observations infer weak selection indicating variants are acquired much earlier in life with unexplained initial growth spurts. Here we use fluctuating CpG methylation as a lineage marker to track stem cell clonal dynamics of hematopoiesis. We show, via the shared prenatal circulation of monozygotic twins, that weak selection conferred by stem cell variation created before birth can reliably yield clonal hematopoiesis later in life. Theory indicates weak selection will lead to dominance given enough time and large enough population sizes. Human hematopoiesis satisfies both these conditions. Stochastic loss of weakly selected variants is naturally prevented by the expansion of stem cell lineages during development. The dominance of stem cell clones created before birth is supported by blood fluctuating CpG methylation patterns that exhibit low correlation between unrelated individuals but are highly correlated between many elderly monozygotic twins. Therefore, clonal hematopoiesis driven by weak selection in later life appears to reflect variation created before birth.

Gastrulation-stage gene expression in Nipbl+/- mouse embryos foreshadows the development of syndromic birth defects.

In animal models, Nipbl -deficiency phenocopies gene expression changes and birth defects seen in Cornelia de Lange Syndrome (CdLS), the most common cause of which is Nipbl -haploinsufficiency. Previous studies in Nipbl +/- mice suggested that heart development is abnormal as soon as cardiogenic tissue is formed. To investigate this, we performed single-cell RNA-sequencing on wildtype (WT) and Nipbl +/- mouse embryos at gastrulation and early cardiac crescent stages. Nipbl +/- embryos had fewer mesoderm cells than WT and altered proportions of mesodermal cell subpopulations. These findings were associated with underexpression of genes implicated in driving specific mesodermal lineages. In addition, Nanog was found to be overexpressed in all germ layers, and many gene expression changes observed in Nipbl +/- embryos could be attributed to Nanog overexpression. These findings establish a link between Nipbl -deficiency, Nanog overexpression, and gene expression dysregulation/lineage misallocation, which ultimately manifest as birth defects in Nipbl +/- animals and CdLS. Teaser: Gene expression changes during gastrulation of Nipbl -deficient mice shed light on early origins of structural birth defects.

Gastrulation-stage gene expression in Nipbl+/- mouse embryos foreshadows the development of syndromic birth defects.

In animal models, Nipbl deficiency phenocopies gene expression changes and birth defects seen in Cornelia de Lange syndrome, the most common cause of which is Nipbl haploinsufficiency. Previous studies in Nipbl+/- mice suggested that heart development is abnormal as soon as cardiogenic tissue is formed. To investigate this, we performed single-cell RNA sequencing on wild-type and Nipbl+/- mouse embryos at gastrulation and early cardiac crescent stages. Nipbl+/- embryos had fewer mesoderm cells than wild-type and altered proportions of mesodermal cell subpopulations. These findings were associated with underexpression of genes implicated in driving specific mesodermal lineages. In addition, Nanog was found to be overexpressed in all germ layers, and many gene expression changes observed in Nipbl+/- embryos could be attributed to Nanog overexpression. These findings establish a link between Nipbl deficiency, Nanog overexpression, and gene expression dysregulation/lineage misallocation, which ultimately manifest as birth defects in Nipbl+/- animals and Cornelia de Lange syndrome.