Search

Alternative splicing shapes the transcriptome and contributes to each cell's unique identity, but single-cell RNA sequencing (scRNA-seq) has struggled to capture the impact of alternative splicing. We previously showed that low recovery of mRNAs from single cells led to erroneous conclusions about the cell-to-cell variability of alternative splicing. Here, we present a method, Psix, to confidently identify splicing that changes across a landscape of single cells, using a probabilistic model that is robust against the data limitations of scRNA-seq. Its autocorrelation-inspired approach finds patterns of alternative splicing that correspond to patterns of cell identity, such as cell type or developmental stage, without the need for explicit cell clustering, labeling, or trajectory inference. Applying Psix to data that follow the trajectory of mouse brain development, we identify exons whose alternative splicing patterns cluster into modules of coregulation. We show that the exons in these modules are enriched for binding by distinct neuronal splicing factors and that their changes in splicing correspond to changes in expression of these splicing factors. Thus, Psix reveals cell type–dependent splicing patterns and the wiring of the splicing regulatory networks that control them. Our new method will enable scRNA-seq analysis to go beyond transcription to understand the roles of post-transcriptional regulation in determining cell identity.

Single-cell RNA sequencing (scRNA-seq) enables molecular characterization of complex biological tissues at high resolution. The requirement of single-cell extraction, however, makes it challenging for profiling tissues such as adipose tissue, for which collection of intact single adipocytes is complicated by their fragile nature. For such tissues, single-nucleus extraction is often much more efficient and therefore single-nucleus RNA sequencing (snRNA-seq) presents an alternative to scRNA-seq. However, nuclear transcripts represent only a fraction of the transcriptome in a single cell, with snRNA-seq marked with inherent transcript enrichment and detection biases. Therefore, snRNA-seq may be inadequate for mapping important transcriptional signatures in adipose tissue. In this study, we compare the transcriptomic landscape of single nuclei isolated from preadipocytes and mature adipocytes across human white and brown adipocyte lineages, with whole-cell transcriptome. We show that snRNA-seq is capable of identifying the broad cell types present in scRNA-seq at all states of adipogenesis. However, we also explore how and why the nuclear transcriptome is biased and limited, as well as how it can be advantageous. We robustly characterize the enrichment of nuclear-localized transcripts and adipogenic regulatory lncRNAs in snRNA-seq, while also providing a detailed understanding for the preferential detection of long genes upon using this technique. To remove such technical detection biases, we propose a normalization strategy for a more accurate comparison of nuclear and cellular data. Finally, we show successful integration of scRNA-seq and snRNA-seq data sets with existing bioinformatic tools. Overall, our results illustrate the applicability of snRNA-seq for the characterization of cellular diversity in the adipose tissue.

Contemporary single-cell omics technologies have enabled complex experimental designs incorporating hundreds of samples accompanied by detailed information on sample-level conditions. Current approaches for analyzing condition-level heterogeneity in these experiments often rely on a simplification of the data such as an aggregation at the cell-type or cell-state-neighborhood level. Here we present MrVI, a deep generative model that provides sample-sample comparisons at a single-cell resolution, permitting the discovery of subtle sample-specific effects across cell populations. Additionally, the output of MrVI can be used to quantify the association between sample-level metadata and cell state variation. We benchmarked MrVI against conventional meta-analysis procedures on two synthetic datasets and one real dataset with a well-controlled experimental structure. This work introduces a novel approach to understanding sample-level heterogeneity while leveraging the full resolution of single-cell sequencing data.

RNA velocity has been rapidly adopted to guide the interpretation of transcriptional dynamics in snapshot single-cell transcriptomics data. Current approaches for estimating and analyzing RNA velocity can empirically reveal complex dynamics but lack effective strategies for quantifying the uncertainty of the estimate and its overall applicability to the system of interest. Here, we present veloVI (velocity variational inference), a deep generative modeling framework for estimating RNA velocity. veloVI learns a gene-specific dynamical model of RNA metabolism and provides a transcriptome-wide quantification of velocity uncertainty. We show in a series of examples that veloVI compares favorably to previous approaches for inferring RNA velocity with improvements in fit to the data, consistency across transcriptionally similar cells, and stability across preprocessing pipelines for quantifying RNA abundance. Further, we demonstrate that properties unique to veloVI, such as posterior velocity uncertainty, can be used to assess the appropriateness of analysis with velocity to the data at hand. Finally, we highlight veloVI as a flexible framework for modeling transcriptional dynamics by adapting the underlying dynamical model to use time-dependent transcription rates.

As SARS-CoV-2 continues to spread worldwide, tractable primary airway cell models that accurately recapitulate the cell-intrinsic response to arising viral variants are needed. Here we describe an adult stem cell-derived human airway organoid model overexpressing the ACE2 receptor that supports robust viral replication while maintaining 3D architecture and cellular diversity of the airway epithelium. ACE2-OE organoids were infected with SARS-CoV-2 variants and subjected to single-cell RNA-sequencing. NF-κB inhibitor alpha was consistently upregulated in infected epithelial cells, and its mRNA expression positively correlated with infection levels. Confocal microscopy showed more IκBα expression in infected than bystander cells, but found concurrent nuclear translocation of NF-κB that IκBα usually prevents. Overexpressing a nondegradable IκBα mutant reduced NF-κB translocation and increased viral infection. These data demonstrate the functionality of ACE2-OE organoids in SARS-CoV-2 research and identify an incomplete NF-κB feedback loop as a rheostat of viral infection that may promote inflammation and severe disease.

BackgroundThe use of graphs as a way of abstracting and representing biological systems has provided a powerful analysis paradigm. Specifically, graph optimization algorithms are routinely used to address various connectivity queries, such as finding paths between proteins in a protein-protein interaction network, while maximizing objectives such as parsimony. While present studies in this field mostly concern static graphs, new types of data now motivate the need to account for changes that might occur to the elements (nodes) that are represented by the graph on the relationships (edges) between them.Results and DiscussionWe define the notion of Directed Temporal Graphs as a series of directed subgraphs of an underlying graph, ordered by time, where only a subset of vertices and edges are present. We then build up towards the Time Conditioned Shortest Walk problem on Directed Temporal Graphs: given a series of time ordered directed graphs, find the shortest walk from any given source node at time point 1 to a target node at time T ≥ 1, such that the walk is consistent (monotonically increasing) with the timing of nodes and edges. We show, contrary to the Directed Shortest Walk problem which can be solved in polynomial time, that the Time Conditioned Shortest Walk (TCSW) problem is NP-Hard, and is hard to approximate to factor for T ≥ 3 and ε > 0. Lastly, we develop an integer linear program to solve a generalized version of TCSW, and demonstrate its ability to reach optimality with instances of the human protein interaction network.ConclusionWe demonstrate that when extending the shortest walk problem in computational biology to account for multiple ordered conditions, the problem not only becomes hard to solve, but hard to approximate, a limitation which we address via a new solver. From this narrow definition of TCSW, we relax the constraint of time consistency within the shortest walk, deriving a direct relationship between hardness of approximation and the allowable step size in our walk between time conditioned networks. Lastly we briefly explore a variety of alternative formulations for this problem, providing insight into both tractable and intractable variants.AvailabilityOur solver for the general k-Time Condition Shortest Walk problem is available at https://github.com/YosefLab/temporal_condition_shortest_walk

We developed MEBOCOST, a computational algorithm for quantitatively inferring metabolite-mediated intercellular communications using single cell RNA-seq data. The algorithm identifies cell-cell communications in which metabolites, such as lipids, are secreted by sender cells and traveled to interact with sensor proteins of receiver cells. The sensor proteins on receiver cell might be cell surface receptors, transporters across the cell membrane, or nuclear receptors. MEBOCOST relies on a comprehensive database of metabolite-sensor partners, which we manually curated from the literatures and other public sources. MEBOCOST defines sender and receiver cells for an extracellular metabolite based on the expression levels of the enzymes and sensors, respectively, thus identifies metabolite-sensor communications between the cells. Applying MEBOCOST to mouse brown adipose tissue (BAT) successfully recaptured known metabolite-mediated cell communications and further identified new communications. Additionally, MEBOCOST identified a set of BAT intercellular metabolite-sensor communications that was regulated by cold exposure of the mice. MEBOCOST will be useful to numerous researchers to investigate metabolite-mediated cell-cell communications in many biological and disease models. The MEBOCOST software is freely available at https://github.com/zhengrongbin/MEBOCOST.

A major limitation in single-cell genomics is a lack of ability to conduct cost-effective population-level studies. As a result, much of the current research in single-cell genomics focuses on biological processes that are broadly conserved across individuals, such as cellular organization and tissue development. This limitation prevents us from studying the etiology of experimental or clinical conditions that may be inconsistent across individuals owing to molecular variation and a wide range of effects in the population. In order to address this gap, we developed “kernel of integrated single cells” (Keris), a novel model-based framework to inform the analysis of single-cell gene expression data with population-level effects of a condition of interest. By inferring cell-type-specific moments and their variation across conditions using large tissue-level bulk data representing a population, Keris allows us to generate testable hypotheses at the single-cell level that would otherwise require collecting single-cell data from a large number of donors. Within the Keris framework, we show how the combination of low-resolution, large bulk data with small but high-resolution single-cell data enables the identification of changes in cell-subtype compositions and the characterization of subpopulations of cells that are affected by a condition of interest. Using Keris we estimate linear and non-linear age-associated changes in cell-type expression in large bulk peripheral blood mononuclear cells (PBMC) data. Combining with three independent single-cell PBMC datasets, we demonstrate that Keris can identify changes in cell-subtype composition with age and capture cell-type-specific subpopulations of senescent cells. This demonstrates the promise of enhancing single-cell data with population-level information to study compositional changes and to profile condition-affected subpopulations of cells, and provides a potential resource of targets for future clinical interventions.

Lineage-tracing technologies based on Clustered Regularly Interspaced Short Palindromic Repeats and CRISPR-associated protein 9 (CRISPR-Cas9) genome editing have emerged as a powerful tool for investigating development in single-cell contexts, but exact reconstruction of the underlying clonal relationships in experiment is complicated by features of the data. These complications are functions of the experimental parameters in these systems, such as the Cas9 cutting rate, the diversity of indel outcomes, and the rate of missing data. In this paper, we develop two theoretically grounded algorithms for the reconstruction of the underlying single-cell phylogenetic tree as well as asymptotic bounds for the number of recording sites necessary for exact recapitulation of the ground truth phylogeny at high probability. In doing so, we explore the relationship between the problem difficulty and the experimental parameters, with implications for experimental design. Lastly, we provide simulations showing the empirical performance of these algorithms and showing that the trends in the asymptotic bounds hold empirically. Overall, this work provides a theoretical analysis of phylogenetic reconstruction in single-cell CRISPR-Cas9 lineage-tracing technologies.

SummaryHepatitis C virus (HCV) is the leading cause of death from liver disease. How HCV infection causes lasting liver damage and increases cancer risk beyond viral clearance remains unclear. We identify bipotent liver stem cells as novel targets for HCV infection, and their erroneous differentiation as the potential cause of impaired liver regeneration and cancer development. We show 3D organoids generated from liver stem cells from actively HCV-infected individuals carry replicating virus and maintain low-grade infection over months. Organoids can be infected with a primary HCV isolate. Virus-inclusive single-cell RNA-sequencing uncovered extensive transcriptional reprogramming in HCV+ cells supporting hepatocytic differentiation, cancer stem cell development and viral replication while stem cell proliferation and interferon signaling are disrupted. Our data adds a pathogenesis factor – infection of liver stem cells – to the biology of HCV infection that explains persistent liver damage and enhanced cancer risk through an altered stem cell state.

SUMMARYTumor evolution is driven by the progressive acquisition of genetic and epigenetic alterations that enable uncontrolled growth, expansion to neighboring and distal tissues, and therapeutic resistance. The study of phylogenetic relationships between cancer cells provides key insights into these processes. Here, we introduced an evolving lineage-tracing system with a single-cell RNA-seq readout into a mouse model of Kras;Trp53(KP)-driven lung adenocarcinoma which enabled us to track tumor evolution from single transformed cells to metastatic tumors at unprecedented resolution. We found that loss of the initial, stable alveolar-type2-like state was accompanied by transient increase in plasticity. This was followed by adoption of distinct fitness-associated transcriptional programs which enable rapid expansion and ultimately clonal sweep of rare, stable subclones capable of metastasizing to distant sites. Finally, we showed that tumors develop through stereotypical evolutionary trajectories, and perturbing additional tumor suppressors accelerates tumor progression by creating novel evolutionary paths. Overall, our study elucidates the hierarchical nature of tumor evolution, and more broadly enables the in-depth study of tumor progression.

SUMMARYRecent advances in CRISPR-Cas9 engineering and single-cell assays have enabled the simultaneous measurement of single-cell transcriptomic and phylogenetic profiles. However, there are few computational tools enabling users to integrate and derive insight from a joint analysis of these two modalities. Here, we describe PhyloVision: an open source software for interactively exploring data from both modalities and for identifying and interpreting heritable gene modules whose concerted expression are associated with phylogenetic relationships. PhyloVision provides a feature-rich, interactive, and shareable web-based report for investigating these modules, while also supporting several other data and meta-data exploration capabilities. We demonstrate the utility of PhyloVision using a published dataset of metastatic lung adenocarcinoma cells, whose phylogeny was resolved using a CRISPR/Cas9-based lineage-tracing system. Together, we anticipate that PhyloVision and the methods it implements will be a useful resource for scalable and intuitive data exploration for any assay that simultaneously measures cell state and lineage.

The ability to jointly profile the transcriptional and chromatin land-scape of single-cells has emerged as a powerful technique to identify cellular populations and shed light on their regulation of gene expression. Current computational methods analyze jointly profiled (paired) or individual data modalities (unpaired), but do not offer a principled method to analyze both paired and unpaired samples jointly. Here we present MultiVI, a probabilistic framework that leverages deep neural networks to jointly analyze scRNA, scATAC and multiomic (scRNA + scATAC) data. MultiVI creates an informative low-dimensional latent space that accurately reflects both chromatin and transcriptional properties of the cells even when one of the modalities is missing. MultiVI accounts for technical effects in both scRNA and scATAC-seq while correcting for batch effects in both data modalities. We use public datasets to demonstrate that MultiVI is stable, easy to use, and outperforms current approaches for the joint analysis of paired and unpaired data. MultiVI is available as an open source package, implemented in the scvi-tools frame-work: https://docs.scvi-tools.org/.

Alternative splicing shapes the transcriptome and contributes to each cell’s unique identity, but single-cell RNA sequencing has struggled to capture the impact of alternative splicing. We previously showed that low recovery of mRNAs from single cells led to erroneous conclusions about the cell-to-cell variability of alternative splicing (1). Here, we present a method, Psix, to confidently identify splicing that changes across a landscape of single cells, using a probabilistic model that is robust against the data limitations of scRNA-seq. Its autocorrelation-inspired approach finds patterns of alternative splicing that correspond to patterns of cell identity, such as cell type or developmental stage, without the need for explicit cell clustering, labeling, or trajectory inference. Applying Psix to data that follow the trajectory of mouse brain development, we identify exons whose alternative splicing patterns cluster into modules of co-regulation. We show that the exons in these modules are enriched for binding by distinct neuronal splicing factors, and that their changes in splicing correspond to changes in expression of these splicing factors. Thus, Psix reveals cell-type-dependent splicing patterns and the wiring of the splicing regulatory networks that control them. Our new method will enable scRNA-seq analysis to go beyond transcription to understand the roles of post-transcriptional regulation in determining cell identity.

The development of CD4 and CD8 T cells in the thymus is critical to adaptive immunity and is widely studied as a model of lineage commitment. Recognition of self-MHCI/II by the T cell antigen receptor (TCR) determines the lineage choice, but how distinct TCR signals drive transcriptional programs of lineage commitment remains largely unknown. We applied CITE-seq to measure RNA and surface proteins in thymocytes from wild-type and lineage-restricted mice to generate a comprehensive timeline of cell state for each lineage. These analyses revealed a sequential process whereby all thymocytes initiate CD4 lineage differentiation during an initial wave of TCR signaling, followed by a second TCR signaling wave that coincides with CD8 lineage specification. CITE-seq and pharmaceutical inhibition experiments implicate a TCR/calcineurin/NFAT/GATA3 axis in driving the CD4 fate. Overall, our data suggest that multiple redundant mechanisms contribute to the accuracy and efficiency of the lineage choice.

Whether regeneration is primarily accomplished by re-activating gene regulatory networks used previously during development or by activating novel regeneration-specific transcriptional programs remains a longstanding question. Currently, most genes implicated in regeneration also function during development. Using single-cell transcriptomics in regenerating Drosophila wing discs, we identified two regeneration-specific cell populations within the blastema. They are each composed of cells that upregulate multiple genes encoding secreted proteins that promote regeneration. In this regenerative secretory zone, the transcription factor Ets21C controls the expression of multiple regenerationpromoting genes. While eliminating Ets21C function has no discernible effect on development, it severely compromises regeneration. This Ets21C-dependent gene regulatory network is also activated in blastema-like cells in tumorous discs, suggesting that pro-regenerative mechanisms can be co-opted by tumors to promote aberrant growth.

Novel experimental assays now simultaneously measure lineage relationships and transcriptomic states from single cells, thanks to CRISPR/Cas9-based genome engineering. These multimodal measurements allow researchers not only to build comprehensive phylogenetic models relating all cells but also infer transcriptomic determinants of consequential subclonal behavior. The gene expression data, however, is limited to cells that are currently present (“leaves” of the phylogeny). As a consequence, researchers cannot form hypotheses about unobserved, or “ancestral”, states that gave rise to the observed population. To address this, we introduce TreeVAE: a probabilistic framework for estimating ancestral transcriptional states. TreeVAE uses a variational autoencoder (VAE) to model the observed transcriptomic data while accounting for the phylogenetic relationships between cells. Using simulations, we demonstrate that TreeVAE outperforms benchmarks in reconstructing ancestral states on several metrics. TreeVAE also provides a measure of uncertainty, which we demonstrate to correlate well with its prediction accuracy. This estimate therefore potentially provides a data-driven way to estimate how far back in the ancestor chain predictions could be made. Finally, using real data from lung cancer metastasis, we show that accounting for phylogenetic relationship between cells improves goodness of fit. Together, TreeVAE provides a principled framework for reconstructing unobserved cellular states from single cell lineage tracing data.

Emerging evidence indicates a fundamental role for the epigenome in immunity. Here, we used a systems biology approach to map the epigenomic and transcriptional landscape of immunity to influenza vaccination in humans at the single-cell level. Vaccination against seasonal influenza resulted in persistently reduced H3K27ac in monocytes and myeloid dendritic cells, which was associated with impaired cytokine responses to TLR stimulation. Single cell ATAC-seq analysis of 120,305 single cells revealed an epigenomically distinct subcluster of monocytes with reduced chromatin accessibility at AP-1-targeted loci after vaccination. Similar effects were also observed in response to vaccination with the AS03-adjuvanted H5N1 pandemic influenza vaccine. However, this vaccine also stimulated persistently increased chromatin accessibility at loci targeted by interferon response factors (IRFs). This was associated with elevated expression of antiviral genes and type 1 IFN production and heightened resistance to infection with the heterologous viruses Zika and Dengue. These results demonstrate that influenza vaccines stimulate persistent epigenomic remodeling of the innate immune system. Notably, AS03-adjuvanted vaccination remodeled the epigenome of myeloid cells to confer heightened resistance against heterologous viruses, revealing its potentially unappreciated role as an epigenetic adjuvant.

The function of mammalian cells is largely influenced by their tissue microenvironment. Advances in spatial transcriptomics open the way for studying these important determinants of cellular function by enabling a transcriptome-wide evaluation of gene expression in situ. A critical limitation of the current technologies, however, is that their resolution is limited to niches (spots) of sizes well beyond that of a single cell, thus providing measurements for cell aggregates which may mask critical interactions between neighboring cells of different types. While joint analysis with single-cell RNA-sequencing (scRNA-seq) can be leveraged to alleviate this problem, current analyses are limited to a discrete view of cell type proportion inside every spot. This limitation becomes critical in the common case where, even within a cell type, there is a continuum of cell states that cannot be clearly demarcated but reflects important differences in the way cells function and interact with their surroundings. To address this, we developed Deconvolution of Spatial Transcriptomics profiles using Variational Inference (DestVI), a probabilistic method for multi-resolution analysis for spatial transcriptomics that explicitly models continuous variation within cell types. Using simulations, we demonstrate that DestVI is capable of providing higher resolution compared to the existing methods and that it can estimate gene expression by every cell type inside every spot. We then introduce an automated pipeline that uses DestVI for analysis of single tissue slices and comparison between tissues. We apply this pipeline to study the immune crosstalk within lymph nodes to infection and explore the spatial organization of a mouse tumor model. In both cases, we demonstrate that DestVI can provide a high resolution and accurate spatial characterization of the cellular organization of these tissues, and that it is capable of identifying important cell-type-specific changes in gene expression - between different tissue regions or between conditions. DestVI is available as an open-source software package in the scvi-tools codebase (https://scvi-tools.org).

Single-cell ATAC sequencing (scATAC-seq) is a powerful and increasingly popular technique to explore the regulatory landscape of heterogeneous cellular populations. However, the high noise levels, degree of sparsity, and scale of the generated data make its analysis challenging. Here we present PeakVI, a probabilistic framework that leverages deep neural networks to analyze scATAC-seq data. PeakVI fits an informative latent space that preserves biological heterogeneity while correcting batch effects and accounting for technical effects such as library size and region-specific biases. Additionally, PeakVI provides a technique for identifying differential accessibility at a single region resolution, which can be used for cell-type annotation as well as identification of key cis-regulatory elements. We use public datasets to demonstrate that PeakVI is scalable, stable, robust to low-quality data, and outperforms current analysis methods on a range of critical analysis tasks. PeakVI is publicly available and implemented in the scvi-tools framework: https://docs.scvi-tools.org/.

Gene regulatory elements play a key role in orchestrating gene expression during cellular differentiation, but what determines their function over time remains largely unknown. Here, we performed perturbation-based massively parallel reporter assays at seven early time points of neural differentiation to systematically characterize how regulatory elements and motifs within them guide cellular differentiation. By perturbing over 2,000 putative DNA binding motifs in active regulatory regions, we delineated four categories of functional elements, and observed that activity direction is mostly determined by the sequence itself, while the magnitude of effect depends on the cellular environment. We also find that fine-tuning transcription rates is often achieved by a combined activity of adjacent activating and repressing elements. Our work provides a blueprint for the sequence components needed to induce different transcriptional patterns in general and specifically during neural differentiation.

SummaryThe cytokine receptor IL-23R plays a fundamental role in inflammation and autoimmunity. However, several observations have been difficult to reconcile under the assumption that only Th17 cells critically depend on IL-23 to acquire a pathogenic phenotype. Here, we report that Th1 cells differentiated in vitro with IL-12 + IL-21 show similar levels of IL-23R expression as in pathogenic Th17 cells. We demonstrate that IL-23R is required for Th1 cells to acquire a highly colitogenic phenotype. scRNAseq analysis of intestinal T cells enabled us to identify novel regulators induced by IL-23R-signaling in Th1 cells which differed from those expressed in Th17 cells. The perturbation of one of these regulators (CD160) in Th1 cells inhibited induction of colitis. In this process, we were able to uncouple IL-23R as a purely Th17 cell-specific factor and implicate IL-23R signaling as a pathogenic driver of Th1 cell-mediated tissue inflammation and disease.

SUMMARYThe response of cytotoxic T cells to their cognate antigen involves rapid and broad changes in gene expression that are coupled with extensive chromatin remodeling. Here, we study the mechanisms by which the basic leucine zipper ATF-like transcription factor Batf helps regulate this process. Through genome-scale profiling, we observe critical roles for Batf in inducing transcriptional changes in stimulated naive cells, while affecting the chromatin at several levels, namely binding of key transcription factors, accessibility, and long range contacts. We identify a critical network of transcription factors that cooperate with Batf, including its binding partner Irf4, as well as Runx3 and T-bet, and demonstrate its synergistic activity in initiating aspects of the effector T cells’ transcriptional and epigenetic program in ectopically-induced fibroblasts. Our results provide a comprehensive resource for studying the epigenomic and transcriptomic landscape of effector differentiation of cytotoxic T cells.

The only licensed vaccine for tuberculosis, Mycobacterium bovis Bacille Calmette-Guérin (BCG), is not reliably effective against adult pulmonary tuberculosis. A major hurdle to tuberculosis vaccine development is incomplete understanding of successful immunity against the causative agent Mycobacterium tuberculosis. Recently, we demonstrated that a protein subunit vaccine adjuvanted with STING-activating cyclic-di-nucleotides (CDNs) robustly protects against tuberculosis infection in mice. Here we show mucosal vaccination with this vaccine induces production of T cells that home to lung parenchyma and penetrate lesions in the lung. Protection is partially dependent on IL-17, type I interferon (IFN), and IFN-γ, while the transcription factor STAT-6 is dispensable. Single cell transcriptomics reveals mucosal vaccination with a CDN vaccine increases transcriptional heterogeneity in CD4 cells, including a significant population of non-classical IFN-γ and IL-17 co-expressing Th1-Th17 cells, as well as markers of memory and activation. Th1-Th17 cells in vaccinated mice are enriched for expression of the T cell functional markers Tnfsf8 and Il1r1 relative to more conventional Th1 cells. These data provide critical insight into the immune mediators and diversity of T cell responses that can contribute to vaccine efficacy against M. tuberculosis infection.

In both vertebrates and invertebrates, generating a functional appendage requires interactions between ectoderm-derived epithelia and mesoderm-derived cells. To investigate such interactions, we used single-cell transcriptomics to generate a cell atlas of theDrosophilawing disc at two time points during development. Using these data, we investigate gene expression using a multi-layered model of the wing disc and catalogued ligand-receptor pairs that could mediate signaling between epithelial cells and adult muscle precursors (AMPs). We found that localized expression of the FGF ligands, Thisbe and Pyramus, in the disc epithelium regulates the number and location of the AMPs. In addition, Hedgehog ligand from the epithelium activates a specific transcriptional program within adjacent AMP cells, which is critical for proper formation of a subset of the direct flight muscles. More generally, our annotated atlas provides a global view of potential cell-cell interactions between subpopulations of epithelial and myogenic cells.

The paired measurement of RNA and surface protein abundance in single cells with CITE-seq is a promising approach to connect transcriptional variation with cell phenotypes and functions. However, each data modality exhibits unique technical biases, making it challenging to conduct a joint analysis and combine these two views into a unified representation of cell state. Here we present Total Variational Inference (totalVI), a framework for the joint probabilistic analysis of paired RNA and protein data from single cells. totalVI probabilistically represents the data as a composite of biological and technical factors such as limited sensitivity of the RNA data, background in the protein data, and batch effects. To evaluate totalVI, we performed CITE-seq on immune cells from murine spleen and lymph nodes with biological replicates and with different antibody panels measuring over 100 surface proteins. With this dataset, we demonstrate that totalVI provides a cohesive solution for common analysis tasks like the integration of datasets with matched or unmatched protein panels, dimensionality reduction, clustering, evaluation of correlations between molecules, and differential expression testing. totalVI enables scalable, end-to-end analysis of paired RNA and protein data from single cells and is available as open-source software.

Cancer progression is characterized by rare, transient events which are nonetheless highly consequential to disease etiology and mortality. Detailed cell phylogenies can recount the history and chronology of these critical events – including metastatic seeding. Here, we applied our Cas9-based lineage tracer to study the subclonal dynamics of metastasis in a lung cancer xenograft mouse model, revealing the underlying rates, routes, and drivers of metastasis. We report deeply resolved phylogenies for tens of thousands of metastatically disseminated cancer cells. We observe surprisingly diverse metastatic phenotypes, ranging from metastasis-incompetent to aggressive populations. These phenotypic distinctions result from pre-existing, heritable, and characteristic differences in gene expression, and we demonstrate that these differentially expressed genes can drive invasiveness. Furthermore, metastases transit via diverse, multidirectional tissue routes and seeding topologies. Our work demonstrates the power of tracing cancer progression at unprecedented resolution and scale.One Sentence SummarySingle-cell lineage tracing and RNA-seq capture diverse metastatic behaviors and drivers in lung cancer xenografts in mice.

FoxP3 positive regulatory T cells (Tregs) rely on fatty acid β-oxidation (FAO)-driven OXPHOS for differentiation and function. Recent data have demonstrated a role for Tregs in the maintenance of tissue homeostasis with tissue-resident Tregs possessing tissue-specific transcriptomes. However, specific signals that establish these tissue-resident Tregs programs are largely unknown. As Tregs metabolically rely on FAO, and considering the lipid-rich environments of tissues, we hypothesized that environmental lipids drive Treg homeostasis. Using human adipose tissue as a model for tissue residency, we identify oleic acid as the most prevalent free fatty acid in human adipose tissue. Mechanistically, oleic acid amplifies Treg FAO-driven OXPHOS metabolism, creating a positive feedback mechanism that induces the expression of Foxp3 and enhances phosphorylation of STAT5, which acts to stabilize the Treg lineage and increase suppressive function. Comparing the transcriptomic program induced by oleic acid to that of the pro-inflammatory arachidonic acid, we find that Tregs sorted from peripheral blood and adipose of healthy donors transcriptomically resemble the oleic acid in vitro treated Tregs, whereas Tregs obtained from the adipose tissue of relapsing-remitting MS patients more closely resemble an arachidonic acid profile. Finally, we find that oleic acid concentrations are reduced in the fat tissue of MS patients, and exposure of dysfunctional MS Tregs to oleic acid restores defects in their suppressive function. These data demonstrate the importance of fatty acids in regulating tissue inflammatory signals.

SUMMARYWe report here a central role for polyamines in T cell differentiation and function. Deficiency in ornithine decarboxylase (ODC), a critical enzyme for polyamine synthesis, resulted in a profound failure of CD4+ T cells to adopt correct subset specification, underscored by ectopic expression of multiple cytokines and lineage-defining transcription factors across TH1, TH2, TH17, and Treg polarizing conditions, and enhanced colitogenic potential. T cells deficient in deoxyhypusine synthase (DHPS) or deoxyhypusine hydroxylase (DOHH), which sequentially utilize polyamines to generate hypusine, phenocopied Odc-deficient T cells, and mice in which T cells lacked Dhps or Dohh developed colitis. Polyamine-hypusine pathway enzyme deficiency caused widespread chromatin and transcriptional dysregulation accompanied by alterations in histone methylation, histone acetylation, and TCA cycle metabolites. Epigenetic modulation by 2-hydroxyglutarate, or histone acetyltransferase inhibition, restored CD4+ T cell subset specification. Thus, polyamine synthesis via hypusine is critical for maintaining the epigenome to focus TH cell subset fidelity.

Cellular metabolism, a key regulator of immune responses, is difficult to study with current technologies in individual cells Here, we present Compass, an algorithm to characterize the metabolic state of cells based on single-cell RNA-Seq and flux balance analysis. We applied Compass to associate metabolic states with functional variability (pathogenic potential) amongst Th17 cells and recovered a metabolic switch between glycolysis and fatty acid oxidation, akin to known differences between Th17 and Treg cells, as well as novel targets in amino-acid pathways, which we tested through targeted metabolic assays. Compass further predicted a particular glycolytic reaction (phosphoglycerate mutase — PGAM) that promotes an anti-inflammatory Th17 phenotype, contrary to the common understanding of glycolysis as pro-inflammatory. We demonstrate that PGAM inhibition leads non-pathogenic Th17 cells to adopt a pro-inflammatory transcriptome and induce autoimmunity in vivo. Compass is broadly applicable for characterizing metabolic states of cells and relating metabolic heterogeneity to other cellular phenotypes.

Cellular metabolism can orchestrate immune cell function. We previously demonstrated that lipid biosynthesis represents one such gatekeeper to Th17 cell functional state. Utilizing Compass, a transcriptome-based algorithm for prediction of metabolic flux, we constructed a comprehensive metabolic circuitry for Th17 cell function and identified the polyamine pathway as a candidate metabolic node, the flux of which regulates the inflammatory function of T cells. Testing this prediction, we found that expression and activities of enzymes of the polyamine pathway were enhanced in pathogenic Th17 cells and suppressed in regulatory T cells. Perturbation of the polyamine pathway in Th17 cells suppressed canonical Th17 cell cytokines and promoted the expression of Foxp3, accompanied by dramatic shift in transcriptome and epigenome, transitioning Th17 cells into a Treg-like state. Genetic and chemical perturbation of the polyamine pathway resulted in attenuation of tissue inflammation in an autoimmune disease model of central nervous system, with changes in T cell effector phenotype.

The pairing of CRISPR/Cas9-based gene editing with massively parallel single-cell readouts now enables large-scale lineage tracing. However, the rapid growth in complexity of data from these assays has outpaced our ability to accurately infer phylogenetic relationships. To address this, we provide three resources. First, we introduce Cassiopeia - a suite of scalable and theoretically grounded maximum parsimony approaches for tree reconstruction. Second, we provide a simulation framework for evaluating algorithms and exploring lineage tracer design principles. Finally, we generate the most complex experimental lineage tracing dataset to date - consisting of 34,557 human cells continuously traced over 15 generations, 71% of which are uniquely marked - and use it for benchmarking phylogenetic inference approaches. We show that Cassiopeia outperforms traditional methods by several metrics and under a wide variety of parameter regimes, and provide insight into the principles for the design of improved Cas9-enabled recorders. Together these should broadly enable large-scale mammalian lineage tracing efforts. Cassiopeia and its benchmarking resources are publicly available at www.github.com/YosefLab/Cassiopeia.

In single-cell RNA sequencing data, biological processes or technical factors may induce an overabundance of zero measurements. Existing probabilistic approaches to interpreting these data either model all genes as zero-inflated, or none. But the overabundance of zeros might be gene-specific. Hence, we propose the AutoZI model, which, for each gene, places a spike-and-slab prior on a mixture assignment between a negative binomial (NB) component and a zero-inflated negative binomial (ZINB) component. We approximate the posterior distribution under this model using variational inference, and employ Bayesian decision theory to decide whether each gene is zero-inflated. On simulated data, AutoZI outperforms the alternatives. On negative control data, AutoZI retrieves predictions consistent to a previous study on ERCC spike-ins and recovers similar results on control RNAs. Applied to several datasets and instances of the 10x Chromium protocol, AutoZI allows both biological and technical interpretations of zero-inflation. Finally, AutoZI’s decisions on mouse embyronic stem-cells suggest that zero-inflation might be due to transcriptional bursting.

Cellular indexing of transcriptomes and epitopes by sequencing (CITE-seq) combines unbiased single-cell transcriptome measurements with surface protein quantification comparable to flow cytometry, the gold standard for cell type identification. However, current analysis pipelines cannot address the two primary challenges of CITE-seq data: combining both modalities in a shared latent space that harnesses the power of the paired measurements, and handling the technical artifacts of the protein measurement, which is obscured by non-negligible background noise. Here we present Total Variational Inference (totalVI), a fully probabilistic end-to-end framework for normalizing and analyzing CITE-seq data, based on a hierarchical Bayesian model. In totalVI, the mRNA and protein measurements for each cell are generated from a low-dimensional latent random variable unique to that cell, representing its cellular state. totalVI uses deep neural networks to specify conditional distributions. By leveraging advances in stochastic variational inference, it scales easily to millions of cells. Explicit modeling of nuisance factors enables totalVI to produce denoised data in both domains, as well as a batch-corrected latent representation of cells for downstream analysis tasks.

mRNA levels are determined by the balance between mRNA synthesis and decay. Factors that mediate both processes, including the 5’ to 3’ exonuclease Xrn1, are responsible for the cross talk between the two processes in a manner that buffers steady-state mRNA levels. However, these proteins’ roles in transcription remain elusive and controversial. Applying NET-seq to yeast cells, we show that Xrn1 functions mainly as a transcriptional activator and that its disruption manifests via the reduction of RNA polymerase II (Pol II) occupancy downstream of transcription start sites. We combine our data and novel mathematical modeling of transcription to suggest that transcription initiation and elongation of targeted genes is modulated by Xrn1. Furthermore, Pol II occupancy markedly increases near cleavage and polyadenylation sites in xrn1Δ cells while its activity decreases, a characteristic feature of backtracked Pol II. We also provide indirect evidence that Xrn1 is involved in transcription termination downstream of polyadenylation sites. Two additional decay factors, Dhh1 and Lsm1, seem to function similarly to Xrn1 in transcription, perhaps as a complex, while the decay factors Ccr4 and Rpb4 also perturb transcription in other ways. Interestingly, DFs are capable of differentiating between SAGA- and TFIID-dominated promoters. These two classes of genes respond differently to XRN1 deletion in mRNA synthesis and differentially utilize mRNA decay pathways, raising the possibility that one distinction between the two types of genes lies in the mechanism(s) that balance these processes.

Massively parallel reporter assays (MPRAs) are a technique that enables testing thousands of regulatory DNA sequences and their variants in a single, quantitative experiment. Despite growing popularity, there is lack of statistical methods that account for the different sources of uncertainty inherent to these assays, thus effectively leveraging their promise. Development of such methods could help enhance our ability to identify regulatory sequences in the genome, understand their function under various setting, and ultimately gain a better understanding of how the regulatory code and its alteration lead to phenotypic consequence.Here we present MPRAnalyze: a statistical framework dedicated to analyzing MPRA count data. MPRAnalyze addresses the major questions that are posed in the context of MPRA experiments: estimating the magnitude of the effect of a regulatory sequence in a single condition setting, and comparing differential activity of regulatory sequences across multiple conditions. The framework uses a nested construction of generalized linear models to account for uncertainty in both DNA and RNA observations, controls for various sources of unwanted variation, and incorporates negative controls for robust hypothesis testing, thereby providing clear quantitative answers in complex experimental settings.We demonstrate the robustness, accuracy and applicability of MPR-Analyze on simulated data and published data sets and compare it against the existing analysis methodologies. MPRAnalyze is implemented as an R package and is publicly available through Bioconductor [1].

We present VISION, a tool for annotating the sources of variation in single cell RNA-seq data in an automated, unbiased and scalable manner. VISION operates directly on the manifold of cell-cell similarity and employs a flexible annotation approach that can operate either with or without preconceived stratification of the cells into groups or along a continuum. We demonstrate the utility of VISION using a relatively homogeneous set of B cells from a cohort of lupus patients and healthy controls and show that it can derive important sources of cellular variation and link them to clinical phenotypes in a stratification free manner. VISION produces an interactive, low latency and feature rich web-based report that can be easily shared amongst researchers.

Cerebrospinal fluid (CSF) protects the central nervous system (CNS) and analyzing CSF aids the diagnosis of CNS diseases, but our understanding of CSF leukocytes remains superficial. Here, using single cell transcriptomics, we identify a specific location-associated composition and transcriptome of CSF leukocytes. Multiple sclerosis (MS) – an autoimmune disease of the CNS – increases transcriptional diversity in blood, but increases cell type diversity in CSF including a higher abundance of cytotoxic phenotype T helper cells. An analytical approach, named cell set enrichment analysis (CSEA) identifies a cluster-independent increase of follicular (TFH) cells potentially driving the known expansion of B lineage cells in the CSF in MS. In mice, TFH cells accordingly promote B cell infiltration into the CNS and the severity of MS animal models. Immune mechanisms in MS are thus highly compartmentalized and indicate ongoing local T/B cell interaction., Here the authors provide a single-cell characterization of cerebrospinal fluid and blood of newly diagnosed multiple sclerosis (MS) patients, revealing altered composition of lymphocyte and monocyte subsets, validated by other methods including the interrogation of the TFH subset in mouse models of MS.

RNA-sequencing of single B cells provides simultaneous measurements of the cell state and its binding specificity. However, in order to uncover the latter further reconstruction of the B cell receptor (BCR) sequence is needed. We present BRAPeS, an algorithm for reconstructing BCRs from short-read paired-end single cell RNA-sequencing. BRAPeS is accurate and achieves a high success rate even at very short (25bp) read length, which can decrease the cost and increase the number of cells that can be analyzed compared to long reads. BRAPeS is publicly available in the following link: https://github.com/YosefLab/BRAPeS.

The abundance of new computational methods for processing and interpreting transcriptomes at a single cell level raises the need for in-silico platforms for evaluation and validation. Simulated datasets which resemble the properties of real datasets can aid in method development and prioritization as well as in questions in experimental design by providing an objective ground truth. Here, we present SymSim, a simulator software that explicitly models the processes that give rise to data observed in single cell RNA-Seq experiments. The components of the SymSim pipeline pertain to the three primary sources of variation in single cell RNA-Seq data: noise intrinsic to the process of transcription, extrinsic variation that is indicative of different cell states (both discrete and continuous), and technical variation due to low sensitivity and measurement noise and bias. Unlike other simulators, the parameters that govern the simulation process directly represent meaningful properties such as mRNA capture rate, the number of PCR cycles, sequencing depth, or the use of unique molecular identifiers. We demonstrate how SymSim can be used for benchmarking methods for clustering and differential expression and for examining the effects of various parameters on their performance. We also show how SymSim can be used to evaluate the number of cells required to detect a rare population and how this number deviates from the theoretical lower bound as the quality of the data decreases. SymSim is publicly available as an R package and allows users to simulate datasets with desired properties or matched with experimental data.

The molecular components governing neural induction remain largely unknown. Here, we applied a suite of genomic and computational tools to comprehensively identify these components. We performed RNA-seq, ChIP-seq (H3K27ac, H3K27me3) and ATAC-seq on human embryonic stem cells (hESCs) at seven early neural differentiation time points (0-72 hours) and identified thousands of induced genes and regulatory regions. We analyzed the function of ~2,500 selected regions using massively parallel reporter assays at all time points. We found numerous temporal enhancers that correlated with similarly timed epigenetic marks and gene expression. Development of a prioritization method that incorporated all genomic data identified key transcription factors (TFs) involved in neural induction. Individual overexpression of eleven TFs and several combinations in hESCs found novel neural induction regulators. Combined, our results provide a comprehensive map of genes and functional regulatory elements involved in neural induction and identify master regulator TFs that are instrumental for this process.One Sentence SummaryUsing numerous genomic assays and computational tools we characterized the dynamic changes that take place during neural induction.

SummaryThe decreasing cost of DNA sequencing over the past decade has led to an explosion of available sequencing datasets, leaving us with terabytes to petabytes of data to explore and analyze. It is critical for analysts in research and clinical settings to be able to develop new data-driven hypotheses from these datasets through bias identification, analysis of data quality, and testing different algorithms and parameter settings. However, current interactive tools for sequence analysis are designed to run on single machines that do not scale to the size of modern genomic datasets, and rely on precomputed static views, rather than allowing direct interaction with the primary dataset. Mango is a genomic sequence visualization and analysis platform that removes these constraints regarding scalability and staticity by leveraging the power of multi-node compute clusters in the cloud to allow interactive analysis over terabytes of sequencing data. Mango provides both a genome browser graphical user interface and programmable notebook form factor to allow users of varying analytical experience to explore large sequencing datasets on both private clusters and in the cloud. These tools provide a flexible environment for interactive exploration of genomic datasets, while surpassing the computational limits of single-node genomic visualization tools.

Transcriptome profiles of individual cells reflect true and often unexplored biological diversity, but are also affected by noise of biological and technical nature. This raises the need to explicitly model the resulting uncertainty and take it into account in any downstream analysis, such as dimensionality reduction, clustering, and differential expression. Here, we introduce Single-cell Variational Inference (scVI), a scalable framework for probabilistic representation and analysis of gene expression in single cells. Our model uses variational inference and stochastic optimization of deep neural networks to approximate the parameters that govern the distribution of expression values of each gene in every cell, using a non-linear mapping between the observations and a low-dimensional latent space.By doing so, scVI pools information between similar cells or genes while taking nuisance factors of variation such as batch effects and limited sensitivity into account. To evaluate scVI, we conducted a comprehensive comparative analysis to existing methods for distributional modeling and dimensionality reduction, all of which rely on generalized linear models. We first show that scVI scales to over one million cells, whereas competing algorithms can process at most tens of thousands of cells. Next, we show that scVI fits unseen data more closely and can impute missing data more accurately, both indicative of a better generalization capacity. We then utilize scVI to conduct a set of fundamental analysis tasks – including batch correction, visualization, clustering and differential expression – and demonstrate its accuracy in comparison to the state-of-the-art tools in each task. scVI is publicly available, and can be readily used as a principled and inclusive solution for multiple tasks of single-cell RNA sequencing data analysis.

Deciphering the potential of non-coding loci to influence gene regulation has been the subject of intense research, with important implications in understanding genetic underpinnings of human diseases. Massively parallel reporter assays (MPRAs) can measure regulatory activity of thousands of DNA-sequences and their variants in a single experiment. With increasing number of publically available MPRA datasets, one can now develop data-driven models which, given a DNA-sequence, predict its regulatory activity. Here, we performed a comprehensive meta-analysis of several MPRA datasets in a variety of cellular contexts. We first applied an ensemble of methods to predict MPRA output in each context and observed that the most predictive features are consistent across datasets. We then demonstrate that predictive models trained in one cellular context can be used to predict MPRA output in another, with loss of accuracy attributed to cell-type specific features. Finally, we show that our approach achieves top performance in the Fifth Critical Assessment of Genome Interpretation “Regulation Saturation” Challenge for predicting effects of single nucleotide variants. Overall, our analysis provides insights into how MPRA data can be leveraged to highlight functional regulatory regions throughout the genome and can guide effective design of future experiments by better prioritizing regions of interest.

The recent advent of methods for high-throughput single-cell molecular profiling has catalyzed a growing sense in the scientific community that the time is ripe to complete the 150-year-old effort to identify all cell types in the human body, by undertaking a Human Cell Atlas Project as an international collaborative effort. The aim would be to define all human cell types in terms of distinctive molecular profiles (e.g., gene expression) and connect this information with classical cellular descriptions (e.g., location and morphology). A comprehensive reference map of the molecular state of cells in healthy human tissues would propel the systematic study of physiological states, developmental trajectories, regulatory circuitry and interactions of cells, as well as provide a framework for understanding cellular dysregulation in human disease. Here we describe the idea, its potential utility, early proofs-of-concept, and some design considerations for the Human Cell Atlas.

The T cell compartment must contain diversity in both TCR repertoire and cell state to provide effective immunity against pathogens1,2. However, it remains unclear how differences in the TCR contribute to heterogeneity in T cell state at the single cell level because most analysis of the TCR repertoire has, to date, aggregated information from populations of cells. Single cell RNA-sequencing (scRNA-seq) can allow simultaneous measurement of TCR sequence and global transcriptional profile from single cells. However, current protocols to directly sequence the TCR require the use of long sequencing reads, increasing the cost and decreasing the number of cells that can be feasibly analyzed. Here we present a tool that can efficiently extract TCR sequence information from standard, short-read scRNA-seq libraries of T cells: TCR Reconstruction Algorithm for Paired-End Single cell (TRAPeS). We apply it to investigate heterogeneity in the CD8+T cell response in humans and mice, and show that it is accurate and more sensitive than previous approaches3,4. We applied TRAPeS to single cell RNA-seq of CD8+T cells specific for a single epitope from Yellow Fever Virus5. We show that the recently-described "naive-like" memory population of YFV-specific CD8+T cells have significantly longer CDR3 regions and greater divergence from germline sequence than do effector-memory phenotype CD8+T cells specific for YFV. This suggests that TCR usage contributes to heterogeneity in the differentiation state of the CD8+T cell response to YFV. TRAPeS is publicly available, and can be readily used to investigate the relationship between the TCR repertoire and cellular phenotype.