Search

The major genetic risk factors for Hirschsprung disease (HSCR) are three common polymorphisms within cis-regulatory elements (CREs) of the receptor tyrosine kinase gene RET, which reduce its expression during enteric nervous system (ENS) development. These risk variants attenuate binding of the transcription factors RARB, GATA2, and SOX10 to their cognate CREs, reduce RET gene expression, and dysregulate other ENS and HSCR genes in the RET–EDNRB gene regulatory network (GRN). Here, we use siRNA, ChIP, and CRISPR-Cas9 deletion analyses in the SK-N-SH cell line to ask how many additional HSCR-associated risk variants reside in RET CREs that affect its gene expression. We identify 22 HSCR-associated variants in candidate RET CREs, of which seven have differential allele-specific in vitro enhancer activity, and four of these seven affect RET gene expression; of these, two enhancers are bound by the transcription factor PAX3. We also show that deleting multiple variant-containing enhancers leads to synergistic effects on RET gene expression. These, coupled with our prior results, show that common sequence variants in at least 10 RET enhancers affect HSCR risk, seven with experimental evidence of affecting RET gene expression, extending the known RET–EDNRB GRN to reveal an extensive regulatory code modulating disease risk at a single gene.

We used capped analysis of gene expression with sequencing (CAGE-seq) to profile eRNA expression and enhancer activity during embryogenesis of a model echinoderm: the sea urchin, Strongylocentrotus purpuratus. We identified more than 18,000 enhancers that were active in mature oocytes and developing embryos and documented a burst of enhancer activation during cleavage and early blastula stages. We found that a large fraction (73.8%) of all enhancers active during the first 48 h of embryogenesis were hyperaccessible no later than the 128-cell stage and possibly even earlier. Most enhancers were located near gene bodies, and temporal patterns of eRNA expression tended to parallel those of nearby genes. Furthermore, enhancers near lineage-specific genes contained signatures of inputs from developmental gene regulatory networks deployed in those lineages. A large fraction (60%) of sea urchin enhancers previously shown to be active in transgenic reporter assays was associated with eRNA expression. Moreover, a large fraction (50%) of a representative subset of enhancers identified by eRNA profiling drove tissue-specific gene expression in isolation when tested by reporter assays. Our findings provide an atlas of developmental enhancers in a model sea urchin and support the utility of eRNA profiling as a tool for enhancer discovery and regulatory biology. The data generated in this study are available at Echinobase, the public database of information related to echinoderm genomics.

Single-cell RNA sequencing (scRNA-seq) has emerged as a powerful experimental approach to study cellular heterogeneity. One of the challenges in scRNA-seq data analysis is integrating different types of biological data to consistently recognize discrete biological functions and regulatory mechanisms of cells, such as transcription factor activities and gene regulatory networks in distinct cell populations. We have developed an approach to infer transcription factor activities from scRNA-seq data that leverages existing biological data on transcription factor binding sites. The Bayesian inference transcription factor activity model (BITFAM) integrates ChIP-seq transcription factor binding information into scRNA-seq data analysis. We show that the inferred transcription factor activities for key cell types identify regulatory transcription factors that are known to mechanistically control cell function and cell fate. The BITFAM approach not only identifies biologically meaningful transcription factor activities, but also provides valuable insights into underlying transcription factor regulatory mechanisms.

Multicellular organisms, such as plants, are characterized by highly specialized and tightly regulated cell populations, establishing specific morphological structures and executing distinct functions. Gene regulatory networks (GRNs) describe condition-specific interactions of transcription factor (TF) regulating the expression of target genes, underpinning these specific functions. As efficient and validated methods to identify cell-type specific GRNs from single-cell data in plants are lacking, limiting our understanding of the organization of specific cell-types in both model species and crops, we developed MINI-EX (Motif-Informed Network Inference based on single-cell Expression data), an integrative approach to infer cell-type specific networks in plants. MINI-EX uses single-cell transcriptomic data to define expression-based networks and integrates TF motif information to filter the inferred regulons, resulting in networks with increased accuracy. Next, regulons are assigned to different cell-types, leveraging cell-specific expression, and candidate regulators are prioritized using network centrality measures, functional annotations, and expression specificity. This embedded prioritization strategy offers a unique and efficient means to unravel signaling cascades in specific cell-types controlling a biological process of interest. We demonstrate MINI-EX’s stability towards input data sets with low number of cells and its robustness towards missing data, and we show it infers state-of-the-art networks with a better performance compared to related single-cell network tools. MINI-EX successfully identifies key regulators controlling root development in Arabidopsis and rice, Arabidopsis leaf development, and governing ear development in maize, enhancing our understanding of cell-type specific regulation and unraveling the role of different regulators controlling the development of specific cell-types in plants.

Single-cell RNA sequencing (scRNA-seq) is a widely used method for identifying cell types and trajectories in biologically heterogeneous samples, but it is limited in its detection and quantification of lowly expressed genes. This results in missing important biological signals, such as the expression of key transcription factors (TFs) driving cellular differentiation. We show that targeted sequencing of ∼1000 TFs (scCapture-seq) in iPSC-derived neuronal cultures greatly improves the biological information garnered from scRNA-seq. Increased TF resolution enhanced cell type identification, developmental trajectories, and gene regulatory networks. This allowed us to resolve differences among neuronal populations, which were generated in two different laboratories using the same differentiation protocol. ScCapture-seq improved TF-gene regulatory network inference and thus identified divergent patterns of neurogenesis into either excitatory cortical neurons or inhibitory interneurons. Furthermore, scCapture-seq revealed a role for of retinoic acid signaling in the developmental divergence between these different neuronal populations. Our results show that TF targeting improves the characterization of human cellular models and allows identification of the essential differences between cellular populations, which would otherwise be missed in traditional scRNA-seq. scCapture-seq TF targeting represents a cost-effective enhancement of scRNA-seq, which could be broadly applied to improve scRNA-seq resolution.

RNA-binding proteins (RNA-BPs) play critical roles in development and disease to regulate gene expression. However, genome-wide identification of their targets in primary human cells has been challenging. Here, we applied a modified CLIP-seq strategy to identify genome-wide targets of the FMRP translational regulator 1 (FMR1), a brain-enriched RNA-BP, whose deficiency leads to Fragile X Syndrome (FXS), the most prevalent inherited intellectual disability. We identified FMR1 targets in human dorsal and ventral forebrain neural progenitors and excitatory and inhibitory neurons differentiated from human pluripotent stem cells. In parallel, we measured the transcriptomes of the same four cell types upon FMR1 gene deletion. We discovered that FMR1 preferentially binds long transcripts in human neural cells. FMR1 targets include genes unique to human neural cells and associated with clinical phenotypes of FXS and autism. Integrative network analysis using graph diffusion and multitask clustering of FMR1 CLIP-seq and transcriptional targets reveals critical pathways regulated by FMR1 in human neural development. Our results demonstrate that FMR1 regulates a common set of targets among different neural cell types but also operates in a cell type–specific manner targeting distinct sets of genes in human excitatory and inhibitory neural progenitors and neurons. By defining molecular subnetworks and validating specific high-priority genes, we identify novel components of the FMR1 regulation program. Our results provide new insights into gene regulation by a critical neuronal RNA-BP in human neurodevelopment.

With the wide availability of single cell RNA-seq (scRNA-seq) technology, population-scale scRNA-seq datasets across multiple individuals and time points are emerging. Although the immediate analysis of these datasets tend to focus on standard analysis of clustering and differential expression, leveraging the power of scRNA-seq data at the personalized dynamic gene co-expression network-level has the potential to unlock subject and/or time-specific network-level variation critical for phenotypic differences. Community detection from co-expression networks of multiple time points or conditions has been well-studied; however, none of the existing settings included networks from multiple subjects and multiple time points simultaneously. To address this, we develop MuDCoD for multi-subject community detection in personalized dynamic gene networks from single cell RNA sequencing. MuDCoD builds on the spectral clustering framework and promotes information sharing among the subjects' networks as well as networks at different time points. It clusters genes in the personalized dynamic gene networks and reveals gene communities that are variable or shared not only across time but also among subjects. Evaluation and benchmarking of MuDCoD against existing approaches reveal that MuDCoD effectively leverages apparent shared signals among networks of the subjects at individual time points, and performs robustly when there is no or little information sharing among the networks. Applications to population-scale scRNA-seq datasets of human-induced pluripotent stem cells during dopaminergic neuron differentiation and CD4+ T cell activation indicate that MuDCoD enables robust inference for identifying time-varying personalized gene modules. Our results illustrate how personalized dynamic community detection can aid in the exploration of subject-specific biological processes that vary across time. MuDCoD is publicly available at https://github.com/bo1929/MuDCoD as a Python package. Implementation includes simulation and real-data experiments together with extensive documentation.

Gestational diabetes mellitus (GDM) is a metabolic disorder during pregnancy. Numerous biomarkers have been identified that are linked with the occurrence and development of GDM. The aim of this investigation was to identify differentially expressed genes (DEGs) in GDM using a bioinformatics approach to elucidate their molecular pathogenesis. GDM associated expression profiling by high throughput sequencing dataset (GSE154377) was obtained from Gene Expression Omnibus (GEO) database including 28 normal pregnancy samples and 33 GDM samples. DEGs were identified using DESeq2. The gene ontology (GO) and REACTOME pathway enrichments of DEGs were performed by g:Profiler. Protein-protein interaction (PPI) networks were assembled with Cytoscape software and separated into modules using the PEWCC1 algorithm. MiRNA-hub gene regulatory network and TF-hub gene regulatory network were performed with the miRNet database and NetworkAnalyst database. Receiver Operating Characteristic (ROC) analyses was conducted to validate the hub genes. A total of 953 DEGs were identified, of which 478 DEGs were up regulated and 475 DEGs were down regulated. Furthermore, GO and REACTOME pathway enrichment analysis demonstrated that these DEGs were mainly enriched in multicellular organismal process, cell activation, formation of the cornified envelope and hemostasis. TRIM54, ELAVL2, PTN, KIT, BMPR1B, APP, SRC, ITGA4, RPA1 and ACTB were identified as key genes in the PPI network, miRNA-hub gene regulatory network and TF-hub gene regulatory network. TRIM54, ELAVL2, PTN, KIT, BMPR1B, APP, SRC, ITGA4, RPA1 and ACTB in GDM were validated using ROC analysis. This investigation provides further insights into the molecular pathogenesis of GDM, which might facilitate the diagnosis and treatment of GDM.

BackgroundAlthough the widespread expression of circular RNAs (circRNAs) has only been recognized recently, increasing evidence has suggested their important roles in health and disease. To identify clinically relevant circRNAs with potential for wound diagnosis and therapy, an in-depth characterization of circRNA expression in human healing and non-healing wounds is a prerequisite that has not been attained yet.MethodsWe collected wound-edge biopsies through the healing process of healthy donors and in chronic non-healing venous ulcers (VU). Paired total RNA- and small RNA-sequencing were performed to profile circRNAs, protein-coding mRNAs, and microRNA expression. We analyzed the co-expression relationship between circRNAs and mRNAs with weighted correlation network analysis (WGCNA) and constructed circRNA-microRNA-mRNA networks. For the circRNAs surfaced in the in-silico analysis, after validating their expression with RT-PCR and sequencing, we silenced hsa-CHST15_0003 and hsa-TNFRSF21_0001 expression in keratinocytes with siRNAs and studied their function with transcriptomic profiling and live-cell monitoring.ResultsOur study unravels the dynamically changed expression patterns of circRNAs during human skin wound healing and their abnormal expression signature in VU, which are presented as a searchable web resource (http://130.229.28.87/shiny/circRNA_wholebiopsy-shinyApp/). In silico analysis deciphers the circRNA-miRNAs-mRNA networks specific to the inflammatory and proliferative phases of wound repair and VU, the biological processes that circRNAs are involved, and the circRNAs that could act as miRNAs sponge in human wounds. Importantly, we found that hsa-CHST15_0003 and hsa-TNFRSF21_0001, two circRNAs upregulated in VU, hampered keratinocyte migration while promoting proliferation through modulating gene networks underpinning these cellular processes.ConclusionBy integrating circRNA, mRNA, and miRNA expression profiles in a unique collection of clinical samples, we identify the circRNAs that are relevant to human wound healing physiology and pathology. This study paves the way to decipher the functional significance of circRNAs in tissue repair.

While single-cell RNA sequencing (scRNA-seq) is widely used to profile gene expression, few methods are available to infer gene regulatory networks using scRNA-seq data. Here, we developed and extended IReNA (Integrated Regulatory Network Analysis) to perform regulatory network analysis using scRNA-seq profiles. Four features are developed for IReNA. First, regulatory networks are divided into different modules which represent distinct biological functions. Second, transcription factors significantly regulating each gene module can be identified. Third, regulatory relationships among modules can be inferred. Fourth, IReNA can integrate ATAC-seq data into regulatory network analysis. If ATAC-seq data is available, both transcription factor footprints and binding motifs are used to refine regulatory relationships among co-expressed genes. Using public datasets, we showed that integrated network analysis of scRNA-seq data with ATAC-seq data identified a higher fraction of known regulators than scRNA-seq data alone. Moreover, IReNA provided a better performance of network analysis than currently available methods. Beyond the reconstruction of regulatory networks, IReNA can modularize regulatory networks, and identify key regulators and significant regulatory relationships for modules, facilitating the systems-level understanding of biological regulatory mechanisms. The R package IReNA is available at https://github.com/jiang-junyao/IReNA. Key pointsO_LIWe have developed IReNA, a new method for reconstructing regulatory networks using either scRNA-seq and ATAC-seq data or scRNA-seq data alone. C_LIO_LIIReNA can establish modular regulatory networks to identify key regulatory factors and statistically significant regulatory relationships among modules. C_LIO_LIThrough analyzing three public scRNA-seq datasets, IReNA shows better performance on identifying known regulators than Rcistarget, the most widely used alternative method. C_LI

Heterosis describes the phenomenon whereby a hybrid population has higher fitness than an inbred population, and has previously been explained by either Mendelian dominance or overdominance, where it is generally assumed that one gene controls one trait. How-ever, recent studies have demonstrated that genes interact through a complex gene regulatory network (GRN). Furthermore, phenotypic variance due to noise is reportedly lower for heterozygotes, whereas the origin of such variance-related heterosis remains elusive. There-fore, a theoretical analysis linking heterosis to GRN evolution and stochastic gene expression dynamics is required. Here, we investigate heterosis related to fitness and phenotypic variance in a system with interacting genes, by numerically evolving diploid GRNs. According to the results, the heterozygote population exhibited higher fitness than the homozygote population, that is, fitness-related heterosis resulting from evolution. In addition, the heterozygote population expressed lower noise-related phenotypic variance in expression levels than the homozygous population, implying that the heterozygote population is more robust to noise. Furthermore, the distribution of the ratio of heterozygote phenotypic variance to homozygote phenotypic variance exhibited quantitative agreement with previous experimental results. By applying dominance and over-dominance to the gene expression pattern rather than only a single gene expression, we confirmed the correlation between heterosis and overdominance. We explain our results by proposing that the convex high-fitness region is evolutionarily shaped in the genetic space to gain noise robustness under genetic mixing through sexual reproduction.Significance StatementHeterosis, that is higher fitness in hybrid populations than inbred populations, is a long-standing problem in genetics, evolution, and breeding studies. Heterosis is not necessarily the result of a trait uniquely determined by a single gene, but likely involves stochasticity and interactions among genes. Through numerical evolution of the gene regulatory network, we demonstrate that heterosis is shaped by evolution to achieve robustness against noise and genetic mixing by sexual recombination in gene expression dynamics. That is, a mixed population is more robust to noise than inbred populations, as revealed experimentally by reduced phenotypic variance. The observed link between heterosis and phenotypic robustness, Mendelian dominance, and the convex single-humped fitness landscape represents a novel avenue in evolution and genetics research.

Predicting the quantitative regulatory function of a TF based on factors such as binding sequence, binding location and promoter type is not possible. The interconnected nature of gene networks and the difficulty in tuning individual TF concentrations makes the isolated study of TF function challenging. Here we present a library of E. coli strains designed to allow for precise control of the concentration of individual TFs enabling the study of the role of TF concentration on physiology and regulation. We demonstrate the usefulness of this resource by measuring the regulatory function of the zinc responsive TF, ZntR and the paralogous TF pair, GalR/GalS. For ZntR, we find that zinc alters ZntR regulatory function in a way that enables activation of the regulated gene to be robust with respect to ZntR concentration. For GalR and GalS, we are able to demonstrate that these parlogous TFs have fundamentally distinct regulatory roles beyond differences in binding affinity.

Adipocytes contribute to metabolic disorders such as obesity, diabetes, and atherosclerosis. Prior characterizations of the transcriptional network driving adipogenesis have overlooked transiently acting transcription factors (TFs), genes, and regulatory elements that are essential for proper differentiation. Moreover, traditional gene regulatory networks provide neither mechanistic details about individual regulatory element–gene relationships nor temporal information needed to define a regulatory hierarchy that prioritizes key regulatory factors. To address these shortcomings, we integrate kinetic chromatin accessibility (ATAC-seq) and nascent transcription (PRO-seq) data to generate temporally resolved networks that describe TF binding events and resultant effects on target gene expression. Our data indicate which TF families cooperate with and antagonize each other to regulate adipogenesis. Compartment modeling of RNA polymerase density quantifies how individual TFs mechanistically contribute to distinct steps in transcription. The glucocorticoid receptor activates transcription by inducing RNA polymerase pause release, whereas SP and AP-1 factors affect RNA polymerase initiation. We identifyTwist2as a previously unappreciated effector of adipocyte differentiation. We find that TWIST2 acts as a negative regulator of 3T3-L1 and primary preadipocyte differentiation. We confirm thatTwist2knockout mice have compromised lipid storage within subcutaneous and brown adipose tissue. Previous phenotyping ofTwist2knockout mice and Setleis syndromeTwist2−/−patients noted deficiencies in subcutaneous adipose tissue. This network inference framework is a powerful and general approach for interpreting complex biological phenomena and can be applied to a wide range of cellular processes.

Significant phenotypic differences exist between the vascular endothelium of different organs, including cell-cell junctions, paracellular fluid transport, shape, and mural cell coverage. These organ-specific morphological features ultimately manifest as different functional capacities, as demonstrated by the dramatic differences in capillary permeability between the leaky vessels of the liver compared to the almost impermeable vasculature found in the brain. While these morphological and functional differences have been long appreciated, the molecular basis of endothelial organ specialization remains unclear. To determine the epigenetic and transcriptional mechanisms driving this functional heterogeneity, we profiled accessible chromatin, as well as gene expression, in six different organs, across three distinct time points, during murine development and in adulthood. After identifying both common, and organ-specific DNA motif usage and transcriptional signatures, we then focused our studies on the endothelium of the central nervous system. Using single cell RNA-seq, we identified key gene regulatory networks governing brain blood vessel maturation, including TCF/LEF and FOX transcription factors. Critically, these unique regulatory regions and gene expression signatures are evolutionarily conserved in humans. Collectively, this work provides a valuable resource for identifying the transcriptional regulators controlling organ-specific endothelial specialization and provides novel insight into the gene regulatory networks governing the maturation and maintenance of the cerebrovasculature.

Boundary domains delimit and organize organ growth throughout plant development almost relentlessly building plant architecture and morphogenesis. Boundary domains display reduced growth and orchestrate development of adjacent tissues in a non-cell autonomous manner. How these two functions are achieved remains elusive despite the identification of several boundary-specific genes. Here, we show using morphometrics at the organ and cellular levels that leaf boundary domain development requires SPINDLY (SPY), an O-fucosyltransferase, to act as cell growth repressor. Further we show that SPY acts redundantly with the CUP-SHAPED COTYLEDON transcription factors (CUC2 and CUC3), which are major determinants of boundaries development. Accordingly at the molecular level, CUC2 and SPY repress a common set of genes involved in cell wall loosening providing a molecular framework for the growth repression associated with boundary domains. Atomic force microscopy (AFM) confirmed that young leaf boundary domain cells have stiffer cell walls than marginal outgrowth. This differential cell wall stiffness was reduced in spy mutant. Taken together our data reveal a concealed CUC2 cell wall associated gene network linking tissue patterning with cell growth and mechanics.Summary statementDecreased cell-wall loosening gene expression contributes to the coordination of cell growth and mechanics with tissue patterning thus driving boundary development.

BackgroundBiologists long recognized that the genetic information encoded in DNA leads to trait innovation via gene regulatory network (GRN) in development.ResultsHere, we generated paired expression and chromatin accessibility data during rumen and esophagus development in sheep and revealed 1,601 active ruminant-specific conserved non-coding elements (active-RSCNEs). To interpret the function of these active-RSCNEs, we developed a Conserved Non-coding Element interpretation method by gene Regulatory network (CNEReg) to define toolkit transcription factors (TTF) and model its regulation on rumen specific gene via batteries of active-RSCNEs during development. Our developmental GRN reveals 18 TTFs and 313 active-RSCNEs regulating the functional modules of the rumen and identifies OTX1, SOX21, HOXC8, SOX2, TP63, PPARG and 16 active-RSCNEs that functionally distinguish the rumen from the esophagus.ConclusionsWe argue that CNEReg is an attractive systematic approach to integrate evo-devo concepts with omics data to understand how gene regulation evolves and shapes complex traits.

Characterizing the underlying topology of gene regulatory networks is one of the fundamental problems of systems biology. Ongoing developments in high throughput sequencing technologies has made it possible to capture the expression of thousands of genes at the single cell resolution. However, inherent cellular heterogeneity and high sparsity of the single cell datasets render void the application of regular Gaussian assumptions for constructing gene regulatory networks. Additionally, most algorithms aimed at single cell gene regulatory network reconstruction, estimate a single network ignoring group-level (cell-type) information present within the datasets. To better characterize single cell gene regulatory networks under different but related conditions we propose the joint estimation of multiple networks using multiview graph learning (mvGL). The proposed method is developed based on recent works in graph signal processing (GSP) for graph learning, where graph signals are assumed to be smooth over the unknown graph structure. Graphs corresponding to the different datasets are regularized to be similar to each other through a learned consensus graph. We further kernelize mvGL with the kernel selected to suit the structure of single cell data. An efficient algorithm based on prox-linear block coordinate descent is used to optimize mvGL. We study the performance of mvGL using synthetic data generated with a diverse set of parameters. We further show that mvGL successfully identifies well-established regulators in a mouse embryonic stem cell differentiation study and a cancer clinical study of medulloblastoma.

BackgroundA central problem in developmental and synthetic biology is understanding the mechanisms by which cells in a tissue or a Petri dish process external cues and transform such information into a coherent response, e.g., a terminal differentiation state. It was long believed that this type of positional information could be entirely attributed to a gradient of concentration of a specific signaling molecule (i.e., a morphogen). However, advances in experimental methodologies and computer modeling have demonstrated the crucial role of the dynamics of a cell’s gene regulatory network (GRN) in decoding the information carried by the morphogen, which is eventually translated into a spatial pattern. This morphogen interpretation mechanism has gained much attention in systems biology as a tractable system to investigate the emergent properties of complex genotype-phenotype maps.MethodsIn this study, we apply a Markov chain Monte Carlo (MCMC)-like algorithm to probe the design space of three-node GRNs with the ability to generate a band-like expression pattern (target phenotype) in the middle of an arrangement of 30 cells, which resemble a simple (1-D) morphogenetic field in a developing embryo. Unlike most modeling studies published so far, here we explore the space of GRN topologies with nodes having the potential to perceive the same input signal differently. This allows for a lot more flexibility during the search space process, and thus enables us to identify a larger set of potentially interesting and realizable morphogen interpretation mechanisms.ResultsOut of 2061 GRNs selected using the search space algorithm, we found 714 classes of network topologies that could correctly interpret the morphogen. Notably, the main network motif that generated the target phenotype in response to the input signal was the type 3 Incoherent Feed-Forward Loop (I3-FFL), which agrees with previous theoretical expectations and experimental observations. Particularly, compared to a previously reported pattern forming GRN topologies, we have uncovered a great variety of novel network designs, some of which might be worth inquiring through synthetic biology methodologies to test for the ability of network design with minimal regulatory complexity to interpret a developmental cue robustly.Author summarySystems biology is a fast growing field largely powered by advances in high-performance computing and sophisticated mathematical modeling of biological systems. Based on these advances, we are now in a position to mechanistically understand and accurately predict the behavior of complex biological processes, including cell differentiation and spatial pattern formation during embryogenesis. In this article, we use an in silico approach to probe the design space of multi-input, three-node Gene Regulatory Networks (GRNs) capable of generating a striped gene expression pattern in the context of a simplified 1-D morphogenetic field.

We present DeepMAPS (Deep learning-based Multi-omics Analysis Platform for Single-cell data) for biological network inference from single-cell multi-omics (scMulti-omics). DeepMAPS includes both cells and genes in a heterogeneous graph to simultaneously infer cell-cell, cell-gene, and gene-gene relations. The multi-head attention mechanism in a graph transformer considers the heterogeneous relation among cells and genes within both local and global context, making DeepMAPS robust to data noise and scale. We benchmarked DeepMAPS on 18 scMulti-omics datasets for cell clustering and biological network inference, and the results showed that our method outperformed various existing tools. We further applied DeepMAPS on lung tumor leukocyte CITE-seq data and matched diffuse small lymphocytic lymphoma scRNA-seq and scATAC-seq data. In both cases, DeepMAPS showed competitive performance in cell clustering and predicted biologically meaningful cell-cell communication pathways based on the inferred gene networks. Note that we deployed a webserver using DeepMAPS implementation equipped with multiple functions and visualizations to improve the feasibility and reproducibility of scMulti-omics data analysis. Overall, DeepMAPS represents a heterogeneous graph transformer for single-cell study and may benefit the use of scMulti-omics data in various biological systems.

Expression quantitative trait loci (eQTLs), or single nucleotide polymorphisms (SNPs) that affect average gene expression levels, provide important insights into context-specific gene regulation. Classic eQTL analyses use one-to-one association tests, which test gene-variant pairs individually and ignore correlations induced by gene regulatory networks and linkage disequilibrium. Probabilistic topic models, such as latent Dirichlet allocation, estimate latent topics for a collection of count observations. Prior multi-modal frameworks that bridge genotype and expression data assume matched sample numbers between modalities. However, many data sets have a nested structure where one individual has several associated gene expression samples and a single germline genotype vector. Here, we build a telescoping bimodal latent Dirichlet allocation (TBLDA) framework to learn shared topics across gene expression and genotype data that allows multiple RNA-sequencing samples to correspond to a single individual9s genotype. By using raw count data, our model avoids possible adulteration via normalization procedures. Ancestral structure is captured in a genotype-specific latent space, effectively removing it from shared components. Using GTEx v8 expression data across ten tissues and genotype data, we show that the estimated topics capture meaningful and robust biological signal in both modalities, and identify associations within and across tissue types. We identify 53,358 cis-eQTLs and 1,173 trans-eQTLs by conducting eQTL mapping between the most informative features in each topic. Our TBLDA model is able to identify associations using raw sequencing count data when the samples in two separate data modalities are matched one-to-many, as is often the case in biological data. All software is available at: https://github.com/gewirtz/TBLDA

Recent advances in single-cell sequencing technologies have enabled simultaneous measurement of multiple cellular modalities, including various combinations of transcriptome, genome and epigenome. However, comprehensive profiling of the histone post-translational modifications that influence gene expression at single-cell resolution has remained limited. Here, we introduce EpiDamID, an experimental approach to target a diverse set of chromatin types by leveraging the binding specificities of genetically engineered proteins. By fusing Dam to single-chain variable fragment antibodies, engineered chromatin reader domains, or endogenous chromatin-binding proteins, we render the DamID technology and all its implementations compatible with the genome-wide identification of histone post-translational modifications. Importantly, this enables the joint analysis of chromatin marks and transcriptome in a variety of biological systems at the single-cell level. In this study, we use EpiDamID to profile single-cell Polycomb occupancy in mouse embryoid bodies and provide evidence for hierarchical gene regulatory networks. We further demonstrate the applicability of this method to in vivo systems by mapping H3K9me3 in early zebrafish embryogenesis, and detect striking heterochromatic regions specifically in the notochord. Overall, EpiDamID is a new addition to a vast existing toolbox for obtaining systematic insights into the role of chromatin states during dynamic cellular processes.

Genetic variants identified through genome-wide association studies (GWAS) are typically non-coding and exert small regulatory effects on downstream genes, but which downstream genes are ultimately impacted and how they confer risk remains mostly unclear. Conversely, variants that cause rare Mendelian diseases are often coding and have a more direct impact on disease development. We demonstrate that common and rare genetic diseases can be linked by studying the gene regulatory networks impacted by common disease-associated variants. We implemented this in the ‘Downstreamer’ method and applied it to 44 GWAS traits and find that predicted downstream “key genes” are enriched with Mendelian disease genes, e.g. key genes for height are enriched for genes that cause skeletal abnormalities and Ehlers-Danlos syndromes. We find that 82% of these key genes are located outside of GWAS loci, suggesting that they result from complex trans regulation rather than being impacted by disease-associated variants in cis. Finally, we discuss the challenges in reconstructing gene regulatory networks and provide a roadmap to improve identification of these highly connected genes for common traits and diseases.

SUMMARYDrosophila has long been a successful model organism in multiple fields such as genetics and developmental biology. Drosophila genome is relatively smaller and less redundant, yet largely conserved with mammals, making it a productive model in studies of embryogenesis, cell signaling, disease mechanisms, etc. Spatial gene expression pattern is critical for understanding of complex signaling pathways and cell-cell interactions, whereas temporal gene expression changes need to be tracked during highly dynamic activities such as tissue development and disease progression. Systematic studies in Drosophila as a whole are still impeded by lack of these spatiotemporal transcriptomic information. Drosophila embryos and tissues are of relatively small size, limiting the application of current technologies to comprehensively resolve their spatiotemporal gene expression patterns. Here, utilizing SpaTial Enhanced REsolution Omics-sequencing (Stereo-seq), we dissected the spatiotemporal transcriptomic changes of developing Drosophila with high resolution and sensitivity. Our data recapitulated the spatial transcriptomes of embryonic and larval development in Drosophila. With these data, we identified known and previously undetected subregions in several tissues during development, and revealed known and potential gene regulatory networks of transcription factors within their topographic background. We further demonstrated that Stereo-seq data can be used for 3D reconstruction of Drosophila embryo spatial transcriptomes. Our data provides Drosophila research community with useful resources of spatiotemporally resolved transcriptomic information across developmental stages.

Reliable operation of cellular programs depends crucially on the specificity of biomolecular interactions. In gene regulatory networks, the appropriate expression of genes is determined through the specific binding of transcription factors (TFs) to their cognate DNA sequences. However, the large genomic background likely contains many DNA sequences showing similarity to TF target motifs, potentially allowing for substantial non-cognate TF binding with low specificity. Whether and how non-cognate TF binding impacts cellular function and fitness remains unclear. We show that increased expression of different transcriptional regulators in Escherichia coli and Salmonella enterica can significantly inhibit population growth across multiple environments. This effect depends upon (i) TF binding to a large number of DNA sequences with low specificity, (ii) TF cooperativity, and (iii) the ratio of TF to DNA. DNA binding due to the limited specificity of promiscuous or non-native TFs can thus severely impact fitness, giving rise to a fundamental biophysical constraint on gene regulatory design and evolution.

Biotic stresses are environmental factors that cause a variety of crop diseases and damages. In contrast, crops trigger specific transduction signaling pathways that the hormones are the central players. Integrative OMICS for systems genetic engineering approach contributes in the understanding of molecular mechanisms. In this research, the system biology approaches were applied to discover particular molecular interactions between biotic stresses and hormonal signaling in barley. The meta-analysis of the data identified a total of 1232 and 304 differentially expressed genes (DEGs) respectively so that were significantly involved in defense processes and hormone signaling. A total of 24 TFs belonged to 15 conserved families and 6 TFs belonged to 6 conserved families were identified for biotic and hormonal data respectively, whereas NF-YC, GNAT, and whirly families were the most abundant groups. The functional analysis of the upstream regions for over-represented cis-acting elements revealed that were involved activation of transcription factors in response to pathogens and hormones. Based on the co-expression analysis, 6 and 7 distinct co-expression modules related to biotic stresses and hormonal signaling were respectively uncovered. The gene network analysis also identified novel hub genes such as TIM10, DRT101, ADG1, and TRA2 which may be involved in regulating defense responses to biotic stresses. In addition, many new genes with unknown function were obtained. Since this study represents a first preliminary curated system biology analysis of barley transcriptomic responses to biotic stresses and hormone treatments, introduces important candidate genes that may be beneficial to crop biotechnologists to accelerate genetic engineering programs.

SummaryMacrophages are driven to form distinct functional phenotypes in response to different immunological stimuli, in a process widely referred to as macrophage polarization. Transcriptional regulators that guide macrophage polarization in response to a given trigger remain largely unknown. In this study, we interrogate the programmable landscape in macrophages to find regulatory panels that determine the precise polarization state that a macrophage is driven to. Towards this, we configure an integrative network analysis pipeline that utilizes macrophage transcriptomes in response to 28 different stimuli and reconstructs contextualized human gene regulatory networks, and identifies epicentres of perturbations in each case. We find that these contextualized regulatory networks form a spectrum of thirteen distinct clusters with M1 and M2 at the two ends. Using our computational pipeline, we identify combinatorial panels of epicentric transcription factors (TFs) for each polarization state. We demonstrate that a set of three TFs i.e., NFE2L2, BCL3 and CEBPB, is sufficient to change the polarization destination from M1 to M2. siRNA-mediated knockdown of the 3-TF set in THP1 derived M0 cells, despite exposure to an M1 stimulant, significantly attenuated the shift to M1 phenotype, and instead switched to an M2-like phenotype. Further, the siRNA-mediated knockdown of the 3-TF set rendered the macrophages hyper-susceptible to S. aureus infection, clearly showing the effect of the knockdown.

The medical and scientific response to emerging pathogens is often severely hampered by ignorance of the genetic determinants of virulence, drug resistance, and clinical outcomes that could be used to identify therapeutic drug targets and forecast patient trajectories1–5. Taking the newly emergent multidrug-resistant bacteriaMycobacterium abscessusas an example6, we show that combining high dimensional phenotyping with whole genome sequencing in a phenogenomic analysis can rapidly reveal actionable systems-level insights into bacterial pathobiology. Usingin vitroandin vivophenotyping, we discovered three distinct clusters of isolates, each associated with a different clinical outcome. We combined genome-wide association studies (GWAS) with proteome-wide computational structural modelling7to define likely causal variants, and employed direct coupling analysis (DCA)8to identify co-evolving, and therefore potentially epistatic, gene networks. We then usedin vivoCRISPR-based silencing to validate our findings, defining a novel secretion system controlling virulence inM. abscessus, and illustrating how phenogenomics can reveal critical pathways within emerging pathogenic bacteria.

Gene expression, the production of protein from DNA and mRNA in the biological cell, is inherently stochastic. Cells with identical DNA exhibit fluctuations or ‘noise’ in gene expression. This noise propagates over gene regulatory networks (GRNs), which encode gene-gene interactions. The propagated ‘extrinsic’ noise interacts and combines with ‘intrinsic’ noise to affect biological decisions. Consequently, it is essential to understand how GRN topology affects total noise. Recently, uncertainty principles were established for noise propagation over GRN. In particular, in ring GRNs, exactly one node can have noise reduction below the intrinsic limit. We establish necessary and sufficient conditions for noise reduction in ring GRN. Specifically, for two- and three-node rings, an odd number of negative regulations is necessary for noise reduction. Further, sufficiency is ensured if sensitivities to input for feedforward and feedback regulations are bounded from below and above, respectively. These constraints are valid even if the ring GRN are regulated by an upstream gene. Finally, we use graph theory to decompose noise propagation in a general directed network over its strongly connected components. The combination of graph theory and stochastic processes may be a general framework for studying noise propagation.

In stem cell biology, cellular pluripotency describes the capacity of a given cell to differentiate into multiple cell types. From a statistical physics perspective, entropy provides a statistical measure of randomness and has been demonstrated as a way to quantitate pluripotency when considering biological gene networks. Furthermore, recent theoretical work has established a relationship between Ricci curvature (a geometric measure of “flatness”) and entropy (also related to robustness), which one can exploit to link the geometric quantity of curvature to the statistical quantity of entropy. Therefore, this study seeks to explore Ricci curvature in biological gene networks as a descriptor of pluripotency and robustness among gene pathways. Here, we investigate Forman-Ricci curvature, a combinatorial discretization of Ricci curvature, along with network entropy, to explore the relationship of the two quantities as they occur in gene networks. First, we demonstrate our approach on an experiment of stem cell gene expression data. As expected, we find Ricci curvature directly correlates with network entropy, suggesting Ricci curvature could serve as an indicator for cellular pluripotency much like entropy. Second, we measure Forman-Ricci curvature in a dataset of cancer and non-cancer cells from melanoma patients. We again find Ricci curvature is increased in the cancer state, reflecting increased pluripotency or “stemness”. Further, we locally examine curvature on the gene level to identify several genes and gene pathways with known relevance to melanoma. In turn, we conclude Forman-Ricci curvature provides valuable biological information related to pluripotency and pathway functionality. In particular, the advantages of this geometric approach are promising for extension to higher-order topological structures in order to represent more complex features of biological systems.

Despite its high and direct impact on nearly all biological processes, the underlying structure of gene-gene interaction networks is investigated so far according to pair connections. To address this, we explore the gene interaction networks of the yeast Saccharomyces cerevisiae beyond pairwise interaction using the structural balance theory (SBT). Specifically, we ask whether essential and nonessential gene interaction networks are structurally balanced. We study triadic interactions in the weighted signed undirected gene networks and observe that balanced and unbalanced triads are over and underrepresented in both networks, thus beautifully in line with the strong notion of balance. Moreover, we note that the energy distribution of triads is significantly different in both essential and nonessential networks compared with the shuffled networks. Yet, this difference is greater in the essential network regarding the frequency as well as the energy of triads. Additionally, results demonstrate that triads in the essential gene network are more interconnected through sharing common links, while in the nonessential network they tend to be isolated. Last but not least, we investigate the contribution of all-length signed walks and its impact on the degree of balance. Our findings reveal that interestingly when considering longer cycles the nonessential gene network is more balanced compared to the essential network.

Evolution optimizes cellular behavior throughout sequential generations by selecting the successful individual cells in a given context. As gene regulatory networks (GRNs) determine the behavior of single cells by ruling the activation of different processes - such as cell differentiation and death - how GRNs change from one generation to the other might have a relevant impact on the course of evolution. It is not clear, however, which mechanisms that affect GRNs effectively favor evolution and how. Here, we use a population of computational robotic models controlled by artificial gene regulatory networks (AGRNs) to evaluate the impact of different genetic modification strategies in the course of evolution. The virtual agent senses the ambient and acts on it as a bacteria in different phototaxis-like tasks - orientation to light, phototaxis, and phototaxis with obstacles. We studied how the strategies of gradual and abrupt changes on the AGRNs impact evolution considering multiple levels of task complexity. The results indicated that a gradual increase in the complexity of the performed tasks is beneficial for the evolution of the model. Furthermore, we have seen that larger gene regulatory networks are needed for more complex tasks, with single-gene duplication being an excellent evolutionary strategy for growing these networks, as opposed to full-genome duplication. Studying how GRNs evolved in a biological environment allows us to improve the computational models produced and provide insights into aspects and events that influenced the development of life on earth.

The aging eye experiences physiological changes that include decreased visual function and increased risk of retinal degeneration. Although there are transcriptomic signatures in the aging retina that correlate with these physiological changes, the gene regulatory mechanisms that contribute to cellular homeostasis during aging remain to be determined. Here, we integrated ATAC-seq and RNA-seq data to identify 61 transcription factors that showed differential activity in agingDrosophilaphotoreceptors. These 61 age-regulated transcription factors include two circadian regulators, Clock and cycle, that showed sustained increases in activity during aging. When we disrupted Clock activity in adult photoreceptors, we observed changes in expression of 15 – 20% of genes including key components of the phototransduction machinery and many eye-specific transcription factors. Using ATAC-seq, we showed that loss of Clock activity leads to changes in activity of 31 transcription factors and causes a progressive decrease in global levels of chromatin accessibility in photoreceptors. Supporting a key role for Clock-dependent transcription in the eye, disruption of Clock activity in photoreceptors also induced light-dependent retinal degeneration and increased oxidative stress, independent of light exposure. Together, our data suggests that the circadian regulators Clock and cycle act as neuroprotective factors in the aging eye by directing gene regulatory networks that maintain expression of the phototransduction machinery and counteract oxidative stress.

Severe acute respiratory syndrome corona virus 2 (SARS-CoV-2)/ coronavirus disease 2019 (COVID-19) infection is the leading cause of respiratory tract infection associated mortality worldwide. The aim of the current investigation was to identify the differentially expressed genes (DEGs) and enriched pathways in COVID-19 infection and its associated complications by bioinformatics analysis, and to provide potential targets for diagnosis and treatment. Valid next-generation sequencing (NGS) data of 93 COVID 19 samples and 100 non COVID 19 samples (GSE156063) were obtained from the Gene Expression Omnibus database. Gene ontology (GO) and REACTOME pathway enrichment analysis was conducted to identify the biological role of DEGs. In addition, a protein-protein interaction network, modules, miRNA-hub gene regulatory network, TF-hub gene regulatory network and receiver operating characteristic curve (ROC) analysis were used to identify the key genes. A total of 738 DEGs were identified, including 415 up regulated genes and 323 down regulated genes. Most of the DEGs were significantly enriched in immune system process, cell communication, immune system and signaling by NTRK1 (TRKA). Through PPI, modules, miRNA-hub gene regulatory network, TF-hub gene regulatory network analysis, ESR1, UBD, FYN, STAT1, ISG15, EGR1, ARRB2, UBE2D1, PRKDC and FOS were selected as hub genes, which were expressed in COVID-19 samples relative to those in non COVID-19 samples, respectively. Among them, ESR1, UBD, FYN, STAT1, ISG15, EGR1, ARRB2, UBE2D1, PRKDC and FOS were suggested to be diagonstic factors for COVID-19. The findings from this bioinformatics analysis study identified molecular mechanisms and the key hub genes that might contribute to COVID-19 infection and its associated complications.

Understanding changes in gene expression under the effects of a perturbation is a key goal of systems biology. A powerful approach to address this goal uses gene networks and describes the perturbation’s effects as a rewiring of each gene’s connections. This approach is known as differential network (DN) analysis. Here, we used DNs to analyze RNA-sequencing time series datasets, focusing on expression changes: (i) In the saliva of a human subject after vaccination with a pneumococcal vaccine (PPSV23), and (ii) in B cells treatedex vivowith a monoclonal antibody drug (Rituximab). Using network community detection, we revealed the collective behavior of clusters of genes, and detected communities of genes based on their longitudinal behavior, and corresponding pathway activations. We identified biological pathways consistent with the mechanism of action of the vaccine and with Rituximab’s targets. The approach may be useful in drug development by providing an effective analysis of expressing changes in response to a drug.

SUMMARYGlioblastoma (GBM) is a complex disease that is difficult to treat. Establishing the complex genetic interactions regulating cell fate decisions in GBM can help to shed light on disease aggressivity and improved treatments. Networks and data science offer alternative approaches to classical bioinformatics pipelines to study gene expression patterns from single-cell datasets, helping to distinguish genes associated with control of differentiation and thus aggressivity. Here, we applied a host of data theoretic techniques, including clustering algorithms, Waddington landscape reconstruction, trajectory inference algorithms, and network approaches, to compare gene expression patterns between pediatric and adult GBM, and those of adult glioma-derived stem cells (GSCs) to identify the key molecular regulators of the complex networks driving GBM/GSC and predict their cell fate dynamics. Using these tools, we identified critical genes and transcription factors coordinating cell state transitions from stem-like to mature GBM phenotypes, including eight transcription factors (OLIG1/2, TAZ, GATA2, FOXG1, SOX6, SATB2, YY1) and four signaling genes (ATL3, MTSS1, EMP1, and TPT1) as clinically targetable novel putative function interactions differentiating pediatric and adult GBMs from adult GSCs. Our study provides strong evidence of the applicability of complex systems approaches for reverse-engineering gene networks from patient-derived single-cell datasets and inferring their complex dynamics, bolstering the search for new clinically relevant targets in GBM.

The properties of random Boolean networks as models of gene regulation have been investigated extensively by the statistical physics community. In the past two decades, there has been a dramatic increase in the reconstruction and analysis of Boolean models of biological networks. In such models, neither network topology nor Boolean functions (or logical update rules) should be expected to be random. In this contribution, we focus on biologically meaningful types of Boolean functions, and perform a systematic study of their preponderance in gene regulatory networks. By applying the k[P] classification based on number of inputs k and bias P of functions, we find that most Boolean functions astonishingly have odd bias in a reference biological dataset of 2687 functions compiled from published models. Subsequently, we are able to explain this observation along with the enrichment of read-once functions (RoFs) and its subset, nested canalyzing functions (NCFs), in the reference dataset in terms of two complexity measures: Boolean complexity based on string lengths in formal logic which is yet unexplored in the biological context, and the average sensitivity. Minimizing the Boolean complexity naturally sifts out a subset of odd-biased Boolean functions which happen to be the RoFs. Finally, we provide an analytical proof that NCFs minimize not only the Boolean complexity, but also the average sensitivity in their k[P] set.

Despite holding great therapeutic potential, existing protocols for in vitro chondrogenesis and hyaline cartilage production from human induced pluripotent stem cells (hiPSC) are laborious and complex with unclear long-term consequences. Here, we developed a simple xeno- and feeder-free protocol for human hyaline cartilage production in vitro using hydrogel-cultured multi-tissue organoids (MTOs). We investigate gene regulatory networks during spontaneous hiPSC-MTO differentiation using RNA sequencing and bioinformatic analyses. We find the interplays between BMPs and neural FGF pathways are associated with the phenotype transition of MTOs. We recognize TGF-beta/BMP and Wnt signaling likely contribute to the long-term maintenance of MTO cartilage growth and further adoption of articular cartilage development. By comparing the MTO transcriptome with human lower limb chondrocytes, we observe that the expression of chondrocyte-specific genes in MTO shows a strong correlation with fetal lower limb chondrocytes. Collectively, our findings describe the self-organized emergence of hyaline cartilage in MTO, its associated molecular pathways, and its spontaneous adoption of articular cartilage development trajectory.

Deep neural networks implementing generative models for dimensionality reduction have been extensively used for the visualization and analysis of genomic data. One of their key limitations is lack of interpretability: it is challenging to quantitatively identify which input features are used to construct the embedding dimensions, thus preventing insight into why cells are organized in a particular data visualization, for example. Here we present a scalable, interpretable variational autoencoder (siVAE) that is interpretable by design: it learns feature embeddings that guide the interpretation of the cell embeddings in a manner analogous to factor loadings of factor analysis. siVAE is as powerful and nearly as fast to train as the standard VAE but achieves full interpretability of the embedding dimensions. Using siVAE, we exploit a number of connections between dimensionality reduction and gene network inference to identify gene neighborhoods and gene hubs, without the explicit need for gene network inference. We observe a systematic difference in the gene neighborhoods identified by dimensionality reduction methods and gene network inference algorithms in general, suggesting they provide complementary information about the underlying structure of the gene co-expression network. Finally, we apply siVAE to implicitly learn gene networks for individual iPSC lines and uncover a correlation between neuronal differentiation efficiency and loss of co-expression of several mitochondrial complexes, including NADH dehydrogenase, cytochrome C oxidase, and cytochrome b.

CRISPR-based programmable transcriptional activators (PTAs) are used in plants for rewiring gene networks. Better tuning of their activity in a time and dose-dependent manner should allow precise control of gene expression. Here, we report the optimization of a Copper Inducible system called CI-switch for conditional gene activation in Nicotiana benthamiana. In the presence of copper, the copper-responsive factor CUP2 undergoes a conformational change and binds a DNA motif named copper-binding site (CBS). In this study, we tested several activation domains fused to CUP2 and found that the non-viral Gal4 domain results in strong activation of a reporter gene equipped with a minimal promoter, offering advantages over previous designs. To connect copper regulation with downstream programable elements, several copper-dependent configurations of the strong dCasEV2.1 PTA were assayed, aiming at maximizing activation range, while minimizing undesired background expression. The best configuration involved a dual copper regulation of the two protein components of the PTA, namely dCas9:EDLL and MS2:VPR, and a constitutive RNA pol III-driven expression of the third component, a guide RNA with anchoring sites for the MS2 RNA-binding domain. With these optimizations in place, the CI/dCasEV2.1 system resulted in copper-dependent activation rates of 2,600-fold for the endogenous N. benthamiana DFR gene, with negligible expression in the absence of the trigger. The tight regulation of copper over CI/dCasEV2.1 makes this system ideal for the conditional production of plant-derived metabolites and recombinant proteins in the field.

Induction of meiosis requires exogenous signals that activate internal gene regulatory networks. Meiotic commitment ensures the irreversible continuation of meiosis, even upon withdrawal of the meiosis inducing signals. Budding yeast cells have a unique property in that cells enter meiosis when starved, but if given nutrient-rich medium prior to the commitment point, they exit meiosis and enter mitosis. After the meiotic commitment point in prometaphase I, cells remain in meiosis even with addition of nutrients. Despite the importance of meiotic commitment in ensuring the production of gametes, only a few genes involved in the commitment process are known. We performed a genome-scale screen in budding yeast and discovered new regulators of meiotic commitment including Bcy1, which is involved in nutrient sensing, the meiosis-specific kinase Ime2, Polo kinase Cdc5, and the 14-3-3 proteins Bmh1 and Bmh2. Importantly, we found that Bmh1 and Bmh2 are involved in multiple processes throughout meiosis including the maintenance of the middle meiosis transcription factor Ndt80, activation of Cdc5, and interaction with an RNA-binding protein Pes4, which is important for regulating the timing of translation of several mRNAs in meiosis II. This study identifies a meiotic commitment regulatory network with the 14-3-3 proteins functioning as central regulators.

Single cell sequencing has recently allowed the generation of exhaustive root cell atlases. However, some cell types are elusive and remain underrepresented. Here, we use a second- generation single cell approach, where we zoom in on the root transcriptome sorting with specific markers to profile the phloem poles at an unprecedented resolution. Our data highlight the similarities among the developmental trajectories and gene regulatory networks communal to protophloem sieve element (PSE) adjacent lineages in relation to PSE enucleation, a key event in phloem biology.As a signature for early PSE-adjacent lineages, we have identified a set of DNA-binding with one finger (DOF) transcription factors, the PINEAPPLEs (PAPL), that act downstream of PHLOEM EARLY DOF (PEAR) genes, and are important to guarantee a proper root nutrition in the transition to autotrophy.Our data provide a holistic view of the phloem poles that act as a functional unit in root development.

Understanding the impact of regulatory variants on complex phenotypes is a significant challenge because the genes and pathways that are targeted by such variants are typically unknown. Furthermore, a regulatory variant might influence a particular gene’s expression in a cell type or tissue-specific manner. Cell-type specific long-range regulatory interactions that occur between a distal regulatory sequence and a gene offers a powerful framework for understanding the impact of regulatory variants on complex phenotypes. However, high-resolution maps of such long-range interactions are available only for a handful of model cell lines. To address this challenge, we have developed L-HiC-Reg, a Random Forests based regression method to predict high- resolution contact counts in new cell lines, and a network-based framework to identify candidate cell line-specific gene networks targeted by a set of variants from a Genome-wide association study (GWAS). We applied our approach to predict interactions in 55 Roadmap Epigenome Consortium cell lines, which we used to interpret regulatory SNPs in the NHGRI GWAS catalogue. Using our approach, we performed an in-depth characterization of fifteen different phenotypes including Schizophrenia, Coronary Artery Disease (CAD) and Crohn’s disease. In CAD, we found differentially wired subnetworks consisting of known as well as novel gene targets of regulatory SNPs. Taken together, our compendium of interactions and associated network-based analysis pipeline offers a powerful resource to leverage long-range regulatory interactions to examine the context-specific impact of regulatory variation in complex phenotypes.

Time-course single-cell RNA sequencing (scRNA-seq) data have been widely applied to reconstruct the cell-type-specific gene regulatory networks by exploring the dynamic changes of gene expression between transcription factors (TFs) and their target genes. The existing algorithms were commonly designed to analyze bulk gene expression data and could not deal with the dropouts and cell heterogeneity in scRNA-seq data. In this paper, we developed dynDeepDRIM that represents gene pair joint expression as images and considers the neighborhood context to eliminate the transitive interactions. dynDeepDRIM integrated the primary image, neighbor images with time-course into a four-dimensional tensor and trained a convolutional neural network to predict the direct regulatory interactions between TFs and genes. We evaluated the performance of dynDeepDRIM on five time-course gene expression datasets. dynDeepDRIM outperformed the state-of-the-art methods for predicting TF-gene direct interactions and gene functions. We also observed gene functions could be better performed if more neighbor images were involved.

The gene regulatory network (GRN) that underlies echinoderm skeletogenesis is a prominent model of GRN architecture and evolution. KirrelL is an essential downstream effector gene in this network and encodes an Ig-superfamily protein required for the fusion of skeletogenic cells and the formation of the skeleton. In this study, we dissected the transcriptional control region of the kirrelL gene of the purple sea urchin, Strongylocentrotus purpuratus. Using plasmid- and BAC-based transgenic reporter assays, we identified key cis-regulatory elements (CREs) and transcription factor inputs that regulate Sp-kirrelL, including direct, positive inputs from two key transcription factors in the skeletogenic GRN, Alx1 and Ets1. We next identified kirrelL cis-regulatory regions from seven other echinoderm species that together represent all classes within the phylum. By introducing these heterologous regulatory regions into developing sea urchin embryos we provide evidence of their remarkable conservation across ~500 million years of evolution. We dissected in detail the kirrelL regulatory region of the sea star, Patiria miniata, and demonstrated that it also receives direct inputs from Alx1 and Ets1. Our findings identify kirrelL as a component of the ancestral echinoderm skeletogenic GRN. They support the view that GRN sub-circuits, including specific transcription factor-CRE interactions, can remain stable over vast periods of evolutionary history. Lastly, our analysis of kirrelL establishes direct linkages between a developmental GRN and an effector gene that controls a key morphogenetic cell behavior, cell-cell fusion, providing a paradigm for extending the explanatory power of GRNs.

Transcriptional regulatory network (TRN) orchestrates spatio-temporal expression of genes to generate cellular responses for survival. The transcription factors (TF) regulating expression of their target genes (TG) are the fundamental units of TRN. Several databases have been developed to catalogue human TRN based on low- and high-throughput experimental and computational studies considering their importance in understanding cellular physiology. However, literature lacks comparative assessment on the strength and weakness of each database. In this study, we compared over 2.2 million regulatory pairs between 1,379 TF and 22,518 TG assembled from 14 data resources. Our study reveals that the TF and TG were common across data resources but not their regulatory pairs. We observed that the TF and TG of the regulatory pairs showed weak expression correlation, significant gene ontology overlap, co-citations in PubMed and low numbers of TF-TG pairs representing transcriptional repression relationships. Furthermore, each TF-TG regulatory pair assigned a combined confidence score reflecting its reliability based on its presence in multiple databases and co-expression. The TRN containing 2,246,598 TF-TG pairs, of which, 44,284 with the information on TF’s activating or repressing effects on their TG is available upon request. This study brings the information about transcriptional regulation scattered over the literature and databases at one place in the form of one of the most comprehensive and complete human TRNs assembled to date, which will be a valuable resource for benchmarking TRN prediction tools, and to the scientific community working in functional genomics, gene expression and gene regulation analysis.

Myocardial infarction (MI) is the leading cardiovascular diseases in worldwide, yet relatively little is known about the genes and signaling pathways involved in MI progression. The present investigation aimed to elucidate potential crucial candidate genes and pathways in MI. expression profiling by high throughput sequencing dataset (GSE132143) was downloaded from the Gene Expression Omnibus (GEO) database, which included data from 20 MI samples and 12 normal control samples. Differentially expressed genes (DEGs) were identified using t-tests in the DESeq2 R package. These DEGs were subsequently investigated by Gene Ontology (GO) and pathway enrichment analysis, a protein-protein interaction (PPI) network, modules, miRNA-hub gene regulatory network and TF-hub gene regulatory network were constructed and analyzed. Hub genes were validated by receiver operating characteristic curve (ROC) analysis. In total, 958 DEGs were identified, of which 480 were up regulated and 478 were down regulated. GO and pathway enrichment analysis results revealed that the DEGs were mainly enriched in, immune system, neuronal system, response to stimulus, and multicellular organismal process. A PPI network, modules, miRNA-hub gene regulatory network and TF-hub gene regulatory network was constructed by using Cytoscape software, and CFTR, CDK1, RPS13, RPS15A, RPS27, NOTCH1, MRPL12, NOS2, CCDC85B and ATN1 were identified as the hub genes. Our results highlight the important roles of the genes including CFTR, CDK1, RPS13, RPS15A, RPS27, NOTCH1, MRPL12, NOS2, CCDC85B and ATN1 in MI pathogenesis or therapeutic management.

The reproducibility of embryonic development is a remarkable feat of biological organization, but the underlying mechanisms are poorly understood. Clearly, gene regulatory networks are central to the orderly progression of development, but noisy molecular and cellular processes should reduce reproducibility. Here, we identify ergodicity, a type of dynamical stability, as underlying the reproducibility of development. In ergodic systems, a single timepoint measurement equals a time average. Focusing on the zebrafish tailbud, we define gene expression and cell motion states using a parallel statistical analyses of single cell RNA sequencing data and in vivo timelapse cell tracking data and a change point detection algorithm. Strikingly, the cell motion state transitions in each embryo exhibit the same patterns for both a single timepoint and a 2-3 hour time average. Both the cell motion and gene expression cell states exhibit balanced influx and outflux rates reflecting a spatiotemporal stability. Stated simply, these data indicate the pattern of changes in the tailbud doesn’t change. This ergodic pattern of cell state transitions may represent an emergent meta-state that links gene networks to the reproducible progression of embryogenesis.

Development is driven by gene regulatory networks (GRNs) that progressively dictate specification and differentiation of cell fates. The architecture of GRNs directly determines the specificity and accuracy of developmental outcomes. We report here that the core regulatory circuitry for endoderm development in C. elegans is comprised of a recursive series of interlocked feedforward modules linking a cascade of six sequentially expressed GATA-type transcription factors. This structure results in a reiterated sequential redundancy, in which removal of a single factor or alternate factors in the cascade results in no, or a mild, effect on endoderm development and gut differentiation, while elimination of any two factors that are sequentially deployed in the cascade invariably results in a strong phenotype. The strength of the observed phenotypes is successfully predicted by a computational model based on the timing and levels of transcriptional states. The feedforward regulatory logic in the GRN appears to ensure timely onset of terminal differentiation genes and allows rapid and robust lockdown of cell fate during early embryogenesis. We further found that specification-to-differentiation transition is linked through a common regulator, the END-1 GATA factor that straddles the two processes. Finally, we revealed roles for key GATA factors in establishing spatial regulatory state domains by acting as transcriptional repressors that appear to define the boundaries of the digestive tract. Our findings support a comprehensive model of the core gene network that describes how robust endoderm development is achieved during C. elegans embryogenesis.Graphic abstract

SummaryCalpains are cysteine proteases that control cell fate transitions. Although calpains are viewed as modulatory proteases displaying severe, pleiotropic phenotypes in eukaryotes, human calpain targets are also directed to the N-end rule degradatory pathway. Several of these destabilized targets are transcription factors, hinting at a gene regulatory role. Here, we analyze the gene regulatory networks of Physcomitrium patens and characterize the regulons that are deregulated in DEK1 calpain mutants. Predicted cleavage patterns of regulatory hierarchies in the five DEK1-controlled subnetworks are consistent with the gene’s pleiotropy and the regulatory role in cell fate transitions targeting a broad spectrum of functions. Network structure suggests DEK1-gated sequential transition between cell fates in 2D to 3D development. We anticipate that both our method combining phenotyping, transcriptomics and data science to dissect phenotypic traits and our model explaining the calpain’s role as a switch gatekeeping cell fate transitions will inform biology beyond plant development.

In this work we analyze the flow of information through the Fibroblast Growth Factor Receptor (FGFR) communication channel when different types of signals are transmitted by the MAPK cascade to the gene regulatory network (GRN) formed by the genes C-Myc, DUSP, and Cdc25A, which control fibroblast proliferation. We used the canonical mathematical model of the MAPK cascade coupled to a stochastic model for the activation of the gene regulatory network, subject to different types of FGF inputs (step, quadratic pulses, Dirac delta, and white noise), in order to analyze the response of the gene regulatory network to each type of signal, and determine the temporal variation of the value of its Shannon entropy in each case. Our model suggests that the sustained activation of the FGFR communication channel with a step of FGF > 1 nM is required for cell cycle progression and that during the G1/S transition the amount of uncertainty of the GRN remains at a steady value of ∼ 2.75 bits, indicating that while the fibroblast stimulation with FGF continues the G1/S transition does not require an additional interchange of information between the emitter and the gene regulatory network to be completed. We also found that either low frequency pulses of FGF or low frequency noise, both with a frequency f ≤ 2.77 Hz, are not filtered by the MAPK cascade and can modify the output of the communication channel, i.e., the amount of the effector proteins c-myc, cdc25A and DUSP. An additional effect suggested by our model is that o low frequency periodic signals and noise possibly blockage cell cycle progression because the threshold value concentration of cdc25A for the G1/S transition is not sustained in the in the nucleus during the 10 hours that this process lasts. Finally, from our model we can estimate the capacity of this communication channel in 0.96 bits/min.