A large-scale Boolean model of the rheumatoid arthritis fibroblast-like synoviocytes predicts drug synergies in the arthritic joint.

Rheumatoid arthritis (RA) is a complex autoimmune disease with an unknown aetiology. However, rheumatoid arthritis fibroblast-like synoviocytes (RA-FLS) play a significant role in initiating and perpetuating destructive joint inflammation by expressing immuno-modulating cytokines, adhesion molecules, and matrix remodelling enzymes. In addition, RA-FLS are primary drivers of inflammation, displaying high proliferative rates and an apoptosis-resistant phenotype. Thus, RA-FLS-directed therapies could become a complementary approach to immune-directed therapies by predicting the optimal conditions that would favour RA-FLS apoptosis, limit inflammation, slow the proliferation rate and minimise bone erosion and cartilage destruction. In this paper, we present a large-scale Boolean model for RA-FLS that consists of five submodels focusing on apoptosis, cell proliferation, matrix degradation, bone erosion and inflammation. The five-phenotype-specific submodels can be simulated independently or as a global model. In silico simulations and perturbations reproduced the expected biological behaviour of the system under defined initial conditions and input values. The model was then used to mimic the effect of mono or combined therapeutic treatments and predict novel targets and drug candidates through drug repurposing analysis.

Drug-target identification in COVID-19 disease mechanisms using computational systems biology approaches.

Introduction: The COVID-19 Disease Map project is a large-scale community effort uniting 277 scientists from 130 Institutions around the globe. We use high-quality, mechanistic content describing SARS-CoV-2-host interactions and develop interoperable bioinformatic pipelines for novel target identification and drug repurposing. Methods: Extensive community work allowed an impressive step forward in building interfaces between Systems Biology tools and platforms. Our framework can link biomolecules from omics data analysis and computational modelling to dysregulated pathways in a cell-, tissue- or patient-specific manner. Drug repurposing using text mining and AI-assisted analysis identified potential drugs, chemicals and microRNAs that could target the identified key factors. Results: Results revealed drugs already tested for anti-COVID-19 efficacy, providing a mechanistic context for their mode of action, and drugs already in clinical trials for treating other diseases, never tested against COVID-19. Discussion: The key advance is that the proposed framework is versatile and expandable, offering a significant upgrade in the arsenal for virus-host interactions and other complex pathologies.

BioSimulators: a central registry of simulation engines and services for recommending specific tools.

Computational models have great potential to accelerate bioscience, bioengineering, and medicine. However, it remains challenging to reproduce and reuse simulations, in part, because the numerous formats and methods for simulating various subsystems and scales remain siloed by different software tools. For example, each tool must be executed through a distinct interface. To help investigators find and use simulation tools, we developed BioSimulators (https://biosimulators.org), a central registry of the capabilities of simulation tools and consistent Python, command-line and containerized interfaces to each version of each tool. The foundation of BioSimulators is standards, such as CellML, SBML, SED-ML and the COMBINE archive format, and validation tools for simulation projects and simulation tools that ensure these standards are used consistently. To help modelers find tools for particular projects, we have also used the registry to develop recommendation services. We anticipate that BioSimulators will help modelers exchange, reproduce, and combine simulations.

UPMaBoSS: A Novel Framework for Dynamic Cell Population Modeling.

Mathematical modeling aims at understanding the effects of biological perturbations, suggesting ways to intervene and to reestablish proper cell functioning in diseases such as cancer or in autoimmune disorders. This is a difficult task for obvious reasons: the level of details needed to describe the intra-cellular processes involved, the numerous interactions between cells and cell types, and the complex dynamical properties of such populations where cells die, divide and interact constantly, to cite a few. Another important difficulty comes from the spatial distribution of these cells, their diffusion and motility. All of these aspects cannot be easily resolved in a unique mathematical model or with a unique formalism. To cope with some of these issues, we introduce here a novel framework, UPMaBoSS (for Update Population MaBoSS), dedicated to modeling dynamic populations of interacting cells. We rely on the preexisting tool MaBoSS, which enables probabilistic simulations of cellular networks. A novel software layer is added to account for cell interactions and population dynamics, but without considering the spatial dimension. This modeling approach can be seen as an intermediate step towards more complex spatial descriptions. We illustrate our methodology by means of a case study dealing with TNF-induced cell death. Interestingly, the simulation of cell population dynamics with UPMaBoSS reveals a mechanism of resistance triggered by TNF treatment. Relatively easy to encode, UPMaBoSS simulations require only moderate computational power and execution time. To ease the reproduction of simulations, we provide several Jupyter notebooks that can be accessed within the CoLoMoTo Docker image, which contains all software and models used for this study.

Sex Dimorphism of Nonalcoholic Fatty Liver Disease (NAFLD) in Pparg-Null Mice.

Men with nonalcoholic fatty liver disease (NAFLD) are more exposed to nonalcoholic steatohepatitis (NASH) and liver fibrosis than women. However, the underlying molecular mechanisms of NALFD sex dimorphism are unclear. We combined gene expression, histological and lipidomic analyses to systematically compare male and female liver steatosis. We characterized hepatosteatosis in three independent mouse models of NAFLD, ob/ob and lipodystrophic fat-specific (PpargFΔ/Δ) and whole-body PPARγ-null (PpargΔ/Δ) mice. We identified a clear sex dimorphism occurring only in PpargΔ/Δ mice, with females showing macro- and microvesicular hepatosteatosis throughout their entire life, while males had fewer lipid droplets starting from 20 weeks. This sex dimorphism in hepatosteatosis was lost in gonadectomized PpargΔ/Δ mice. Lipidomics revealed hepatic accumulation of short and highly saturated TGs in females, while TGs were enriched in long and unsaturated hydrocarbon chains in males. Strikingly, sex-biased genes were particularly perturbed in both sexes, affecting lipid metabolism, drug metabolism, inflammatory and cellular stress response pathways. Most importantly, we found that the expression of key sex-biased genes was severely affected in all the NAFLD models we tested. Thus, hepatosteatosis strongly affects hepatic sex-biased gene expression. With NAFLD increasing in prevalence, this emphasizes the urgent need to specifically address the consequences of this deregulation in humans.

Logical modelling of in vitro differentiation of human monocytes into dendritic cells unravels novel transcriptional regulatory interactions.

Dendritic cells (DCs) are the major specialized antigen-presenting cells, thereby connecting innate and adaptive immunity. Because of their role in establishing adaptive immunity, they constitute promising targets for immunotherapy. Monocytes can differentiate into DCs in vitro in the presence of colony-stimulating factor 2 (CSF2) and interleukin 4 (IL4), activating four signalling pathways (MAPK, JAK/STAT, NFKB and PI3K). However, the downstream transcriptional programme responsible for DC differentiation from monocytes (moDCs) remains unknown. By analysing the scientific literature on moDC differentiation, we established a preliminary logical model that helped us identify missing information regarding the activation of genes responsible for this differentiation, including missing targets for key transcription factors (TFs). Using ChIP-seq and RNA-seq data from the Blueprint consortium, we defined active and inactive promoters, together with differentially expressed genes in monocytes, moDCs and macrophages, which correspond to an alternative cell fate. We then used this functional genomic information to predict novel targets for previously identified TFs. By integrating this information, we refined our model and recapitulated the main established facts regarding moDC differentiation. Prospectively, the resulting model should be useful to develop novel immunotherapies targeting moDCs.

Comparison of SYK Signaling Networks Reveals the Potential Molecular Determinants of Its Tumor-Promoting and Suppressing Functions.

Spleen tyrosine kinase (SYK) can behave as an oncogene or a tumor suppressor, depending on the cell and tissue type. As pharmacological SYK inhibitors are currently evaluated in clinical trials, it is important to gain more information on the molecular mechanisms underpinning these opposite roles. To this aim, we reconstructed and compared its signaling networks using phosphoproteomic data from breast cancer and Burkitt lymphoma cell lines where SYK behaves as a tumor suppressor and promoter. Bioinformatic analyses allowed for unveiling the main differences in signaling pathways, network topology and signal propagation from SYK to its potential effectors. In breast cancer cells, the SYK target-enriched signaling pathways included intercellular adhesion and Hippo signaling components that are often linked to tumor suppression. In Burkitt lymphoma cells, the SYK target-enriched signaling pathways included molecules that could play a role in SYK pro-oncogenic function in B-cell lymphomas. Several protein interactions were profoundly rewired in the breast cancer network compared with the Burkitt lymphoma network. These data demonstrate that proteomic profiling combined with mathematical network modeling allows untangling complex pathway interplays and revealing difficult to discern interactions among the SYK pathways that positively and negatively affect tumor formation and progression.

Benchmarking joint multi-omics dimensionality reduction approaches for the study of cancer.

High-dimensional multi-omics data are now standard in biology. They can greatly enhance our understanding of biological systems when effectively integrated. To achieve proper integration, joint Dimensionality Reduction (jDR) methods are among the most efficient approaches. However, several jDR methods are available, urging the need for a comprehensive benchmark with practical guidelines. We perform a systematic evaluation of nine representative jDR methods using three complementary benchmarks. First, we evaluate their performances in retrieving ground-truth sample clustering from simulated multi-omics datasets. Second, we use TCGA cancer data to assess their strengths in predicting survival, clinical annotations and known pathways/biological processes. Finally, we assess their classification of multi-omics single-cell data. From these in-depth comparisons, we observe that intNMF performs best in clustering, while MCIA offers an effective behavior across many contexts. The code developed for this benchmark study is implemented in a Jupyter notebook-multi-omics mix (momix)-to foster reproducibility, and support users and future developers.

Interplay between SMAD2 and STAT5A is a critical determinant of IL-17A/IL-17F differential expression.

Interleukins (IL)-17A and F are critical cytokines in anti-microbial immunity but also contribute to auto-immune pathologies. Recent evidence suggests that they may be differentially produced by T-helper (Th) cells, but the underlying mechanisms remain unknown. To address this question, we built a regulatory graph integrating all reported upstream regulators of IL-17A and F, completed by ChIP-seq data analyses. The resulting regulatory graph encompasses 82 components and 136 regulatory links. The graph was then supplemented by logical rules calibrated with original flow cytometry data using naive CD4+ T cells, in conditions inducing IL-17A or IL-17F. The model displays specific stable states corresponding to virtual phenotypes explaining IL-17A and IL-17F differential regulation across eight cytokine stimulatory conditions. Our model analysis points to the transcription factors NFAT2A, STAT5A and SMAD2 as key regulators of the differential expression of IL-17A and IL-17F, with STAT5A controlling IL-17F expression, and an interplay of NFAT2A, STAT5A and SMAD2 controlling IL-17A expression. We experimentally observed that the production of IL-17A was correlated with an increase of SMAD2 transcription, and the expression of IL-17F correlated with an increase of BLIMP-1 transcription, together with an increase of STAT5A expression (mRNA), as predicted by our model. Interestingly, RORγt presumably plays a more determinant role in IL-17A expression as compared to IL-17F expression. In conclusion, we propose the first mechanistic model accounting for the differential expression of IL-17A and F in Th cells, providing a basis to design novel therapeutic interventions in auto-immune and inflammatory diseases.

Setting the basis of best practices and standards for curation and annotation of logical models in biology-highlights of the [BC]2 2019 CoLoMoTo/SysMod Workshop.

The fast accumulation of biological data calls for their integration, analysis and exploitation through more systematic approaches. The generation of novel, relevant hypotheses from this enormous quantity of data remains challenging. Logical models have long been used to answer a variety of questions regarding the dynamical behaviours of regulatory networks. As the number of published logical models increases, there is a pressing need for systematic model annotation, referencing and curation in community-supported and standardised formats. This article summarises the key topics and future directions of a meeting entitled 'Annotation and curation of computational models in biology', organised as part of the 2019 [BC]2 conference. The purpose of the meeting was to develop and drive forward a plan towards the standardised annotation of logical models, review and connect various ongoing projects of experts from different communities involved in the modelling and annotation of molecular biological entities, interactions, pathways and models. This article defines a roadmap towards the annotation and curation of logical models, including milestones for best practices and minimum standard requirements.

Dynamical Boolean Modeling of Immunogenic Cell Death.

As opposed to the standard tolerogenic apoptosis, immunogenic cell death (ICD) constitutes a type of cellular demise that elicits an adaptive immune response. ICD has been characterized in malignant cells following cytotoxic interventions, such as chemotherapy or radiotherapy. Briefly, ICD of cancer cells releases some stress/danger signals that attract and activate dendritic cells (DCs). The latter can then engulf and cross-present tumor antigens to T lymphocytes, thus priming a cancer-specific immunity. This series of reactions works as a positive feedback loop where the antitumor immunity further improves the therapeutic efficacy by targeting cancer cells spared by the cytotoxic agent. However, not all chemotherapeutic drugs currently approved for cancer treatment are able to stimulate bona fide ICD: some commonly used agents, such as cisplatin or 5-fluorouracil, are unable to activate all features of ICD. Therefore, a better characterization of the process could help identify some gene or protein candidates to target pharmacologically and suggest combinations of drugs that would favor/increase antitumor immune response. To this end, we have built a mathematical model of the major cell types that intervene in ICD, namely cancer cells, DCs, CD8+ and CD4+ T cells. Our model not only integrates intracellular mechanisms within each individual cell entity, but also incorporates intercellular communications between them. The resulting cell population model recapitulates key features of the dynamics of ICD after an initial treatment, in particular the time-dependent size of the different cell types. The model is based on a discrete Boolean formalism and is simulated by means of a software tool, UPMaBoSS, which performs stochastic simulations with continuous time, considering the dynamics of the system at the cell population level with appropriate timing of events, and accounting for death and division of each cell type. With this model, the time scales of some of the processes involved in ICD, which are challenging to measure experimentally, have been predicted. In addition, our model analysis led to the identification of actionable targets for boosting ICD-induced antitumor response. All computational analyses and results are compiled in interactive notebooks which cover the presentation of the network structure, model simulations, and parameter sensitivity analyses.

Computational Verification of Large Logical Models-Application to the Prediction of T Cell Response to Checkpoint Inhibitors.

At the crossroad between biology and mathematical modeling, computational systems biology can contribute to a mechanistic understanding of high-level biological phenomenon. But as knowledge accumulates, the size and complexity of mathematical models increase, calling for the development of efficient dynamical analysis methods. Here, we propose the use of two approaches for the development and analysis of complex cellular network models. A first approach, called "model verification" and inspired by unitary testing in software development, enables the formalization and automated verification of validation criteria for whole models or selected sub-parts. When combined with efficient analysis methods, this approach is suitable for continuous testing, thereby greatly facilitating model development. A second approach, called "value propagation," enables efficient analytical computation of the impact of specific environmental or genetic conditions on the dynamical behavior of some models. We apply these two approaches to the delineation and the analysis of a comprehensive model for T cell activation, taking into account CTLA4 and PD-1 checkpoint inhibitory pathways. While model verification greatly eases the delineation of logical rules complying with a set of dynamical specifications, propagation provides interesting insights into the different potential of CTLA4 and PD-1 immunotherapies. Both methods are implemented and made available in the all-inclusive CoLoMoTo Docker image, while the different steps of the model analysis are fully reported in two companion interactive jupyter notebooks, thereby ensuring the reproduction of our results.

SBML Level 3: an extensible format for the exchange and reuse of biological models.

Systems biology has experienced dramatic growth in the number, size, and complexity of computational models. To reproduce simulation results and reuse models, researchers must exchange unambiguous model descriptions. We review the latest edition of the Systems Biology Markup Language (SBML), a format designed for this purpose. A community of modelers and software authors developed SBML Level 3 over the past decade. Its modular form consists of a core suited to representing reaction-based models and packages that extend the core with features suited to other model types including constraint-based models, reaction-diffusion models, logical network models, and rule-based models. The format leverages two decades of SBML and a rich software ecosystem that transformed how systems biologists build and interact with models. More recently, the rise of multiscale models of whole cells and organs, and new data sources such as single-cell measurements and live imaging, has precipitated new ways of integrating data with models. We provide our perspectives on the challenges presented by these developments and how SBML Level 3 provides the foundation needed to support this evolution.

Automated inference of Boolean models from molecular interaction maps using CaSQ.

MOTIVATION: Molecular interaction maps have emerged as a meaningful way of representing biological mechanisms in a comprehensive and systematic manner. However, their static nature provides limited insights to the emerging behaviour of the described biological system under different conditions. Computational modelling provides the means to study dynamic properties through in silico simulations and perturbations. We aim to bridge the gap between static and dynamic representations of biological systems with CaSQ, a software tool that infers Boolean rules based on the topology and semantics of molecular interaction maps built with CellDesigner. RESULTS: We developed CaSQ by defining conversion rules and logical formulas for inferred Boolean models according to the topology and the annotations of the starting molecular interaction maps. We used CaSQ to produce executable files of existing molecular maps that differ in size, complexity and the use of Systems Biology Graphical Notation (SBGN) standards. We also compared, where possible, the manually built logical models corresponding to a molecular map to the ones inferred by CaSQ. The tool is able to process large and complex maps built with CellDesigner (either following SBGN standards or not) and produce Boolean models in a standard output format, Systems Biology Marked Up Language-qualitative (SBML-qual), that can be further analyzed using popular modelling tools. References, annotations and layout of the CellDesigner molecular map are retained in the obtained model, facilitating interoperability and model reusability. AVAILABILITY AND IMPLEMENTATION: The present tool is available online: https://lifeware.inria.fr/∼soliman/post/casq/ and distributed as a Python package under the GNU GPLv3 license. The code can be accessed here: https://gitlab.inria.fr/soliman/casq. SUPPLEMENTARY INFORMATION: Supplementary data are available at Bioinformatics online.

PTPN13 induces cell junction stabilization and inhibits mammary tumor invasiveness.

Clinical data suggest that the protein tyrosine phosphatase PTPN13 exerts an anti-oncogenic effect. Its exact role in tumorigenesis remains, however, unclear due to its negative impact on FAS receptor-induced apoptosis. Methods: We crossed transgenic mice deleted for PTPN13 phosphatase activity with mice that overexpress human HER2 to assess the exact role of PTPN13 in tumor development and aggressiveness. To determine the molecular mechanism underlying the PTPN13 tumor suppressor activity we developed isogenic clones of the aggressive human breast cancer cell line MDA-MB-231 overexpressing either wild type or a catalytically-inactive mutant PTPN13 and subjected these to phosphoproteomic and gene ontology analyses. We investigated the PTPN13 consequences on cell aggressiveness using wound healing and Boyden chamber assays, on intercellular adhesion using videomicroscopy, cell aggregation assay and immunofluorescence. Results: The development, growth and invasiveness of breast tumors were strongly increased by deletion of the PTPN13 phosphatase activity in transgenic mice. We observed that PTPN13 phosphatase activity is required to inhibit cell motility and invasion in the MDA-MB-231 cell line overexpressing PTPN13. In vivo, the negative PTPN13 effect on tumor invasiveness was associated with a mesenchymal-to-epithelial transition phenotype in athymic mice xenografted with PTPN13-overexpressing MDA-MB-231 cells, as well as in HER2-overexpressing mice with wild type PTPN13, compared to HER2-overexpressing mice that lack PTPN13 phosphatase activity. Phosphoproteomic and gene ontology analyses indicated a role of PTPN13 in the regulation of intercellular junction-related proteins. Finally, protein localization studies in MDA-MB-231 cells and HER2-overexpressing mice tumors confirmed that PTPN13 stabilizes intercellular adhesion and promotes desmosome formation. Conclusions: These data provide the first evidence for the negative role of PTPN13 in breast tumor invasiveness and highlight its involvement in cell junction stabilization.

Contribution of ROS and metabolic status to neonatal and adult CD8+ T cell activation.

In neonatal T cells, a low response to infection contributes to a high incidence of morbidity and mortality of neonates. Here we have evaluated the impact of the cytoplasmic and mitochondrial levels of Reactive Oxygen Species of adult and neonatal CD8+ T cells on their activation potential. We have also constructed a logical model connecting metabolism and ROS with T cell signaling. Our model indicates the interplay between antigen recognition, ROS and metabolic status in T cell responses. This model displays alternative stable states corresponding to different cell fates, i.e. quiescent, activated and anergic states, depending on ROS levels. Stochastic simulations with this model further indicate that differences in ROS status at the cell population level contribute to the lower activation rate of neonatal, compared to adult, CD8+ T cells upon TCR engagement. These results are relevant for neonatal health care. Our model can serve to analyze the impact of metabolic shift during cancer in which, similar to neonatal cells, a high glycolytic rate and low concentrations of glutamine and arginine promote tumor tolerance.

Network Reconstruction and Significant Pathway Extraction Using Phosphoproteomic Data from Cancer Cells.

Protein phosphorylation acts as an efficient switch controlling deregulated key signaling pathway in cancer. Computational biology aims to address the complexity of reconstructed networks but overrepresents well-known proteins and lacks information on less-studied proteins. A bioinformatic tool to reconstruct and select relatively small networks that connect signaling proteins to their targets in specific contexts is developed. It enables to propose and validate new signaling axes of the Syk kinase. To validate the potency of the tool, it is applied to two phosphoproteomic studies on oncogenic mutants of the well-known phosphatidyl-inositol 3-kinase (PIK3CA) and the unfamiliar Src-related tyrosine kinase lacking C-terminal regulatory tyrosine and N-terminal myristoylation sites (SRMS) kinase. By combining network reconstruction and signal propagation, comprehensive signaling networks from large-scale experimental data are built and multiple molecular paths from these kinases to their targets are extracted. Specific paths from two distinct PIK3CA mutants are retrieved, and their differential impact on the HER3 receptor kinase is explained. In addition, to address the missing connectivities of the SRMS kinase to its targets in interaction pathway databases, phospho-tyrosine and phospho-serine/threonine proteomic data are integrated. The resulting SRMS-signaling network comprises casein kinase 2, thereby validating its currently suggested role downstream of SRMS. The computational pipeline is publicly available, and contains a user-friendly graphical interface (http://doi.org/10.5281/zenodo.3333687).

System analysis of cross-talk between nuclear receptors reveals an opposite regulation of the cell cycle by LXR and FXR in human HepaRG liver cells.

Transcriptional regulations exert a critical control of metabolic homeostasis. In particular, the nuclear receptors (NRs) are involved in regulating numerous pathways of the intermediate metabolism. The purpose of the present study was to explore in liver cells the interconnectedness between three of them, LXR, FXR, and PPARα, all three known to act on lipid and glucose metabolism, and also on inflammation. The human cell line HepaRG was selected for its best proximity to human primary hepatocytes. Global gene expression of differentiated HepaRG cells was assessed after 4 hours and 24 hours of exposure to GW3965 (LXR agonist), GW7647 (PPARα agonist), and GW4064 and CDCA (FXR synthetic and natural agonist, respectively). Our work revealed that, contrary to our expectations, NR specificity is largely present at the level of target genes, with a smaller than expected overlap of the set of genes targeted by the different NRs. It also highlighted the much broader activity of the synthetic FXR ligand compared to CDCA. More importantly, our results revealed that activation of FXR has a pro-proliferative effect and decreases the number of tetraploid (or binucleated) hepatocytes, while LXR inhibits the cell cycle progression, inducing hepatocyte differentiation and an increase in tetraploidism. Conclusion: these results highlight the importance of analyzing the different NR activities in a context allowing a direct confrontation of each receptor outcome, and reveals the opposite role of FXR and LXR in hepatocyte cells division and maturation.

BioLQM: A Java Toolkit for the Manipulation and Conversion of Logical Qualitative Models of Biological Networks.

Here we introduce bioLQM, a new Java software toolkit for the conversion, modification, and analysis of Logical Qualitative Models of biological regulatory networks. BioLQM provides core modeling operations as building blocks for the development of integrated modeling software, or for the assembly of heterogeneous analysis workflows involving several complementary tools. Based on the definition of multi-valued logical models, bioLQM implements import and export facilities, notably for the recent SBML qual exchange format, as well as for formats used by several popular tools, facilitating the design of workflows combining these tools. Model modifications enable the definition of various perturbations, as well as model reduction, easing the analysis of large models. Another modification enables the study of multi-valued models with tools limited to the Boolean case. Finally, bioLQM provides a framework for the development of novel analysis tools. The current version implements various updating modes for model simulation (notably synchronous, asynchronous, and random asynchronous), as well as some static analysis features for the identification of attractors. The bioLQM software can be integrated into analysis workflows through command line and scripting interfaces. As a Java library, it further provides core data structures to the GINsim and EpiLog interactive tools, which supply graphical interfaces and additional analysis methods for cellular and multi-cellular qualitative models.

GOATOOLS: A Python library for Gene Ontology analyses.

The biological interpretation of gene lists with interesting shared properties, such as up- or down-regulation in a particular experiment, is typically accomplished using gene ontology enrichment analysis tools. Given a list of genes, a gene ontology (GO) enrichment analysis may return hundreds of statistically significant GO results in a "flat" list, which can be challenging to summarize. It can also be difficult to keep pace with rapidly expanding biological knowledge, which often results in daily changes to any of the over 47,000 gene ontologies that describe biological knowledge. GOATOOLS, a Python-based library, makes it more efficient to stay current with the latest ontologies and annotations, perform gene ontology enrichment analyses to determine over- and under-represented terms, and organize results for greater clarity and easier interpretation using a novel GOATOOLS GO grouping method. We performed functional analyses on both stochastic simulation data and real data from a published RNA-seq study to compare the enrichment results from GOATOOLS to two other popular tools: DAVID and GOstats. GOATOOLS is freely available through GitHub: https://github.com/tanghaibao/goatools .