Demultiplexing of single-cell RNA-sequencing data using interindividual variation in gene expression.

Motivation: Pooled designs for single-cell RNA sequencing, where many cells from distinct samples are processed jointly, offer increased throughput and reduced batch variation. This study describes expression-aware demultiplexing (EAD), a computational method that employs differential co-expression patterns between individuals to demultiplex pooled samples without any extra experimental steps. Results: We use synthetic sample pools and show that the top interindividual differentially co-expressed genes provide a distinct cluster of cells per individual, significantly enriching the regulation of metabolism. Our application of EAD to samples of six isogenic inbred mice demonstrated that controlling genetic and environmental effects can solve interindividual variations related to metabolic pathways. We utilized 30 samples from both sepsis and healthy individuals in six batches to assess the performance of classification approaches. The results indicate that combining genetic and EAD results can enhance the accuracy of assignments (Min. 0.94, Mean 0.98, Max. 1). The results were enhanced by an average of 1.4% when EAD and barcoding techniques were combined (Min. 1.25%, Median 1.33%, Max. 1.74%). Furthermore, we demonstrate that interindividual differential co-expression analysis within the same cell type can be used to identify cells from the same donor in different activation states. By analysing single-nuclei transcriptome profiles from the brain, we demonstrate that our method can be applied to nonimmune cells. Availability and implementation: EAD workflow is available at https://isarnassiri.github.io/scDIV/ as an R package called scDIV (acronym for single-cell RNA-sequencing data demultiplexing using interindividual variations).

An isoform quantitative trait locus in SBNO2 links genetic susceptibility to Crohn's disease with defective antimicrobial activity.

Despite major advances in linking single genetic variants to single causal genes, the significance of genetic variation on transcript-level regulation of expression, transcript-specific functions, and relevance to human disease has been poorly investigated. Strawberry notch homolog 2 (SBNO2) is a candidate gene in a susceptibility locus with different variants associated with Crohn's disease and bone mineral density. The SBNO2 locus is also differentially methylated in Crohn's disease but the functional mechanisms are unknown. Here we show that the isoforms of SBNO2 are differentially regulated by lipopolysaccharide and IL-10. We identify Crohn's disease associated isoform quantitative trait loci that negatively regulate the expression of the noncanonical isoform 2 corresponding with the methylation signals at the isoform 2 promoter in IBD and CD. The two isoforms of SBNO2 drive differential gene networks with isoform 2 dominantly impacting antimicrobial activity in macrophages. Our data highlight the role of isoform quantitative trait loci to understand disease susceptibility and resolve underlying mechanisms of disease.

Characterization of the genetic determinants of context-specific DNA methylation in primary monocytes.

To better understand inter-individual variation in sensitivity of DNA methylation (DNAm) to immune activity, we characterized effects of inflammatory stimuli on primary monocyte DNAm (n = 190). We find that monocyte DNAm is site-dependently sensitive to lipopolysaccharide (LPS), with LPS-induced demethylation occurring following hydroxymethylation. We identify 7,359 high-confidence immune-modulated CpGs (imCpGs) that differ in genomic localization and transcription factor usage according to whether they represent a gain or loss in DNAm. Demethylated imCpGs are profoundly enriched for enhancers and colocalize to genes enriched for disease associations, especially cancer. DNAm is age associated, and we find that 24-h LPS exposure triggers approximately 6 months of gain in epigenetic age, directly linking epigenetic aging with innate immune activity. By integrating LPS-induced changes in DNAm with genetic variation, we identify 234 imCpGs under local genetic control. Exploring shared causal loci between LPS-induced DNAm responses and human disease traits highlights examples of disease-associated loci that modulate imCpG formation.

A role for SETD2 loss in tumorigenesis through DNA methylation dysregulation.

SETD2-dependent H3 Lysine-36 trimethylation (H3K36me3) has been recently linked to the deposition of de-novo DNA methylation. SETD2 is frequently mutated in cancer, however, the functional impact of SETD2 loss and depletion on DNA methylation across cancer types and tumorigenesis is currently unknown. Here, we perform a pan-cancer analysis and show that both SETD2 mutation and reduced expression are associated with DNA methylation dysregulation across 21 out of the 24 cancer types tested. In renal cancer, these DNA methylation changes are associated with altered gene expression of oncogenes, tumour suppressors, and genes involved in neoplasm invasiveness, including TP53, FOXO1, and CDK4. This suggests a new role for SETD2 loss in tumorigenesis and cancer aggressiveness through DNA methylation dysregulation. Moreover, using a robust machine learning methodology, we develop and validate a 3-CpG methylation signature which is sufficient to predict SETD2 mutation status with high accuracy and correlates with patient prognosis.

scQCEA: a framework for annotation and quality control report of single-cell RNA-sequencing data.

BACKGROUND: Systematic description of library quality and sequencing performance of single-cell RNA sequencing (scRNA-seq) data is imperative for subsequent downstream modules, including re-pooling libraries. While several packages have been developed to visualise quality control (QC) metrics for scRNA-seq data, they do not include expression-based QC to discriminate between true variation and background noise. RESULTS: We present scQCEA (acronym of the single-cell RNA sequencing Quality Control and Enrichment Analysis), an R package to generate reports of process optimisation metrics for comparing sets of samples and visual evaluation of quality scores. scQCEA can import data from 10X or other single-cell platforms and includes functions for generating an interactive report of QC metrics for multi-omics data. In addition, scQCEA provides automated cell type annotation on scRNA-seq data using differential gene expression patterns for expression-based quality control. We provide a repository of reference gene sets, including 2348 marker genes, which are exclusively expressed in 95 human and mouse cell types. Using scRNA-seq data from 56 gene expressions and V(D)J T cell replicates, we show how scQCEA can be applied for the visual evaluation of quality scores for sets of samples. In addition, we use the summary of QC measures from 342 human and mouse shallow-sequenced gene expression profiles to specify optimal sequencing requirements to run a cell-type enrichment analysis function. CONCLUSIONS: The open-source R tool will allow examining biases and outliers over biological and technical measures, and objective selection of optimal cluster numbers before downstream analysis. scQCEA is available at https://isarnassiri.github.io/scQCEA/ as an R package. Full documentation, including an example, is provided on the package website.

GWAS and meta-analysis identifies 49 genetic variants underlying critical COVID-19.

Critical illness in COVID-19 is an extreme and clinically homogeneous disease phenotype that we have previously shown1 to be highly efficient for discovery of genetic associations2. Despite the advanced stage of illness at presentation, we have shown that host genetics in patients who are critically ill with COVID-19 can identify immunomodulatory therapies with strong beneficial effects in this group3. Here we analyse 24,202 cases of COVID-19 with critical illness comprising a combination of microarray genotype and whole-genome sequencing data from cases of critical illness in the international GenOMICC (11,440 cases) study, combined with other studies recruiting hospitalized patients with a strong focus on severe and critical disease: ISARIC4C (676 cases) and the SCOURGE consortium (5,934 cases). To put these results in the context of existing work, we conduct a meta-analysis of the new GenOMICC genome-wide association study (GWAS) results with previously published data. We find 49 genome-wide significant associations, of which 16 have not been reported previously. To investigate the therapeutic implications of these findings, we infer the structural consequences of protein-coding variants, and combine our GWAS results with gene expression data using a monocyte transcriptome-wide association study (TWAS) model, as well as gene and protein expression using Mendelian randomization. We identify potentially druggable targets in multiple systems, including inflammatory signalling (JAK1), monocyte-macrophage activation and endothelial permeability (PDE4A), immunometabolism (SLC2A5 and AK5), and host factors required for viral entry and replication (TMPRSS2 and RAB2A).

IL7 genetic variation and toxicity to immune checkpoint blockade in patients with melanoma.

Treatment with immune checkpoint blockade (ICB) frequently triggers immune-related adverse events (irAEs), causing considerable morbidity. In 214 patients receiving ICB for melanoma, we observed increased severe irAE risk in minor allele carriers of rs16906115, intronic to IL7. We found that rs16906115 forms a B cell-specific expression quantitative trait locus (eQTL) to IL7 in patients. Patients carrying the risk allele demonstrate increased pre-treatment B cell IL7 expression, which independently associates with irAE risk, divergent immunoglobulin expression and more B cell receptor mutations. Consistent with the role of IL-7 in T cell development, risk allele carriers have distinct ICB-induced CD8+ T cell subset responses, skewing of T cell clonality and greater proportional repertoire occupancy by large clones. Finally, analysis of TCGA data suggests that risk allele carriers independently have improved melanoma survival. These observations highlight key roles for B cells and IL-7 in both ICB response and toxicity and clinical outcomes in melanoma.

An immunodominant NP105-113-B*07:02 cytotoxic T cell response controls viral replication and is associated with less severe COVID-19 disease.

NP105-113-B*07:02-specific CD8+ T cell responses are considered among the most dominant in SARS-CoV-2-infected individuals. We found strong association of this response with mild disease. Analysis of NP105-113-B*07:02-specific T cell clones and single-cell sequencing were performed concurrently, with functional avidity and antiviral efficacy assessed using an in vitro SARS-CoV-2 infection system, and were correlated with T cell receptor usage, transcriptome signature and disease severity (acute n = 77, convalescent n = 52). We demonstrated a beneficial association of NP105-113-B*07:02-specific T cells in COVID-19 disease progression, linked with expansion of T cell precursors, high functional avidity and antiviral effector function. Broad immune memory pools were narrowed postinfection but NP105-113-B*07:02-specific T cells were maintained 6 months after infection with preserved antiviral efficacy to the SARS-CoV-2 Victoria strain, as well as Alpha, Beta, Gamma and Delta variants. Our data show that NP105-113-B*07:02-specific T cell responses associate with mild disease and high antiviral efficacy, pointing to inclusion for future vaccine design.

Immune checkpoint blockade sensitivity and progression-free survival associates with baseline CD8+ T cell clone size and cytotoxicity.

The antitumor action of immune checkpoint blockade (ICB) is primarily mediated by CD8+ T cells. How sensitivity to ICB varies across CD8+ T cell subsets and clonotypes and the relationship of these with clinical outcome is unclear. To explore this, we used single-cell V(D)J and RNA-sequencing to track gene expression changes elicited by ICB across individual peripheral CD8+ T cell clones, identify baseline markers of CD8+ T cell clonal sensitivity, and chart how CD8+ T cell transcriptional changes vary according to phenotypic subset and clonal size. We identified seven subsets of CD8+ T cells with divergent reactivity to ICB and found that the cytotoxic effector subset showed the greatest number of differentially expressed genes while remaining stable in clonal size after ICB. At the level of CD8+ T cell clonotypes, we found a relationship between transcriptional changes and clone size, with large clones showing a greater number of differentially regulated genes enriched for pathways including T cell receptor (TCR) signaling. Cytotoxic CD8+ effector clones were more likely to persist following ICB and were more likely to correspond with public tumor-infiltrating lymphocyte clonotypes. Last, we demonstrated that individuals whose CD8+ T cell pretreatment showed low cytotoxicity and had fewer expanded clones typically had worse outcomes after ICB treatment. This work further advances understanding of the molecular determinants of ICB response, assisting in the search for peripheral prognostic biomarkers and highlighting the importance of the baseline CD8+ immune landscape in determining ICB response in metastatic melanoma.

Interferon-Gamma-Producing CD8+ Tissue Resident Memory T Cells Are a Targetable Hallmark of Immune Checkpoint Inhibitor-Colitis.

BACKGROUND & AIMS: The pathogenesis of immune checkpoint inhibitor (ICI)-colitis remains incompletely understood. We sought to identify key cellular drivers of ICI-colitis and their similarities to idiopathic ulcerative colitis, and to determine potential novel therapeutic targets. METHODS: We used a cross-sectional approach to study patients with ICI-colitis, those receiving ICI without the development of colitis, idiopathic ulcerative colitis, and healthy controls. A subset of patients with ICI-colitis were studied longitudinally. We applied a range of methods, including multiparameter and spectral flow cytometry, spectral immunofluorescence microscopy, targeted gene panels, and bulk and single-cell RNA sequencing. RESULTS: We demonstrate CD8+ tissue resident memory T (TRM) cells are the dominant activated T cell subset in ICI-colitis. The pattern of gastrointestinal immunopathology is distinct from ulcerative colitis at both the immune and epithelial-signaling levels. CD8+ TRM cell activation correlates with clinical and endoscopic ICI-colitis severity. Single-cell RNA sequencing analysis confirms activated CD8+ TRM cells express high levels of transcripts for checkpoint inhibitors and interferon-gamma in ICI-colitis. We demonstrate similar findings in both anti-CTLA-4/PD-1 combination therapy and in anti-PD-1 inhibitor-associated colitis. On the basis of our data, we successfully targeted this pathway in a patient with refractory ICI-colitis, using the JAK inhibitor tofacitinib. CONCLUSIONS: Interferon gamma-producing CD8+ TRM cells are a pathological hallmark of ICI-colitis and a novel target for therapy.

EPISPOT: An epigenome-driven approach for detecting and interpreting hotspots in molecular QTL studies.

We present EPISPOT, a fully joint framework which exploits large panels of epigenetic annotations as variant-level information to enhance molecular quantitative trait locus (QTL) mapping. Thanks to a purpose-built Bayesian inferential algorithm, EPISPOT accommodates functional information for both cis and trans actions, including QTL hotspot effects. It effectively couples simultaneous QTL analysis of thousands of genetic variants and molecular traits with hypothesis-free selection of biologically interpretable annotations which directly contribute to the QTL effects. This unified, epigenome-aided learning boosts statistical power and sheds light on the regulatory basis of the uncovered hits; EPISPOT therefore marks an essential step toward improving the challenging detection and functional interpretation of trans-acting genetic variants and hotspots. We illustrate the advantages of EPISPOT in simulations emulating real-data conditions and in a monocyte expression QTL study, which confirms known hotspots and finds other signals, as well as plausible mechanisms of action. In particular, by highlighting the role of monocyte DNase-I sensitivity sites from >150 epigenetic annotations, we clarify the mediation effects and cell-type specificity of major hotspots close to the lysozyme gene. Our approach forgoes the daunting and underpowered task of one-annotation-at-a-time enrichment analyses for prioritizing cis and trans QTL hits and is tailored to any transcriptomic, proteomic, or metabolomic QTL problem. By enabling principled epigenome-driven QTL mapping transcriptome-wide, EPISPOT helps progress toward a better functional understanding of genetic regulation.

Checkpoint-blocker-induced autoimmunity is associated with favourable outcome in metastatic melanoma and distinct T-cell expression profiles.

BACKGROUND: Immune checkpoint blockers (ICBs) activate CD8+ T cells, eliciting both anti-cancer activity and immune-related adverse events (irAEs). The relationship of irAEs with baseline parameters and clinical outcome is unclear. METHODS: Retrospective evaluation of irAEs on survival was performed across primary (N = 144) and secondary (N = 211) independent cohorts of patients with metastatic melanoma receiving single agent (pembrolizumab/nivolumab-sICB) or combination (nivolumab and ipilimumab-cICB) checkpoint blockade. RNA from pre-treatment and post-treatment CD8+ T cells was sequenced and differential gene expression according to irAE development assessed. RESULTS: 58.3% of patients developed early irAEs and this was associated with longer progression-free (PFS) and overall survival (OS) across both cohorts (log-rank test, OS: P < 0.0001). Median survival for patients without irAEs was 16.6 months (95% CI: 10.9-33.4) versus not-reached (P = 2.8 × 10-6). Pre-treatment monocyte and neutrophil counts, but not BMI, were additional predictors of clinical outcome. Differential expression of numerous gene pathway members was observed in CD8+ T cells according to irAE development, and patients not developing irAEs demonstrating upregulated CXCR1 pre- and post-treatment. CONCLUSIONS: Early irAE development post-ICB is associated with favourable survival in MM. Development of irAEs is coupled to expression of numerous gene pathways, suggesting irAE development in-part reflects baseline immune activation.

Peripheral CD8+ T cell characteristics associated with durable responses to immune checkpoint blockade in patients with metastatic melanoma.

Immune checkpoint blockade (ICB) of PD-1 and CTLA-4 to treat metastatic melanoma (MM) has variable therapeutic benefit. To explore this in peripheral samples, we characterized CD8+ T cell gene expression across a cohort of patients with MM receiving anti-PD-1 alone (sICB) or in combination with anti-CTLA-4 (cICB). Whereas CD8+ transcriptional responses to sICB and cICB involve a shared gene set, the magnitude of cICB response is over fourfold greater, with preferential induction of mitosis- and interferon-related genes. Early samples from patients with durable clinical benefit demonstrated overexpression of T cell receptor-encoding genes. By mapping T cell receptor clonality, we find that responding patients have more large clones (those occupying >0.5% of repertoire) post-treatment than non-responding patients or controls, and this correlates with effector memory T cell percentage. Single-cell RNA-sequencing of eight post-treatment samples demonstrates that large clones overexpress genes implicated in cytotoxicity and characteristic of effector memory T cells, including CCL4, GNLY and NKG7. The 6-month clinical response to ICB in patients with MM is associated with the large CD8+ T cell clone count 21 d after treatment and agnostic to clonal specificity, suggesting that post-ICB peripheral CD8+ clonality can provide information regarding long-term treatment response and, potentially, facilitate treatment stratification.

Context-specific regulation of surface and soluble IL7R expression by an autoimmune risk allele.

IL-7 is a key factor in T cell immunity and common variants at IL7R, encoding its receptor, are associated with autoimmune disease susceptibility. IL7R mRNA is induced in stimulated monocytes, yet a function for IL7R in monocyte biology remains unexplored. Here we characterize genetic regulation of IL7R at the protein level in healthy individuals, and find that monocyte surface and soluble IL7R (sIL7R) are markedly induced by lipopolysaccharide. In monocytes, both surface IL7R and sIL7R expression strongly associate with allelic carriage of rs6897932, a disease-associated IL7R polymorphism. Monocytes produce more sIL7R than CD4 + T cells, and the amount is additionally correlated with the expression of DDX39A, encoding a splicing factor. Synovial fluid-derived monocytes from patients with spondyloarthritis are enriched for IL7R+ cells with a unique transcriptional profile that overlaps with IL-7-induced gene sets. Our data thus suggest a previously unappreciated function for monocytes in IL-7 biology and IL7R-associated diseases.

Regulatory Crosstalk of Doxorubicin, Estradiol and TNFα Combined Treatment in Breast Cancer-derived Cell Lines.

We present a new model of ESR1 network regulation based on analysis of Doxorubicin, Estradiol, and TNFα combination treatment in MCF-7. We used Doxorubicin as a therapeutic agent, TNFα as marker and mediator of an inflammatory microenvironment and 17ß-Estradiol (E2) as an agonist of Estrogen Receptors, known predisposing factor for hormone-driven breast cancer, whose pharmacological inhibition reduces the risk of breast cancer recurrence. Based on the results of transcriptomics analysis, we found 71 differentially expressed genes that are specific for the combination treatment with Doxorubicin + Estradiol + TNFα in comparison with single or double treatments. The responsiveness to the triple treatment was examined for seven genes by qPCR, of which six were validated, and then extended to four additional cell lines differing for p53 and/or ER status. The results of differential regulation enrichment analysis highlight the role of the ESR1 network that included 36 of 71 specific differentially expressed genes. We propose that the combined activation of p53 and NF-kB transcription factors significantly influences ligand-dependent, ER-driven transcriptional responses, also of the ESR1 gene itself. These results provide a model of coordinated interaction of TFs to explain the Doxorubicin, E2 and TNFα induced repression mechanisms.

Systematic exploration of cell morphological phenotypes associated with a transcriptomic query.

Cell morphological phenotypes, including shape, size, intensity, and texture of cellular compartments have been shown to change in response to perturbation with small molecule compounds. Image-based cell profiling or cell morphological profiling has been used to associate changes of cell morphological features with alterations in cellular function and to infer molecular mechanisms of action. Recently, the Library of Integrated Network-based Cellular Signatures (LINCS) Project has measured gene expression and performed image-based cell profiling on cell lines treated with 9515 unique compounds. These data provide an opportunity to study the interdependence between transcription and cell morphology. Previous methods to investigate cell phenotypes have focused on targeting candidate genes as components of known pathways, RNAi morphological profiling, and cataloging morphological defects; however, these methods do not provide an explicit model to link transcriptomic changes with corresponding alterations in morphology. To address this, we propose a cell morphology enrichment analysis to assess the association between transcriptomic alterations and changes in cell morphology. Additionally, for a new transcriptomic query, our approach can be used to predict associated changes in cellular morphology. We demonstrate the utility of our method by applying it to cell morphological changes in a human bone osteosarcoma cell line.

Systems view of adipogenesis via novel omics-driven and tissue-specific activity scoring of network functional modules.

The investigation of the complex processes involved in cellular differentiation must be based on unbiased, high throughput data processing methods to identify relevant biological pathways. A number of bioinformatics tools are available that can generate lists of pathways ranked by statistical significance (i.e. by p-value), while ideally it would be desirable to functionally score the pathways relative to each other or to other interacting parts of the system or process. We describe a new computational method (Network Activity Score Finder - NASFinder) to identify tissue-specific, omics-determined sub-networks and the connections with their upstream regulator receptors to obtain a systems view of the differentiation of human adipocytes. Adipogenesis of human SBGS pre-adipocyte cells in vitro was monitored with a transcriptomic data set comprising six time points (0, 6, 48, 96, 192, 384 hours). To elucidate the mechanisms of adipogenesis, NASFinder was used to perform time-point analysis by comparing each time point against the control (0 h) and time-lapse analysis by comparing each time point with the previous one. NASFinder identified the coordinated activity of seemingly unrelated processes between each comparison, providing the first systems view of adipogenesis in culture. NASFinder has been implemented into a web-based, freely available resource associated with novel, easy to read visualization of omics data sets and network modules.

Discovering dominant pathways and signal-response relationships in signaling networks through nonparametric approaches.

A signaling pathway is a sequence of proteins and passenger molecules that transmits information from the cell surface to target molecules. Understanding signal transduction process requires detailed description of the involved pathways. Several methods and tools resolved this problem by incorporating genomic and proteomic data. However, the difficulty of obtaining prior knowledge of complex signaling networks limited the applicability of these tools. In this study, based on the simulation of signal flow in signaling network, we introduce a method for determining dominant pathways and signal response to stimulations. The model uses topology-weighted transit compartment approach and comprises four main steps which include weighting the edges, simulating signal transduction in the network (weighting the nodes), finding paths between initial and target nodes, and assigning a significance score to each path. We applied the proposed model to eighty-three signaling networks by using biologically derived source and sink molecules. The recovered dominant paths matched many known signaling pathways and suggesting a promising index to analyze the phenotype essentiality of molecule encoding paths. We also modeled the stimulus-response relations in long and short-term synaptic plasticity based on the dominant signaling pathway concept. We showed that the proposed method not only accurately determines dominant signaling pathways, but also identifies effective points of intervention in signal transduction.

Nonparametric simulation of signal transduction networks with semi-synchronized update.

Simulating signal transduction in cellular signaling networks provides predictions of network dynamics by quantifying the changes in concentration and activity-level of the individual proteins. Since numerical values of kinetic parameters might be difficult to obtain, it is imperative to develop non-parametric approaches that combine the connectivity of a network with the response of individual proteins to signals which travel through the network. The activity levels of signaling proteins computed through existing non-parametric modeling tools do not show significant correlations with the observed values in experimental results. In this work we developed a non-parametric computational framework to describe the profile of the evolving process and the time course of the proportion of active form of molecules in the signal transduction networks. The model is also capable of incorporating perturbations. The model was validated on four signaling networks showing that it can effectively uncover the activity levels and trends of response during signal transduction process.