Non-negative matrix factorization (NMF) is an unsupervised learning method well suited to high-throughput biology. Still, inferring biological processes requires additional post hoc statistics and annotation for interpretation of features learned from software packages developed for NMF implementation. Here, we aim to introduce a suite of computational tools that implement NMF and provide methods for accurate, clear biological interpretation and analysis. A generalized discussion of NMF covering its benefits, limitations, and open questions in the field is followed by three vignettes for the Bayesian NMF algorithm CoGAPS (Coordinated Gene Activity across Pattern Subsets). Each vignette will demonstrate NMF analysis to quantify cell state transitions in public domain single-cell RNA-sequencing (scRNA-seq) data of malignant epithelial cells in 25 pancreatic ductal adenocarcinoma (PDAC) tumors and 11 control samples. The first uses PyCoGAPS, our new Python interface for CoGAPS that we developed to enhance runtime of Bayesian NMF for large datasets. The second vignette steps through the same analysis using our R CoGAPS interface, and the third introduces two new cloud-based, plug-and-play options for running CoGAPS using GenePattern Notebook and Docker. By providing Python support, cloud-based computing options, and relevant example workflows, we facilitate user-friendly interpretation and implementation of NMF for single-cell analyses.


RNA velocity analysis of single cells promises to predict temporal dynamics from gene expression. Indeed, in many systems, it has been observed that RNA velocity produces a vector field that qualitatively reflects known features of the system. Despite this observation, the limitations of RNA velocity estimates are poorly understood. Using real data and simulations, we dissect the impact of different steps in the RNA velocity workflow on the estimated vector field. We find that the process of mapping RNA velocity estimates into a low-dimensional representation, such as those produced by UMAP, has a large impact on the result. The RNA velocity vector field strongly depends on the k-NN graph of the data. This dependence leads to significant estimator errors when the k-NN graph is not a faithful representation of the true data structure, a feature that cannot be known for most real datasets. Finally, we establish that RNA velocity estimates expression speed neither at the gene nor cellular level. We propose that RNA velocity is best considered a smoothed interpolation of the observed k-NN structure, as opposed to an extrapolation of future cellular states, and that the use of RNA velocity as a validation of latent space embedding structures is circular.


The receptor tyrosine kinase gene RET plays a critical role in the fate specification of enteric neural crest-derived cells (ENCDCs) during enteric nervous system (ENS) development. Pathogenic RET loss of function (LoF) alleles are associated with Hirschsprung disease (HSCR), which is marked by aganglionosis of the gastrointestinal (GI) tract. ENCDCs invade the developing GI tract, proliferate, migrate caudally, and differentiate into all of the major ENS cell types. Although the major phenotypic consequences and the underlying transcriptional changes from Ret LoF in the developing ENS have been described, its cell type and state-specific effects are unknown. Consequently, we performed single- cell RNA sequencing (scRNA-seq) on an enriched population of ENCDCs isolated from the developing GI tract of Ret null heterozygous and homozygous mouse embryos at embryonic day (E)12.5 and E14.5. We demonstrate four significant findings: (1) Ret-expressing ENCDCs are a heterogeneous population composed of ENS progenitors as well as glial and neuronal committed cells; (2) neurons committed to a predominantly inhibitory motor neuron developmental trajectory are not produced under Ret LoF, leaving behind a mostly excitatory motor neuron developmental program; (3) HSCR-associated and Ret gene regulatory network genes exhibit distinct expression patterns across Ret-expressing ENCDC with their expression impacted by Ret LoF; and (4) Ret deficiency leads to precocious differentiation and reduction in the number of proliferating ENS precursors. Our results support a model in which Ret contributes to multiple distinct cellular phenotypes associated with the proper development of the ENS, including the specification of inhibitory neuron subtypes, cell cycle dynamics of ENS progenitors, and the developmental timing of neuronal and glial commitment. Summary StatementRet LoF affects proper development of the mouse ENS through multiple distinct cellular phenotypes including restriction of neuronal fate potential, disruption of ENCDC migration, and modulation of progenitor proliferation rate.


Recent studies suggested that microglia, the primary brain immune cells, can affect circuit connectivity and neuronal function1-3. Microglia infiltrate the neuroepithelium early in embryonic development and are maintained in the brain throughout adulthood4,5. Several maternal environmental factors, such as aberrant microbiome, immune activation, and poor nutrition, can influence prenatal brain development6-8. Nevertheless, it is unknown how changes in the prenatal environment instruct the developmental trajectory of infiltrating microglia, which in turn affect brain development and function. Here we show that after maternal immune activation (MIA) microglia from the offspring have a long-lived decrease in immune reactivity (blunting) across the developmental trajectory. The blunted immune response was concomitant with changes in the chromatin accessibility and reduced transcription factor occupancy of the open chromatin. Single cell RNA sequencing revealed that MIA does not induce a distinct subpopulation but rather decreases the contribution to inflammatory microglia states. Prenatal replacement of MIA microglia with physiological infiltration of naive microglia ameliorated the immune blunting and restored a decrease in presynaptic vesicle release probability onto dopamine receptor type-two medium spiny neurons, indicating that aberrantly formed microglia due to an adverse prenatal environment impacts the long-term microglia reactivity and proper striatal circuit development.


Latent space techniques have emerged as powerful tools to identify genes and gene sets responsible for cell-type and species-specific differences in single-cell data. Transfer learning methods can compare learned latent spaces across biological systems. However, the robustness that comes from leveraging information across multiple genes in transfer learning is often attained at the sacrifice of gene-wise precision. Thus, methods are needed to identify genes, defined as important within a particular latent space, that significantly differ between contexts. To address this challenge, we have developed a new framework, scProject, and a new metric, projectionDrivers, to quantitatively examine latent space usage across single-cell experimental systems while concurrently extracting the genes driving the differential usage of the latent space between defined contrasts. Here, we demonstrate the efficacy, utility, and scalability of scProject with projectionDrivers and provide experimental validation for predicted species-specific differences between the developing mouse and human retina.


BackgroundThe cell cycle is a highly conserved, continuous process which controls faithful replication and division of cells. Single-cell technologies have enabled increasingly precise measurements of the cell cycle both as a biological process of interest and as a possible confounding factor. Despite its importance and conservation, there is no universally applicable approach to infer position in the cell cycle with high-resolution from single-cell RNA-seq data. ResultsHere, we present tricycle, an R/Bioconductor package, to address this challenge by leveraging key features of the biology of the cell cycle, the mathematical properties of principal component analysis of periodic functions, and the use of transfer learning. We estimate a cell cycle embedding using a fixed reference dataset and project new data into this reference embedding; an approach that overcomes key limitations of learning a dataset dependent embedding. Tricycle then predicts a cell-specific position in the cell cycle based on the data projection. The accuracy of tricycle compares favorably to gold-standard experimental assays, which generally require specialized measurements in specifically constructed in vitro systems. Using internal controls which are available for any dataset, we show that tricycle predictions generalize to datasets with multiple cell types, across tissues, species and even sequencing assays. ConclusionsTricycle generalizes across datasets, is highly scalable and applicable to atlas-level single-cell RNA-seq data.


Radial glial progenitor cells (RGCs) in the dorsal forebrain directly or indirectly produce excitatory projection neurons and macroglia of the neocortex. Recent evidence shows that the pool of RGCs is more heterogeneous than originally thought and that progenitor subpopulations can generate particular neuronal cell types. Using single cell RNA sequencing, we have studied gene expression patterns of two subtypes of RGCs that differ in their neurogenic behavior. One progenitor type rapidly produces postmitotic neurons, whereas the second progenitor remains relatively quiescence before generating neurons. We have identified candidate genes that are differentially expressed between these RGCs progenitor subtypes, including the transcription factor Sox9. Using in utero electroporation, we demonstrate that elevated Sox9 expression in progenitors prevents RGC division and leads to the generation of upper-layer cortical neurons from these progenitors at later ages. Our data thus reveal molecular differences between cortical progenitors with different neurogenic behavior and indicates that Sox9 is critical for the maintenance of RGCs to regulate the generation of upper layer neurons. SIGNIFICANCE STATEMENTThe existence of heterogeneity in the pool of RGCs and its relationship with the generation of cellular diversity in the cerebral cortex has been an interesting topic of debate for many years. Here we describe the existence of a subpopulation of RGCs with reduced neurogenic behavior at early embryonic ages presenting a particular molecular signature. This molecular signature consists of differential expression of some genes including the transcription factor Sox9, found to be a specific master regulator of this subpopulation of progenitor cells. Functional experiments perturbing Sox9 expressions levels reveal its instructive role in the regulation of the neurogenic behavior of RGCs and its relationship with the generation of upper layer projection neurons at later ages.


The enteric nervous system (ENS), a collection of neurons contained in the wall of the gut, is of fundamental importance to gastrointestinal and systemic health. According to the prevailing paradigm, the ENS arises from progenitor cells migrating from the embryonic neural crest and remains largely unchanged thereafter. Here, we show that the composition of maturing ENS changes with time, with a decline in neural-crest derived neurons and their replacement by mesoderm-derived neurons. Single cell transcriptomics and immunochemical approaches establish a distinct expression profile of mesoderm-derived neurons. The dynamic balance between the proportions of neurons from these two different lineages in the post-natal gut is dependent on the availability of their respective trophic signals, GDNF-RET and HGF-MET. With increasing age, the mesoderm-derived neurons become the dominant form of neurons in the ENS, a change associated with significant functional effects on intestinal motility. Normal intestinal function in the adult gastrointestinal tract therefore appears to require an optimal balance between these two distinct lineages within the ENS.


Mosquitoes locate and approach humans ( host-seek) when specific Olfactory Neurons (ORNs) in the olfactory periphery activate a specific combination of glomeruli in the mosquito Antennal Lobe (AL). We hypothesize that dysregulating proper glomerular activation in the presence of human odor will prevent host-seeking behavior. In experiments aimed at ectopically activating most ORNs in the presence of human odor, we made a surprising finding: ectopic expression of an AgOr (AgOr2) in Anopheles gambiae ORNs dampens the activity of the expressing neuron. This contrasts studies in Drosophila melanogaster, the typical insect model of olfaction, in which ectopic expression of non-native ORs in ORNs confers ectopic neuronal responses without interfering with native olfactory physiology. To gain insight into this dysfunction in mosquitoes, RNA-seq analyses were performed comparing wild-type antennae to those ectopically expressing AgOr2 in ORNs. Remarkably, almost all Or transcripts were significantly downregulated (except for AgOr2), and additional experiments suggest that it is AgOR2 protein rather than mRNA that mediates this downregulation. Our study shows that ORNs of Anopheles mosquitoes (in contrast to Drosophila) employ a currently unexplored regulatory mechanism of OR expression, which may be adaptable as a vector-control strategy. SIGNIFICANCE STATEMENTStudies in Drosophila melanogaster suggest that insect Olfactory Receptor Neurons (ORNs) do not contain mechanisms by which Odorant Receptors (ORs) regulate OR expression. This has proved useful in studies where ectopic expression of an OR in Drosophila ORNs confers responses to the odorants that activate the newly expressed OR. In experiments in Anopheles gambiae mosquitoes, we found that ectopic expression of an OR in most Anopheles ORNs dampened the activity of the expressing neurons. RNA-seq analyses demonstrated that ectopic OR expression in Anopheles ORNs leads to downregulation of endogenous Or transcripts. Additional experiments suggest that this downregulation required ectopic expression of a functional OR protein. These findings reveal that Anopheles mosquitoes, in contrast to Drosophila, contain a feedback mechanism to regulate OR expression. Mosquito ORNs might employ regulatory mechanisms of OR expression previously thought to occur only in non-insect olfactory systems.


Parallel processing circuits are thought to dramatically expand the network capabilities of the nervous system. Magnocellular and parvocellular oxytocin neurons have been proposed to subserve two parallel streams of social information processing, which allow a single molecule to encode a diverse array of ethologically distinct behaviors, although to date direct evidence to support this hypothesis is lacking. Here we provide the first comprehensive characterization of magnocellular and parvocellular oxytocin neurons, validated across anatomical, projection target, electrophysiological, and transcriptional criteria. We next used novel multiple feature selection tools in Fmr1 KO mice to provide direct evidence that normal functioning of the parvocellular but not magnocellular oxytocin pathway is required for autism-relevant social reward behavior. Finally, we demonstrate that autism risk genes are uniquely enriched in parvocellular oxytocin neurons. Taken together these results provide the first evidence that oxytocin pathway specific pathogenic mechanisms account for social impairments across a broad range of autism etiologies. One Sentence SummaryPathoclisis of parvocellular oxytocin neurons plays an important role in the pathogenesis of social impairments in autism.


The development of single-cell RNA-Sequencing (scRNA-Seq) has allowed high resolution analysis of cell type diversity and transcriptional networks controlling cell fate specification. To identify the transcriptional networks governing human retinal development, we performed scRNA-Seq over retinal organoid and in vivo retinal development, across 20 timepoints. Using both pseudotemporal and cross-species analyses, we examined the conservation of gene expression across retinal progenitor maturation and specification of all seven major retinal cell types. Furthermore, we examined gene expression differences between developing macula and periphery and between two distinct populations of horizontal cells. We also identify both shared and species-specific patterns of gene expression during human and mouse retinal development. Finally, we identify an unexpected role for ATOH7 expression in regulation of photoreceptor specification during late retinogenesis. These results provide a roadmap to future studies of human retinal development, and may help guide the design of cell-based therapies for treating retinal dystrophies.


MotivationDimension reduction techniques are widely used to interpret high-dimensional biological data. Features learned from these methods are used to discover both technical artifacts and novel biological phenomena. Such feature discovery is critically import to large single-cell datasets, where lack of a ground truth limits validation and interpretation. Transfer learning (TL) can be used to relate the features learned from one source dataset to a new target dataset to perform biologically-driven validation by evaluating their use in or association with additional sample annotations in that independent target dataset.\n\nResultsWe developed an R/Bioconductor package, projectR, to perform TL for analyses of genomics data via TL of clustering, correlation, and factorization methods. We then demonstrate the utility TL for integrated data analysis with an example for spatial single-cell analysis.\n\nAvailabilityprojectR is available on Bioconductor and at\n\;


Chromatin modifiers act to coordinate gene expression changes critical to neuronal differentiation from neural stem/progenitor cells (NSPCs). Lysine-specific methyltransferase 2D (KMT2D) encodes a histone methyltransferase that promotes transcriptional activation, and is frequently mutated in cancers and in the majority (>70%) of patients diagnosed with the congenital, multisystem intellectual disability (ID) disorder Kabuki syndrome 1 (KS1). Critical roles for KMT2D are established in various non-neural tissues, but the effects of KMT2D loss in brain cell development have not been described. We conducted parallel studies of proliferation, differentiation, transcription, and chromatin profiling in KMT2D-deficient human and mouse models to define KMT2D-regulated functions in neurodevelopmental contexts, including adult-born hippocampal NSPCs in vivo and in vitro. We report cell-autonomous defects in proliferation, cell cycle, and survival, accompanied by early NSPC maturation in several KMT2D-deficient model systems. Transcriptional suppression in KMT2D-deficient cells indicated strong perturbation of hypoxia-responsive metabolism pathways. Functional experiments confirmed abnormalities of cellular hypoxia responses in KMT2D-deficient neural cells, and accelerated NSPC maturation in vivo. Together, our findings support a model in which loss of KMT2D function suppresses expression of oxygen-responsive gene programs important to neural progenitor maintenance, resulting in precocious neuronal differentiation in a mouse model of KS1.\n\nGraphical Abstract\n\nO_FIG O_LINKSMALLFIG WIDTH=200 HEIGHT=91 SRC=\"FIGDIR/small/484410v3_ufig1.gif\" ALT=\"Figure 1\">\nView larger version (22K):\norg.highwire.dtl.DTLVardef@55dd1aorg.highwire.dtl.DTLVardef@1270c80org.highwire.dtl.DTLVardef@a5bda2org.highwire.dtl.DTLVardef@144df37_HPS_FORMAT_FIGEXP M_FIG C_FIG


Recent genome-wide association studies (GWAS) identified numerous schizophrenia (SZ) and Alzheimers disease (AD) associated loci, most outside protein-coding regions and hypothesized to affect gene transcription. We used a massively parallel reporter assay (MPRA) to screen, 1,049 SZ and 30 AD variants in 64 and 9 loci respectively for allele differences in driving reporter gene expression. A library of synthetic oligonucleotides assaying each allele 5 times was transfected into K562 chronic myelogenous leukemia lymphoblasts and SK-SY5Y human neuroblastoma cells. 148 variants showed allelic differences in K562 and 53 in SK-SY5Y cells, on average 2.6 variants per locus. Nine showed significant differences in both lines, a modest overlap reflecting different regulatory landscapes of these lines that also differ significantly in chromatin marks. Eight of nine were in the same direction. We observe no preference for risk alleles to increase or decrease expression. We find a positive correlation between the number of SNPs in Linkage Disequilibrium (LD) and the proportion of functional SNPs supporting combinatorial effects that may lead to haplotype selection. Our results prioritize future functional follow up of disease associated SNPs to determine the driver GWAS variant(s), at each locus and enhance our understanding of gene regulation dynamics.


Tumor heterogeneity provides a complex challenge to cancer treatment and is a critical component of therapeutic response, disease recurrence, and patient survival. Single-cell RNA-sequencing (scRNA-seq) technologies reveal the prevalence of intra-and inter-tumor heterogeneity. Computational techniques are essential to quantify the differences in variation of these profiles between distinct cell types, tumor subtypes, and patients to fully characterize intra-and inter-tumor molecular heterogeneity. We devised a new algorithm, Expression Variation Analysis in Single Cells (EVAsc), to perform multivariate statistical analyses of differential variation of expression in gene sets for scRNA-seq. EVAsc has high sensitivity and specificity to detect pathways with true differential heterogeneity in simulated data. We then apply EVAsc to several public domain scRNA-seq tumor datasets to quantify the landscape of tumor heterogeneity in several key applications in cancer genomics, i.e. immunogenicity, cancer subtypes, and metastasis. Immune pathway heterogeneity in hematopoietic cell populations in breast tumors corresponded to the amount diversity present in the T-cell repertoire of each individual. In head and neck squamous cell carcinoma (HNSCC) patients, we found dramatic differences in pathway dysregulation across basal primary tumors. Within the basal primary tumors we also identified increased immune dysregulation in individuals with a high proportion of fibroblasts present in the tumor microenvironment. Moreover, cells in HNSCC primary tumors had significantly more heterogeneity across pathways than cells in metastases, consistent with a model of clonal outgrowth. These results demonstrate the broad utility of EVAsc to quantify inter-and intra-tumor heterogeneity from scRNA-seq data without reliance on low dimensional visualization.


New approaches are urgently needed to glean biological insights from the vast amounts of single cell RNA sequencing (scRNA-Seq) data now being generated. To this end, we propose that cell identity should map to a reduced set of factors which will describe both exclusive and shared biology of individual cells, and that the dimensions which contain these factors reflect biologically meaningful relationships across different platforms, tissues and species. To find a robust set of dependent factors in large-scale scRNA- Seq data, we developed a Bayesian non-negative matrix factorization (NMF) algorithm, scCoGAPS. Application of scCoGAPS to scRNA-Seq data obtained over the course of mouse retinal development identified gene expression signatures for factors associated with specific cell types and continuous biological processes. To test whether these signatures are shared across diverse cellular contexts, we developed projectR to map biologically disparate datasets into the factors learned by scCoGAPS. Because projecting these dimensions preserve relative distances between samples, biologically meaningful relationships/factors will stratify new data consistent with their underlying processes, allowing labels or information from one dataset to be used for annotation of the other--a machine learning concept called transfer learning. Using projectR, data from multiple datasets was used to annotate latent spaces and reveal novel parallels between developmental programs in other tissues, species and cellular assays. Using this approach we are able to transfer cell type and state designations across datasets to rapidly annotate cellular features in a new dataset without a priori knowledge of their type, identify a species-specific signature of microglial cells, and identify a previously undescribed subpopulation of neurosecretory cells within the lung. Together, these algorithms define biologically meaningful dimensions of cellular identity, state, and trajectories that persist across technologies, molecular features, and species.\n\nGRAPHICAL ABSTRACT\n\nO_FIG O_LINKSMALLFIG WIDTH=174 HEIGHT=200 SRC=\"FIGDIR/small/395004_ufig1.gif\" ALT=\"Figure 1\">\nView larger version (81K):\norg.highwire.dtl.DTLVardef@dd1c07org.highwire.dtl.DTLVardef@5b1109org.highwire.dtl.DTLVardef@bb6714org.highwire.dtl.DTLVardef@16c66f0_HPS_FORMAT_FIGEXP M_FIG C_FIG


Precise temporal control of gene expression in neuronal progenitors is necessary for correct regulation of neurogenesis and cell fate specification. However, the extensive cellular heterogeneity of the developing CNS has posed a major obstacle to identifying the gene regulatory networks that control these processes. To address this, we used single cell RNA-sequencing to profile ten developmental stages encompassing the full course of retinal neurogenesis. This allowed us to comprehensively characterize changes in gene expression that occur during initiation of neurogenesis, changes in developmental competence, and specification and differentiation of each of the major retinal cell types. These data identify transitions in gene expression between early and late-stage retinal progenitors, as well as a classification of neurogenic progenitors. We identify here the NFI family of transcription factors (Nfia, Nfib, and Nfix) as genes with enriched expression within late RPCs, and show they are regulators of bipolar interneuron and Muller glia specification and the control of proliferative quiescence.


Omics data contains signal from the molecular, physical, and kinetic inter- and intra-cellular interactions that control biological systems. Matrix factorization techniques can reveal low-dimensional structure from high-dimensional data that reflect these interactions. These techniques can uncover new biological knowledge from diverse high-throughput omics data in topics ranging from pathway discovery to time course analysis. We review exemplary applications of matrix factorization for systems-level analyses. We discuss appropriate application of these methods, their limitations, and focus on analysis of results to facilitate optimal biological interpretation. The inference of biologically relevant features with matrix factorization enables discovery from high-throughput data beyond the limits of current biological knowledge--answering questions from high-dimensional data that we have not yet thought to ask.


Massively parallel reporter assays (MPRAs) have emerged as a popular means for understanding noncoding variation in a variety of conditions. While a large number of experiments have been described in the literature, analysis typically uses ad-hoc methods. There has been little attention to comparing performance of methods across datasets. We present the mpralm method which we show is calibrated and powerful, by analyzing its performance on multiple MPRA datasets. We show that it outperforms existing statistical methods for analysis of this data type, in the first comprehensive evaluation of statistical methods on several datasets. We investigate theoretical and real-data properties of barcode summarization methods and show an unappreciated impact of summarization method for some datasets. Finally, we use our model to conduct a power analysis for this assay and show substantial improvements in power by performing up to 6 replicates per condition, whereas sequencing depth has smaller impact; we recommend to always use at least 4 replicates. Together, these results inform recommendations for differential analysis, general group comparisons, and power analysis and will help improve design and analysis of MPRA experiments. An R package is available from the Bioconductor project at


Parkinsons disease (PD) is caused by the collapse of substantia nigra (SN) dopaminergic (DA) neurons of the midbrain (MB), while other DA populations remain relatively intact. Common variation influencing susceptibility to sporadic PD has been primarily identified through genome wide association studies (GWAS). However, like many other common genetic diseases, the genes impacted by common PD-associated variation remain to be elucidated. Here, we used single-cell RNA-seq to characterize DA neuron populations in the mouse brain at embryonic and early postnatal timepoints. These data allow for the unbiased identification of DA neuron subpopulations, including a novel postnatal neuroblast population and SN DA neurons. Comparison of SN DA neurons with other DA neurons populations in the brain reveals a unique transcriptional profile, novel marker genes, and specific gene regulatory networks. By integrating these cell population specific data with published GWAS, we develop a scoring system for prioritizing candidate genes in PD-associated loci. With this, we prioritize candidate genes in all 32 GWAS intervals implicated in sporadic PD risk, the first such systematically generated list. From this we confirm that the prioritized candidate gene CPLX1 disrupts the nigrostriatal pathway when knocked out in mice. Ultimately, this systematic rationale leads to the identification of biologically pertinent candidates and testable hypotheses for sporadic PD that will inform a new era of PD genetic research.


Single-cell RNA sequencing technologies have generated the first catalogs of transcriptionally defined neuronal subtypes of the brain. However, the biologically informative cellular processes that contribute to neuronal subtype specification and transcriptional heterogeneity remain unclear. By comparing the gene expression profiles of single layer 6 corticothalamic neurons in somatosensory cortex, we show that transcriptional subtypes primarily reflect axonal projection pattern, laminar position within the cortex, and neuronal activity state. Pseudotemporal ordering of 1023 cellular responses to manipulation of sensory input demonstrates that changes in expression of activity-induced genes both reinforced cell-type identity and contributed to increased transcriptional heterogeneity within each cell type. This is due to cell-type specific biases in the choice of transcriptional states following manipulation of neuronal activity. These results reveal that axonal projection pattern, laminar position, and activity state define significant axes of variation that contribute both to the transcriptional identity of individual neurons and to the transcriptional heterogeneity within each neuronal subtype.