Regulatory T cells (Treg) are conventionally viewed as suppressors of endogenous and therapy-induced antitumor immunity; however, their role in modulating responses to immune checkpoint blockade (ICB) is unclear. In this study, we integrated single-cell RNA-seq/T cell receptor sequencing (TCRseq) of >73,000 tumor-infiltrating Treg (TIL-Treg) from anti-PD-1-treated and treatment-naive non-small cell lung cancers (NSCLC) with single-cell analysis of tumor-associated antigen (TAA)-specific Treg derived from a murine tumor model. We identified 10 subsets of human TIL-Treg, most of which have high concordance with murine TIL-Treg subsets. Only one subset selectively expresses high levels of TNFRSF4 (OX40) and TNFRSF18 (GITR), whose engangement by cognate ligand mediated proliferative programs and NF-κB activation, as well as multiple genes involved in Treg suppression, including LAG3. Functionally, the OX40hiGITRhi subset is the most highly suppressive ex vivo, and its higher representation among total TIL-Treg correlated with resistance to PD-1 blockade. Unexpectedly, in the murine tumor model, we found that virtually all TIL-Treg-expressing T cell receptors that are specific for TAA fully develop a distinct TH1-like signature over a 2-week period after entry into the tumor, down-regulating FoxP3 and up-regulating expression of TBX21 (Tbet), IFNG, and certain proinflammatory granzymes. Transfer learning of a gene score from the murine TAA-specific TH1-like Treg subset to the human single-cell dataset revealed a highly analogous subcluster that was enriched in anti-PD-1-responding tumors. These findings demonstrate that TIL-Treg partition into multiple distinct transcriptionally defined subsets with potentially opposing effects on ICB-induced antitumor immunity and suggest that TAA-specific TIL-Treg may positively contribute to antitumor responses.

Pseudotime analysis with single-cell RNA-sequencing (scRNA-seq) data has been widely used to study dynamic gene regulatory programs along continuous biological processes. While many computational methods have been developed to infer the pseudo-temporal trajectories of cells within a biological sample, methods that compare pseudo-temporal patterns with multiple samples (or replicates) across different experimental conditions are lacking. Lamian is a comprehensive and statistically-rigorous computational framework for differential multi-sample pseudotime analysis. It can be used to identify changes in a biological process associated with sample covariates, such as different biological conditions, and also to detect changes in gene expression, cell density, and topology of a pseudotemporal trajectory. Unlike existing methods that ignore sample variability, Lamian draws statistical inference after accounting for cross-sample variability and hence substantially reduces sample-specific false discoveries that are not generalizable to new samples. Using both simulations and real scRNA-seq data, including an analysis of differential immune response programs between COVID-19 patients with different disease severity levels, we demonstrate the advantages of Lamian in decoding cellular gene expression programs in continuous biological processes.

ABSTRACT Promoting myelination capacity of endogenous oligodendrocyte precursor cells (OPCs) is a promising therapeutic approach for central nervous system demyelinating disorders such as Multiple Sclerosis (MS). To aid in the discovery of myelination promoting compounds, we generated an advanced, genome engineered, human pluripotent stem cell (hPSC) line that consist of three reporters (identification-and-purification tag, GFP, and secreted NanoLuc) driven by the endogenous PDGFRα , PLP1 and MBP genes, respectively. Based upon this line, we established a high-throughput drug screening platform and performed a small molecule screen with 2500 bioactive molecules. In addition to a number of previously known pathways, our screening effort identified new pathways whose inhibition enhance oligodendrocyte maturation and myelination. Although further genetic and molecular validation is required, the identified inhibitors could potentially be repurposed to develop remyelination therapy for MS and other demyelinating disorders.

Abstract Natural and induced somatic mutations that accumulate in the genome during development record the phylogenetic relationships of cells; however, whether these lineage barcodes can capture the dynamics of complex progenitor fields remains unclear. Here, we introduce quantitative fate mapping, an approach to simultaneously map the fate and quantify the commitment time, commitment bias, and population size of multiple progenitor groups during development based on a time-scaled phylogeny of their descendants. To reconstruct time-scaled phylogenies from lineage barcodes, we introduce Phylotime, a scalable maximum likelihood clustering approach based on a generalizable barcoding mutagenesis model. We validate these approaches using realistically-simulated barcoding results as well as experimental results from a barcoding stem cell line. We further establish criteria for the minimum number of cells that must be analyzed for robust quantitative fate mapping. Overall, this work demonstrates how lineage barcodes, natural or synthetic, can be used to obtain quantitative fate maps, thus enabling analysis of progenitor dynamics long after embryonic development in any organism.

Abstract Allele-specific expression, which measures the expression of two alleles of a gene in a diploid individual, is a powerful signal to study cis-regulatory effects. Comparing ASE across conditions, or differential ASE, can reveal context-specific gene regulation. Recently, single-cell RNA sequencing (scRNA-seq) has allowed the measurement of ASE at the resolution of individual cells, but there is a lack of statistical methods to analyze such data. We develop DAESC, a statistical method for differential ASE analysis across any condition of interest using scRNA-seq data from multiple individuals. DAESC includes a baseline model based on beta-binomial regression with random effects accounting for multiple cells from the same individual (DAESC-BB), and an extended mixture model that incorporates implicit haplotype phasing (DAESC-Mix). We demonstrate through simulations that DAESC accurately captures differential ASE effects in a wide range of scenarios. Application to scRNA-seq data from 105 induced pluripotent stem cell lines identifies 657 genes that are dynamically regulated during endoderm differentiation. A second application identifies several genes that are differentially regulated in pancreatic endocrine cells between type 2 diabetes patients and controls. In conclusion, DAESC is a powerful method for single-cell differential ASE analysis and can facilitate the discovery of context-specific regulatory effects.

Abstract Spatially resolved transcriptomics involves a set of emerging technologies that enable the transcriptomic profiling of tissues with the physical location of expressions. Although a variety of methods have been developed for data integration, most of them are for single-cell RNA-seq datasets without consideration of spatial information. Thus, methods that can integrate spatial transcriptomics data from multiple tissue slides, possibly from multiple individuals, are needed. Here, we present PRECAST, a data integration method for multiple spatial transcriptomics datasets with complex batch effects and/or biological effects between slides. PRECAST unifies spatial factor analysis simultaneously with spatial clustering and embedding alignment, while requiring only partially shared cell/domain clusters across datasets. Using both simulated and four real datasets, we show improved cell/domain detection with outstanding visualization, and the estimated aligned embeddings and cell/domain labels facilitate many downstream analyses. We demonstrate that PRECAST is computationally scalable and applicable to spatial transcriptomics datasets from different platforms.

Conventional high-throughput technologies for mapping regulatory element activities such as ChIP-seq, DNase-seq and FAIRE-seq cannot analyze samples with small number of cells. The recently developed ATAC-seq allows regulome mapping in small-cell-number samples, but its signal in single cell or samples with ≤500 cells remains discrete or noisy. Compared to these technologies, measuring transcriptome by RNA-seq in single-cell and small-cell-number samples is more mature. Here we show that one can globally predict chromatin accessibility and infer regulome using RNA-seq. Genome-wide chromatin accessibility predicted by RNA-seq from 30 cells is comparable with ATAC-seq from 500 cells. Predictions based on single-cell RNA-seq can more accurately reconstruct bulk chromatin accessibility than using single-cell ATAC-seq by pooling the same number of cells. Integrating ATAC-seq with predictions from RNA-seq increases power of both methods. Thus, transcriptome-based prediction can provide a new tool for decoding gene regulatory programs in small-cell-number samples.

Transcriptional responses to the Hedgehog (HH) signaling pathway are primarily modulated by GLI repression in the mouse limb. Previous studies suggested a role for the BAF chromatin remodeling complex in mediating GLI repression. Consistent with this possibility, the core BAF complex protein SMARCC1 is present at most active limb enhancers including the majority of GLI enhancers. However, in contrast to GLI repression which reduces chromatin accessibility, SMARCC1 maintains chromatin accessibility at most enhancers, including those bound by GLI. Moreover, SMARCC1 binding at GLI-regulated enhancers occurs independently of GLI3. Consistent with previous studies, some individual GLI target genes are mis-regulated in Smarcc1 conditional knockouts, though most GLI target genes are unaffected. Moreover, SMARCC1 is not necessary for mediating constitutive GLI repression in HH mutant limb buds. We conclude that SMARCC1 does not mediate GLI3 repression, which we propose utilizes alternative chromatin remodeling complexes.

Abstract Accurate alignment of transcribed RNA to reference genomes is a critical step in the analysis of gene expression, which in turn has broad applications in biomedical research and in the basic sciences. We have discovered that widely used splice-aware aligners, such as STAR and HISAT2, can introduce erroneous spliced alignments between repeated sequences, leading to the inclusion of falsely spliced transcripts in RNA-seq experiments. In some cases, the “phantom” introns resulting from these errors have made their way into widely-used genome annotation databases. To address this issue, we have developed EASTR (Emending Alignments of Spliced Transcript Reads), a novel software tool that can detect and remove falsely spliced alignments or transcripts from alignment and annotation files. EASTR improves the accuracy of spliced alignments across diverse species, including human, maize, and Arabidopsis thaliana , by detecting sequence similarity between intron-flanking regions. We demonstrate that applying EASTR before transcript assembly substantially reduces false positive introns, exons, and transcripts, improving the overall accuracy of assembled transcripts. Additionally, we show that EASTR’s application to reference annotation databases can detect and correct likely cases of mis-annotated transcripts.

Abstract The diversity of genetic programs and cellular plasticity of glioma-associated myeloid cells, and thus their contribution to tumor growth and immune evasion, is poorly understood. We performed single cell RNA-sequencing of immune and tumor cells from 33 glioma patients of varying tumor grades. We identified two populations characteristic of myeloid derived suppressor cells (MDSC), unique to glioblastoma (GBM) and absent in grades II and III tumors: i) an early progenitor population (E-MDSC) characterized by strong upregulation of multiple catabolic, anabolic, oxidative stress, and hypoxia pathways typically observed within tumor cells themselves, and ii) a monocytic MDSC (M-MDSC) population. The E-MDSCs geographically co-localize with a subset of highly metabolic glioma stem-like tumor cells with a mesenchymal program in the pseudopalisading region, a pathognomonic feature of GBMs associated with poor prognosis. Ligand-receptor interaction analysis revealed symbiotic cross-talk between the stemlike tumor cells and E-MDSCs in GBM, whereby glioma stem cells produce chemokines attracting E-MDSCs, which in turn produce growth and survival factors for the tumor cells. Our large-scale single-cell analysis elucidated unique MDSC populations as key facilitators of GBM progression and mediators of tumor immunosuppression, suggesting that targeting these specific myeloid compartments, including their metabolic programs, may be a promising therapeutic intervention in this deadly cancer. One-Sentence Summary Aggressive glioblastoma harbors two unique myeloid populations capable of promoting stem-like properties of tumor cells and suppressing T cell function in the tumor microenvironment.