Single-Cell Multi-omic Integration Compares and Contrasts Features of Brain Cell Identity.

Defining cell types requires integrating diverse single-cell measurements from multiple experiments and biological contexts. To flexibly model single-cell datasets, we developed LIGER, an algorithm that delineates shared and dataset-specific features of cell identity. We applied it to four diverse and challenging analyses of human and mouse brain cells. First, we defined region-specific and sexually dimorphic gene expression in the mouse bed nucleus of the stria terminalis. Second, we analyzed expression in the human substantia nigra, comparing cell states in specific donors and relating cell types to those in the mouse. Third, we integrated in situ and single-cell expression data to spatially locate fine subtypes of cells present in the mouse frontal cortex. Finally, we jointly defined mouse cortical cell types using single-cell RNA-seq and DNA methylation profiles, revealing putative mechanisms of cell-type-specific epigenomic regulation. Integrative analyses using LIGER promise to accelerate investigations of cell-type definition, gene regulation, and disease states.

Endosteal stem cells at the bone-blood interface: A double-edged sword for rapid bone formation: Bone marrow endosteal stem cells provide a robust source of bone-making osteoblasts both in normal and abnormal bone formation.

Endosteal stem cells are a subclass of bone marrow skeletal stem cell populations that are particularly important for rapid bone formation occurring in growth and regeneration. These stem cells are strategically located near the bone surface in a specialized microenvironment of the endosteal niche. These stem cells are abundant in young stages but eventually depleted and replaced by other stem cell types residing in a non-endosteal perisinusoidal niche. Single-cell molecular profiling and in vivo cell lineage analyses play key roles in discovering endosteal stem cells. Importantly, endosteal stem cells can transform into bone tumor-making cells when deleterious mutations occur in tumor suppressor genes. The emerging hypothesis is that osteoblast-chondrocyte transitional identities confer a special subset of endosteal stromal cells with stem cell-like properties, which may make them susceptible for tumorigenic transformation. Endosteal stem cells are likely to represent an important therapeutic target of bone diseases caused by aberrant bone formation.

PyLiger: scalable single-cell multi-omic data integration in Python.

MOTIVATION: LIGER (Linked Inference of Genomic Experimental Relationships) is a widely used R package for single-cell multi-omic data integration. However, many users prefer to analyze their single-cell datasets in Python, which offers an attractive syntax and highly optimized scientific computing libraries for increased efficiency. RESULTS: We developed PyLiger, a Python package for integrating single-cell multi-omic datasets. PyLiger offers faster performance than the previous R implementation (2-5× speedup), interoperability with AnnData format, flexible on-disk or in-memory analysis capability and new functionality for gene ontology enrichment analysis. The on-disk capability enables analysis of arbitrarily large single-cell datasets using fixed memory. AVAILABILITY AND IMPLEMENTATION: PyLiger is available on Github at https://github.com/welch-lab/pyliger and on the Python Package Index. SUPPLEMENTARY INFORMATION: Supplementary data are available at Bioinformatics online.

Functional coordination of non-myocytes plays a key role in adult zebrafish heart regeneration.

Cardiac regeneration occurs primarily through proliferation of existing cardiomyocytes, but also involves complex interactions between distinct cardiac cell types including non-cardiomyocytes (non-CMs). However, the subpopulations, distinguishing molecular features, cellular functions, and intercellular interactions of non-CMs in heart regeneration remain largely unexplored. Using the LIGER algorithm, we assemble an atlas of cell states from 61,977 individual non-CM scRNA-seq profiles isolated at multiple time points during regeneration. This analysis reveals extensive non-CM cell diversity, including multiple macrophage (MC), fibroblast (FB), and endothelial cell (EC) subpopulations with unique spatiotemporal distributions, and suggests an important role for MC in inducing the activated FB and EC subpopulations. Indeed, pharmacological perturbation of MC function compromises the induction of the unique FB and EC subpopulations. Furthermore, we developed computational algorithm Topologizer to map the topological relationships and dynamic transitions between functional states. We uncover dynamic transitions between MC functional states and identify factors involved in mRNA processing and transcriptional regulation associated with the transition. Together, our single-cell transcriptomic analysis of non-CMs during cardiac regeneration provides a blueprint for interrogating the molecular and cellular basis of this process.

Single-cell transcriptomics reconstructs fate conversion from fibroblast to cardiomyocyte.

Direct lineage conversion offers a new strategy for tissue regeneration and disease modelling. Despite recent success in directly reprogramming fibroblasts into various cell types, the precise changes that occur as fibroblasts progressively convert to the target cell fates remain unclear. The inherent heterogeneity and asynchronous nature of the reprogramming process renders it difficult to study this process using bulk genomic techniques. Here we used single-cell RNA sequencing to overcome this limitation and analysed global transcriptome changes at early stages during the reprogramming of mouse fibroblasts into induced cardiomyocytes (iCMs). Using unsupervised dimensionality reduction and clustering algorithms, we identified molecularly distinct subpopulations of cells during reprogramming. We also constructed routes of iCM formation, and delineated the relationship between cell proliferation and iCM induction. Further analysis of global gene expression changes during reprogramming revealed unexpected downregulation of factors involved in mRNA processing and splicing. Detailed functional analysis of the top candidate splicing factor, Ptbp1, revealed that it is a critical barrier for the acquisition of cardiomyocyte-specific splicing patterns in fibroblasts. Concomitantly, Ptbp1 depletion promoted cardiac transcriptome acquisition and increased iCM reprogramming efficiency. Additional quantitative analysis of our dataset revealed a strong correlation between the expression of each reprogramming factor and the progress of individual cells through the reprogramming process, and led to the discovery of new surface markers for the enrichment of iCMs. In summary, our single-cell transcriptomics approaches enabled us to reconstruct the reprogramming trajectory and to uncover intermediate cell populations, gene pathways and regulators involved in iCM induction.

Deep sequencing shows multiple oligouridylations are required for 3' to 5' degradation of histone mRNAs on polyribosomes.

Histone mRNAs are rapidly degraded when DNA replication is inhibited during S phase with degradation initiating with oligouridylation of the stem loop at the 3' end. We developed a customized RNA sequencing strategy to identify the 3' termini of degradation intermediates of histone mRNAs. Using this strategy, we identified two types of oligouridylated degradation intermediates: RNAs ending at different sites of the 3' side of the stem loop that resulted from initial degradation by 3'hExo and intermediates near the stop codon and within the coding region. Sequencing of polyribosomal histone mRNAs revealed that degradation initiates and proceeds 3' to 5' on translating mRNA and that many intermediates are capped. Knockdown of the exosome-associated exonuclease PM/Scl-100, but not the Dis3L2 exonuclease, slows histone mRNA degradation consistent with 3' to 5' degradation by the exosome containing PM/Scl-100. Knockdown of No-go decay factors also slowed histone mRNA degradation, suggesting a role in removing ribosomes from partially degraded mRNAs.

TUT7 catalyzes the uridylation of the 3' end for rapid degradation of histone mRNA.

The replication-dependent histone mRNAs end in a stem-loop instead of the poly(A) tail present at the 3' end of all other cellular mRNAs. Following processing, the 3' end of histone mRNAs is trimmed to 3 nucleotides (nt) after the stem-loop, and this length is maintained by addition of nontemplated uridines if the mRNA is further trimmed by 3'hExo. These mRNAs are tightly cell-cycle regulated, and a critical regulatory step is rapid degradation of the histone mRNAs when DNA replication is inhibited. An initial step in histone mRNA degradation is digestion 2-4 nt into the stem by 3'hExo and uridylation of this intermediate. The mRNA is then subsequently degraded by the exosome, with stalled intermediates being uridylated. The enzyme(s) responsible for oligouridylation of histone mRNAs have not been definitively identified. Using high-throughput sequencing of histone mRNAs and degradation intermediates, we find that knockdown of TUT7 reduces both the uridylation at the 3' end as well as uridylation of the major degradation intermediate in the stem. In contrast, knockdown of TUT4 did not alter the uridylation pattern at the 3' end and had a small effect on uridylation in the stem-loop during histone mRNA degradation. Knockdown of 3'hExo also altered the uridylation of histone mRNAs, suggesting that TUT7 and 3'hExo function together in trimming and uridylating histone mRNAs.

Robust detection of alternative splicing in a population of single cells.

Single cell RNA-seq experiments provide valuable insight into cellular heterogeneity but suffer from low coverage, 3' bias and technical noise. These unique properties of single cell RNA-seq data make study of alternative splicing difficult, and thus most single cell studies have restricted analysis of transcriptome variation to the gene level. To address these limitations, we developed SingleSplice, which uses a statistical model to detect genes whose isoform usage shows biological variation significantly exceeding technical noise in a population of single cells. Importantly, SingleSplice is tailored to the unique demands of single cell analysis, detecting isoform usage differences without attempting to infer expression levels for full-length transcripts. Using data from spike-in transcripts, we found that our approach detects variation in isoform usage among single cells with high sensitivity and specificity. We also applied SingleSplice to data from mouse embryonic stem cells and discovered a set of genes that show significant biological variation in isoform usage across the set of cells. A subset of these isoform differences are linked to cell cycle stage, suggesting a novel connection between alternative splicing and the cell cycle.

A subset of replication-dependent histone mRNAs are expressed as polyadenylated RNAs in terminally differentiated tissues.

Histone proteins are synthesized in large amounts during S-phase to package the newly replicated DNA, and are among the most stable proteins in the cell. The replication-dependent (RD)-histone mRNAs expressed during S-phase end in a conserved stem-loop rather than a polyA tail. In addition, there are replication-independent (RI)-histone genes that encode histone variants as polyadenylated mRNAs. Most variants have specific functions in chromatin, but H3.3 also serves as a replacement histone for damaged histones in long-lived terminally differentiated cells. There are no reported replacement histone genes for histones H2A, H2B or H4. We report that a subset of RD-histone genes are expressed in terminally differentiated tissues as polyadenylated mRNAs, likely serving as replacement histone genes in long-lived non-dividing cells. Expression of two genes, HIST2H2AA3 and HIST1H2BC, is conserved in mammals. They are expressed as polyadenylated mRNAs in fibroblasts differentiated in vitro, but not in serum starved fibroblasts, suggesting that their expression is part of the terminal differentiation program. There are two histone H4 genes and an H3 gene that encode mRNAs that are polyadenylated and expressed at 5- to 10-fold lower levels than the mRNAs from H2A and H2B genes, which may be replacement genes for the H3.1 and H4 proteins.

Selective single cell isolation for genomics using microraft arrays.

Genomic methods are used increasingly to interrogate the individual cells that compose specific tissues. However, current methods for single cell isolation struggle to phenotypically differentiate specific cells in a heterogeneous population and rely primarily on the use of fluorescent markers. Many cellular phenotypes of interest are too complex to be measured by this approach, making it difficult to connect genotype and phenotype at the level of individual cells. Here we demonstrate that microraft arrays, which are arrays containing thousands of individual cell culture sites, can be used to select single cells based on a variety of phenotypes, such as cell surface markers, cell proliferation and drug response. We then show that a common genomic procedure, RNA-seq, can be readily adapted to the single cells isolated from these rafts. We show that data generated using microrafts and our modified RNA-seq protocol compared favorably with the Fluidigm C1. We then used microraft arrays to select pancreatic cancer cells that proliferate in spite of cytotoxic drug treatment. Our single cell RNA-seq data identified several expected and novel gene expression changes associated with early drug resistance.

A multiprotein occupancy map of the mRNP on the 3' end of histone mRNAs.

The animal replication-dependent (RD) histone mRNAs are coordinately regulated with chromosome replication. The RD-histone mRNAs are the only known cellular mRNAs that are not polyadenylated. Instead, the mature transcripts end in a conserved stem-loop (SL) structure. This SL structure interacts with the stem-loop binding protein (SLBP), which is involved in all aspects of RD-histone mRNA metabolism. We used several genomic methods, including high-throughput sequencing of cross-linked immunoprecipitate (HITS-CLIP) to analyze the RNA-binding landscape of SLBP. SLBP was not bound to any RNAs other than histone mRNAs. We performed bioinformatic analyses of the HITS-CLIP data that included (i) clustering genes by sequencing read coverage using CVCA, (ii) mapping the bound RNA fragment termini, and (iii) mapping cross-linking induced mutation sites (CIMS) using CLIP-PyL software. These analyses allowed us to identify specific sites of molecular contact between SLBP and its RD-histone mRNA ligands. We performed in vitro crosslinking assays to refine the CIMS mapping and found that uracils one and three in the loop of the histone mRNA SL preferentially crosslink to SLBP, whereas uracil two in the loop preferentially crosslinks to a separate component, likely the 3'hExo. We also performed a secondary analysis of an iCLIP data set to map UPF1 occupancy across the RD-histone mRNAs and found that UPF1 is bound adjacent to the SLBP-binding site. Multiple proteins likely bind the 3' end of RD-histone mRNAs together with SLBP.

EnD-Seq and AppEnD: sequencing 3' ends to identify nontemplated tails and degradation intermediates.

Existing methods for detecting RNA intermediates resulting from exonuclease degradation are low-throughput and laborious. In addition, mapping the 3' ends of RNA molecules to the genome after high-throughput sequencing is challenging, particularly if the 3' ends contain post-transcriptional modifications. To address these problems, we developed EnD-Seq, a high-throughput sequencing protocol that preserves the 3' end of RNA molecules, and AppEnD, a computational method for analyzing high-throughput sequencing data. Together these allow determination of the 3' ends of RNA molecules, including nontemplated additions. Applying EnD-Seq and AppEnD to histone mRNAs revealed that a significant fraction of cytoplasmic histone mRNAs end in one or two uridines, which have replaced the 1-2 nt at the 3' end of mature histone mRNA maintaining the length of the histone transcripts. Histone mRNAs in fly embryos and ovaries show the same pattern, but with different tail nucleotide compositions. We increase the sensitivity of EnD-Seq by using cDNA priming to specifically enrich low-abundance tails of known sequence composition allowing identification of degradation intermediates. In addition, we show the broad applicability of our computational approach by using AppEnD to gain insight into 3' additions from diverse types of sequencing data, including data from small capped RNA sequencing and some alternative polyadenylation protocols.

Pseudogenes transcribed in breast invasive carcinoma show subtype-specific expression and ceRNA potential.

BACKGROUND: Recent studies have shown that some pseudogenes are transcribed and contribute to cancer when dysregulated. In particular, pseudogene transcripts can function as competing endogenous RNAs (ceRNAs). The high similarity of gene and pseudogene nucleotide sequence has hindered experimental investigation of these mechanisms using RNA-seq. Furthermore, previous studies of pseudogenes in breast cancer have not integrated miRNA expression data in order to perform large-scale analysis of ceRNA potential. Thus, knowledge of both pseudogene ceRNA function and the role of pseudogene expression in cancer are restricted to isolated examples. RESULTS: To investigate whether transcribed pseudogenes play a pervasive regulatory role in cancer, we developed a novel bioinformatic method for measuring pseudogene transcription from RNA-seq data. We applied this method to 819 breast cancer samples from The Cancer Genome Atlas (TCGA) project. We then clustered the samples using pseudogene expression levels and integrated sample-paired pseudogene, gene and miRNA expression data with miRNA target prediction to determine whether more pseudogenes have ceRNA potential than expected by chance. CONCLUSIONS: Our analysis identifies with high confidence a set of 440 pseudogenes that are transcribed in breast cancer tissue. Of this set, 309 pseudogenes exhibit significant differential expression among breast cancer subtypes. Hierarchical clustering using only pseudogene expression levels accurately separates tumor samples from normal samples and discriminates the Basal subtype from the Luminal and Her2 subtypes. Correlation analysis shows more positively correlated pseudogene-parent gene pairs and negatively correlated pseudogene-miRNA pairs than expected by chance. Furthermore, 177 transcribed pseudogenes possess binding sites for co-expressed miRNAs that are also predicted to target their parent genes. Taken together, these results increase the catalog of putative pseudogene ceRNAs and suggest that pseudogene transcription in breast cancer may play a larger role than previously appreciated.

Integrating single-cell multimodal epigenomic data using 1D-convolutional neural networks.

Recent experimental developments enable single-cell multimodal epigenomic profiling, which measures multiple histone modifications and chromatin accessibility within the same cell. Such parallel measurements provide exciting new opportunities to investigate how epigenomic modalities vary together across cell types and states. A pivotal step in using this type of data is integrating the epigenomic modalities to learn a unified representation of each cell, but existing approaches are not designed to model the unique nature of this data type. Our key insight is to model single-cell multimodal epigenome data as a multi-channel sequential signal. Based on this insight, we developed ConvNet-VAEs, a novel framework that uses 1D-convolutional variational autoencoders (VAEs) for single-cell multimodal epigenomic data integration. We evaluated ConvNet-VAEs on nano-CT and scNTT-seq data generated from juvenile mouse brain and human bone marrow. We found that ConvNet-VAEs can perform dimension reduction and batch correction better than previous architectures while using significantly fewer parameters. Furthermore, the performance gap between convolutional and fully-connected architectures increases with the number of modalities, and deeper convolutional architectures can increase performance while performance degrades for deeper fully-connected architectures. Our results indicate that convolutional autoencoders are a promising method for integrating current and future single-cell multimodal epigenomic datasets.

HIV-1 Vpr combats the PU.1-driven antiviral response in primary human macrophages.

HIV-1 Vpr promotes efficient spread of HIV-1 from macrophages to T cells by transcriptionally downmodulating restriction factors that target HIV-1 Envelope protein (Env). Here we find that Vpr induces broad transcriptomic changes by targeting PU.1, a transcription factor necessary for expression of host innate immune response genes, including those that target Env. Consistent with this, we find silencing PU.1 in infected macrophages lacking Vpr rescues Env. Vpr downmodulates PU.1 through a proteasomal degradation pathway that depends on physical interactions with PU.1 and DCAF1, a component of the Cul4A E3 ubiquitin ligase. The capacity for Vpr to target PU.1 is highly conserved across primate lentiviruses. In addition to impacting infected cells, we find that Vpr suppresses expression of innate immune response genes in uninfected bystander cells, and that virion-associated Vpr can degrade PU.1. Together, we demonstrate Vpr counteracts PU.1 in macrophages to blunt antiviral immune responses and promote viral spread.

HIV-1 Vpr combats the PU.1-driven antiviral response in primary human macrophages.

HIV-1 Vpr promotes efficient spread of HIV-1 from macrophages to T cells by transcriptionally downmodulating restriction factors that target HIV-1 Envelope protein (Env). Here we find that Vpr induces broad transcriptomic changes by targeting PU.1, a transcription factor necessary for expression of host innate immune response genes, including those that target Env. Consistent with this, we find silencing PU.1 in infected macrophages lacking Vpr rescues Env. Vpr downmodulates PU.1 through a proteasomal degradation pathway that depends on physical interactions with PU.1 and DCAF1, a component of the Cul4A E3 ubiquitin ligase. The capacity for Vpr to target PU.1 is highly conserved across primate lentiviruses. In addition to impacting infected cells, we find that Vpr suppresses expression of innate immune response genes in uninfected bystander cells, and that virion-associated Vpr can degrade PU.1. Together, we demonstrate Vpr counteracts PU.1 in macrophages to blunt antiviral immune responses and promote viral spread.

Multi-omic single-cell velocity models epigenome-transcriptome interactions and improves cell fate prediction.

Multi-omic single-cell datasets, in which multiple molecular modalities are profiled within the same cell, offer an opportunity to understand the temporal relationship between epigenome and transcriptome. To realize this potential, we developed MultiVelo, a differential equation model of gene expression that extends the RNA velocity framework to incorporate epigenomic data. MultiVelo uses a probabilistic latent variable model to estimate the switch time and rate parameters of chromatin accessibility and gene expression and improves the accuracy of cell fate prediction compared to velocity estimates from RNA only. Application to multi-omic single-cell datasets from brain, skin and blood cells reveals two distinct classes of genes distinguished by whether chromatin closes before or after transcription ceases. We also find four types of cell states: two states in which epigenome and transcriptome are coupled and two distinct decoupled states. Finally, we identify time lags between transcription factor expression and binding site accessibility and between disease-associated SNP accessibility and expression of the linked genes. MultiVelo is available on PyPI, Bioconda and GitHub ( https://github.com/welch-lab/MultiVelo ).

Predicting the Structural Impact of Human Alternative Splicing.

Protein structure prediction with neural networks is a powerful new method for linking protein sequence, structure, and function, but structures have generally been predicted for only a single isoform of each gene, neglecting splice variants. To investigate the structural implications of alternative splicing, we used AlphaFold2 to predict the structures of more than 11,000 human isoforms. We employed multiple metrics to identify splicing-induced structural alterations, including template matching score, secondary structure composition, surface charge distribution, radius of gyration, accessibility of post-translational modification sites, and structure-based function prediction. We identified examples of how alternative splicing induced clear changes in each of these properties. Structural similarity between isoforms largely correlated with degree of sequence identity, but we identified a subset of isoforms with low structural similarity despite high sequence similarity. Exon skipping and alternative last exons tended to increase the surface charge and radius of gyration. Splicing also buried or exposed numerous post-translational modification sites, most notably among the isoforms of BAX. Functional prediction nominated numerous functional differences among isoforms of the same gene, with loss of function compared to the reference predominating. Finally, we used single-cell RNA-seq data from the Tabula Sapiens to determine the cell types in which each structure is expressed. Our work represents an important resource for studying the structure and function of splice isoforms across the cell types of the human body.

Bone marrow endosteal stem cells dictate active osteogenesis and aggressive tumorigenesis.

The bone marrow contains various populations of skeletal stem cells (SSCs) in the stromal compartment, which are important regulators of bone formation. It is well-described that leptin receptor (LepR)+ perivascular stromal cells provide a major source of bone-forming osteoblasts in adult and aged bone marrow. However, the identity of SSCs in young bone marrow and how they coordinate active bone formation remains unclear. Here we show that bone marrow endosteal SSCs are defined by fibroblast growth factor receptor 3 (Fgfr3) and osteoblast-chondrocyte transitional (OCT) identities with some characteristics of bone osteoblasts and chondrocytes. These Fgfr3-creER-marked endosteal stromal cells contribute to a stem cell fraction in young stages, which is later replaced by Lepr-cre-marked stromal cells in adult stages. Further, Fgfr3+ endosteal stromal cells give rise to aggressive osteosarcoma-like lesions upon loss of p53 tumor suppressor through unregulated self-renewal and aberrant osteogenic fates. Therefore, Fgfr3+ endosteal SSCs are abundant in young bone marrow and provide a robust source of osteoblasts, contributing to both normal and aberrant osteogenesis.

Population-scale skeletal muscle single-nucleus multi-omic profiling reveals extensive context specific genetic regulation.

Skeletal muscle, the largest human organ by weight, is relevant to several polygenic metabolic traits and diseases including type 2 diabetes (T2D). Identifying genetic mechanisms underlying these traits requires pinpointing the relevant cell types, regulatory elements, target genes, and causal variants. Here, we used genetic multiplexing to generate population-scale single nucleus (sn) chromatin accessibility (snATAC-seq) and transcriptome (snRNA-seq) maps across 287 frozen human skeletal muscle biopsies representing 456,880 nuclei. We identified 13 cell types that collectively represented 983,155 ATAC summits. We integrated genetic variation to discover 6,866 expression quantitative trait loci (eQTL) and 100,928 chromatin accessibility QTL (caQTL) (5% FDR) across the five most abundant cell types, cataloging caQTL peaks that atlas-level snATAC maps often miss. We identified 1,973 eGenes colocalized with caQTL and used mediation analyses to construct causal directional maps for chromatin accessibility and gene expression. 3,378 genome-wide association study (GWAS) signals across 43 relevant traits colocalized with sn-e/caQTL, 52% in a cell-specific manner. 77% of GWAS signals colocalized with caQTL and not eQTL, highlighting the critical importance of population-scale chromatin profiling for GWAS functional studies. GWAS-caQTL colocalization showed distinct cell-specific regulatory paradigms. For example, a C2CD4A/B T2D GWAS signal colocalized with caQTL in muscle fibers and multiple chromatin loop models nominated VPS13C, a glucose uptake gene. Sequence of the caQTL peak overlapping caSNP rs7163757 showed allelic regulatory activity differences in a human myocyte cell line massively parallel reporter assay. These results illuminate the genetic regulatory architecture of human skeletal muscle at high-resolution epigenomic, transcriptomic, and cell state scales and serve as a template for population-scale multi-omic mapping in complex tissues and traits.