Abstract Differential expression analysis of scRNA-seq data is central for characterizing how experimental factors affect the distribution of gene expression. However, it remains challenging to distinguish biological and technical sources of cell-cell variability and to assess the statistical significance of quantitative comparisons between groups of cells. We introduce memento to address these limitations and enable accurate and efficient differential expression analysis of the mean, variability, and gene correlation from scRNA-seq. We used memento to analyze 70,000 tracheal epithelial cells to identify interferon response genes with distinct variability and correlation patterns, 160,000 T cells perturbed with CRISPR-Cas9 to reconstruct gene-regulatory networks that control T cell activation, and 1.2 million PMBCs to map cell-type-specific cis expression quantitative trait loci (eQTLs). In all cases, memento identified more significant and reproducible differences in mean expression but also identified differences in variability and gene correlation that suggest distinct modes of transcriptional regulation imparted by cytokines, genetic perturbations, and natural genetic variation. These results demonstrate memento as a first-in-class method for the quantitative comparisons of scRNA-seq data scalable to millions of cells and thousands of samples.

Droplet-based single-cell RNA-sequencing (dscRNA-seq) has enabled rapid, massively parallel profiling of transcriptomes from tens of thousands of cells. Multiplexing samples for single cell capture and library preparation in dscRNA-seq would enable cost-effective designs of differential expression and genetic studies while avoiding technical batch effects, but its implementation remains challenging. Here, we introduce an in-silico algorithm demuxlet that harnesses natural genetic variation to discover the sample identity of each cell and identify droplets containing two cells. These capabilities enable multiplexed dscRNA-seq experiments where cells from unrelated individuals are pooled and captured at higher throughput than standard workflows. To demonstrate the performance of demuxlet, we sequenced 3 pools of peripheral blood mononuclear cells (PBMCs) from 8 lupus patients. Given genotyping data for each individual, demuxlet correctly recovered the sample identity of > 99% of singlets, and identified doublets at rates consistent with previous estimates. In PBMCs, we demonstrate the utility of multiplexed dscRNA-seq in two applications: characterizing cell type specificity and inter-individual variability of cytokine response from 8 lupus patients and mapping genetic variants associated with cell type specific gene expression from 23 donors. Demuxlet is fast, accurate, scalable and could be extended to other single cell datasets that incorporate natural or synthetic DNA barcodes.

Gene regulatory programs controlling the activation and polarization of CD4+ T cells are incompletely mapped and the interindividual variability in these programs remain unknown. We sequenced the transcriptomes of ~160k CD4+ T cells from 9 donors following pooled CRISPR perturbation targeting 140 regulators. We identified 134 regulators that affect T cell functionalization, including IRF2 as a positive regulator of Th2 polarization. Leveraging correlation patterns between cells, we mapped 194 pairs of interacting regulators, including known (e.g. BATF and JUN ) and novel interactions (e.g. ETS1 and STAT6 ). Finally, we identified 80 natural genetic variants with effects on gene expression, 48 of which are modified by a perturbation. In CD4+ T cells, CRISPR perturbations can influence in vitro polarization and modify the effects of trans and cis regulatory elements on gene expression.

The majority of genetic variants associated with common human diseases map to enhancers, non-coding elements that shape cell type-specific transcriptional programs and responses to specific extracellular cues1-3. In order to understand the mechanisms by which non-coding genetic variation contributes to disease, systematic mapping of functional enhancers and their biological contexts is required. Here, we develop an unbiased discovery platform that can identify enhancers for a target gene without prior knowledge of their native functional context. We used tiled CRISPR activation (CRISPRa) to synthetically recruit transcription factors to sites across large genomic regions (>100 kilobases) surrounding two key autoimmunity risk loci, CD69 and IL2RA (interleukin-2 receptor alpha; CD25). We identified several CRISPRa responsive elements (CaREs) with stimulation-dependent enhancer activity, including an IL2RA enhancer that harbors an autoimmunity risk variant. Using engineered mouse models and genome editing of human primary T cells, we found that sequence perturbation of the disease-associated IL2RA enhancer does not block IL2RA expression, but rather delays the timing of gene activation in response to specific extracellular signals. This work develops an approach to rapidly identify functional enhancers within non-coding regions, decodes a key human autoimmunity association, and suggests a general mechanism by which genetic variation can cause immune dysfunction.

Abstract Over 90% of genetic variants associated with complex human traits map to non-coding regions, but little is understood about how they modulate gene regulation in health and disease. One possible mechanism is that genetic variants affect the activity of one or more cis-regulatory elements leading to gene expression variation in specific cell types. To identify such cases, we analyzed Assay for Transposase-Accessible Chromatin sequencing (ATAC-seq) and RNA-seq profiles from activated CD4+ T cells of up to 105 healthy donors. We found that regions of accessible chromatin (ATAC-peaks) are co-accessible at kilobase and megabase resolution, in patterns consistent with the 3D organization of chromosomes measured by in situ Hi-C in T cells. 15% of genetic variants located within ATAC-peaks affected the accessibility of the corresponding peak through disrupting binding sites for transcription factors important for T cell differentiation and activation. These ATAC quantitative trait nucleotides (ATAC-QTNs) have the largest effects on co-accessible peaks, are associated with gene expression from the same aliquot of cells, are rarely affecting core binding motifs, and are enriched for autoimmune disease variants. Our results provide insights into how natural genetic variants modulate cis-regulatory elements, in isolation or in concert, to influence gene expression in primary immune cells that play a key role in many human diseases.

While the impact of common genetic variants on gene expression response to cellular stimuli has been analyzed in depth, less is known about how stimulation modulates the genetic control of isoform usage. Analyzing RNA-seq profiles of monocyte-derived dendritic cells from 243 individuals, we uncovered thousands of unannotated isoforms synthesized in response to viral infection and stimulation with type I interferon. We identified more than a thousand single nucleotide polymorphisms associated with isoform usage (isoQTLs), > 40% of which are independent of expression QTLs for the same gene. Compared to eQTLs, isoQTLs are enriched for splice sites and untranslated regions, and depleted of sequences upstream of annotated transcription start sites. Both eQTLs and isoQTLs in stimulated cells explain a significant proportion of the disease heritability attributed to common genetic variants. At the IRF7 locus, we found alternative promoter usage in response to influenza as a possible mechanism by which DNA variants previously associated with immune-related disorders mediate disease risk. At the ERAP2 locus, we shed light on the function of the major haplotype that has been maintained under long-term balancing selection. At baseline and following type 1 interferon stimulation, the major haplotype is associated with absence of ERAP2 expression while the minor haplotype, known to increase Crohn's disease risk, is associated with high ERAP2 expression. Surprisingly, in response to influenza infection, the major haplotype results in the expression of two uncharacterized, alternatively transcribed, spliced and translated short isoforms. Thus, genetic variants at a single locus could modulate independent gene regulatory processes in the innate immune response, and in the case of ERAP2, may confer a historical fitness advantage in response to virus.