Comparative analysis of differential gene expression analysis tools for single-cell RNA sequencing data

Transcriptomic Evaluation of a Stress Vulnerability Network Using Single-Cell RNA Sequencing in Mouse Prefrontal Cortex

BACKGROUND

Increased vulnerability to stress is a major risk factor for several mood disorders, including major depressive disorder. Although cellular and molecular mechanisms associated with depressive behaviors following stress have been identified, little is known about the mechanisms that confer the vulnerability that predisposes individuals to future damage from chronic stress.

METHODS

We used multisite in vivo neurophysiology in freely behaving male and female C57BL/6 mice (n = 12) to measure electrical brain network activity previously identified as indicating a latent stress vulnerability brain state. We combined this neurophysiological approach with single-cell RNA sequencing of the prefrontal cortex to identify distinct transcriptomic differences between groups of mice with inherent high and low stress vulnerability.

RESULTS

We identified hundreds of differentially expressed genes (padjusted < .05) across 5 major cell types in animals with high and low stress vulnerability brain network activity. This unique analysis revealed that GABAergic (gamma-aminobutyric acidergic) neuron gene expression contributed most to the network activity of the stress vulnerability brain state. Upregulation of mitochondrial and metabolic pathways also distinguished high and low vulnerability brain states, especially in inhibitory neurons. Importantly, genes that were differentially regulated with vulnerability network activity significantly overlapped (above chance) with those identified by genome-wide association studies as having single nucleotide polymorphisms significantly associated with depression as well as genes more highly expressed in postmortem prefrontal cortex of patients with major depressive disorder.

CONCLUSIONS

This is the first study to identify cell types and genes involved in a latent stress vulnerability state in the brain.

Theoretical framework for the difference of two negative binomial distributions and its application in comparative analysis of sequencing data

High-throughput sequencing (HTS) technologies have been instrumental in investigating biological questions at the bulk and single-cell levels. Comparative analysis of two HTS data sets often relies on testing the statistical significance for the difference of two negative binomial distributions (DOTNB). Although negative binomial distributions are well studied, the theoretical results for DOTNB remain largely unexplored. Here, we derive basic analytical results for DOTNB and examine its asymptotic properties. As a state-of-the-art application of DOTNB, we introduce DEGage, a computational method for detecting differentially expressed genes (DEGs) in scRNA-seq data. DEGage calculates the mean of the sample-wise differences of gene expression levels as the test statistic and determines significant differential expression by computing the P-value with DOTNB. Extensive validation using simulated and real scRNA-seq data sets demonstrates that DEGage outperforms five popular DEG analysis tools: DEGseq2, DEsingle, edgeR, Monocle3, and scDD. DEGage is robust against high dropout levels and exhibits superior sensitivity when applied to balanced and imbalanced data sets, even with small sample sizes. We utilize DEGage to analyze prostate cancer scRNA-seq data sets and identify marker genes for 17 cell types. Furthermore, we apply DEGage to scRNA-seq data sets of mouse neurons with and without fear memory and reveal eight potential memory-related genes overlooked in previous analyses. The theoretical results and supporting software for DOTNB can be widely applied to comparative analyses of dispersed count data in HTS and broad research questions.

Diagnostic challenges in complicated case of glioblastoma

Glioblastoma is the commonest primary malignant brain tumor, with a very poor prognosis and short overall survival. It is characterized by its high intra- and intertumoral heterogeneity, in terms of both the level of single-nucleotide variants, copy number alterations, and aneuploidy. Therefore, routine diagnosis can be challenging in some cases. We present a complicated case of glioblastoma, which was characterized with five cytogenomic methods: interphase fluorescence in situ hybridization, multiplex ligation-dependent probe amplification, comparative genomic hybridization array and single-nucleotide polymorphism, targeted gene panel, and whole-genome sequencing. These cytogenomic methods revealed classical findings associated with glioblastoma, such as a lack of IDH and TERT mutations, gain of chromosome 7, and loss of chromosome 10. At least three pathological clones were identified, including one with whole-genome duplication, and one with loss of 1p and suspected loss of 19q. Deletion and mutation of the TP53 gene were detected with numerous breakends on 17p and 20q. Based on these findings, we recommend a combined approach to the diagnosis of glioblastoma involving the detection of copy number alterations, mutations, and aneuploidy. The choice of the best combination of methods is based on cost, time required, staff expertise, and laboratory equipment. This integrated strategy could contribute directly to tangible improvements in the diagnosis, prognosis, and prediction of the therapeutic responses of patients with brain tumors.

IMMUNE AND MOLECULAR CORRELATES OF RESPONSE TO IMMUNOTHERAPY REVEALED BY BRAIN-METASTATIC MELANOMA MODELS

Despite the promising results of immune checkpoint blockade (ICB) therapy, outcomes for patients with brain metastasis (BrM) remain poor. Identifying resistance mechanisms has been hindered by limited access to patient samples and relevant preclinical models. Here, we developed two mouse melanoma BrM models that recapitulate the disparate responses to ICB seen in patients. We demonstrate that these models capture the cellular and molecular complexity of human disease and reveal key factors shaping the tumor microenvironment and influencing ICB response. BR1-responsive tumor cells express inflammatory programs that polarize microglia into reactive states, eliciting robust T cell recruitment. In contrast, BR3-resistant melanoma cells are enriched in neurological programs and exploit tolerance mechanisms to maintain microglia homeostasis and limit T cell infiltration. In humans, BR1 and BR3 expression signatures correlate positively or negatively with T cell infiltration and BrM patient outcomes, respectively. Our study provides clinically relevant models and uncovers mechanistic insights into BrM ICB responses, offering potential biomarkers and therapeutic targets to improve therapy efficacy.

cytoKernel: Robust kernel embeddings for assessing differential expression of single cell data

High-throughput sequencing of single-cell data can be used to rigorously evlauate cell specification and enable intricate variations between groups or conditions. Many popular existing methods for differential expression target differences in aggregate measurements (mean, median, sum) and limit their approaches to detect only global differential changes. We present a robust method for differential expression of single-cell data using a kernel-based score test, cytoKernel. cytoKernel is specifically designed to assess the differential expression of single cell RNA sequencing and high-dimensional flow or mass cytometry data using the full probability distribution pattern. cytoKernel is based on kernel embeddings which employs the probability distributions of the single cell data, by calculating the pairwise divergence/distance between distributions of subjects. It can detect both patterns involving aggregate changes, as well as more elusive variations that are often overlooked due to the multimodal characteristics of single cell data. We performed extensive benchmarks across both simulated and real data sets from mass cytometry data and single-cell RNA sequencing. The cytoKernel procedure effectively controls the False Discovery Rate (FDR) and shows favourable performance compared to existing methods. The method is able to identify more differential patterns than existing approaches. We apply cytoKernel to assess gene expression and protein marker expression differences from cell subpopulations in various publicly available single-cell RNAseq and mass cytometry data sets. The methods described in this paper are implemented in the open-source R package cytoKernel, which is freely available from Bioconductor at http://bioconductor.org/packages/cytoKernel.

Single-cell transcriptome atlases of soybean root and mature nodule reveal new regulatory programs that control the nodulation process

The soybean root system is complex. In addition to being composed of various cell types, the soybean root system includes the primary root, the lateral roots, and the nodule, an organ in which mutualistic symbiosis with N-fixing rhizobia occurs. A mature soybean root nodule is characterized by a central infection zone where atmospheric nitrogen is fixed and assimilated by the symbiont, resulting from the close cooperation between the plant cell and the bacteria. To date, the transcriptome of individual cells isolated from developing soybean nodules has been established, but the transcriptomic signatures of cells from the mature soybean nodule have not yet been characterized. Using single-nucleus RNA-seq and Molecular Cartography technologies, we precisely characterized the transcriptomic signature of soybean root and mature nodule cell types and revealed the co-existence of different sub-populations of B. diazoefficiens-infected cells in the mature soybean nodule, including those actively involved in nitrogen fixation and those engaged in senescence. Mining of the single-cell-resolution nodule transcriptome atlas and the associated gene co-expression network confirmed the role of known nodulation-related genes and identified new genes that control the nodulation process. For instance, we functionally characterized the role of GmFWL3, a plasma membrane microdomain-associated protein that controls rhizobial infection. Our study reveals the unique cellular complexity of the mature soybean nodule and helps redefine the concept of cell types when considering the infection zone of the soybean nodule.

RankCompV3: a differential expression analysis algorithm based on relative expression orderings and applications in single-cell RNA transcriptomics

BACKGROUND

Effective identification of differentially expressed genes (DEGs) has been challenging for single-cell RNA sequencing (scRNA-seq) profiles. Many existing algorithms have high false positive rates (FPRs) and often fail to identify weak biological signals.

RESULTS

We present a novel method for identifying DEGs in scRNA-seq data called RankCompV3. It is based on the comparison of relative expression orderings (REOs) of gene pairs which are determined by comparing the expression levels of a pair of genes in a set of single-cell profiles. The numbers of genes with consistently higher or lower expression levels than the gene of interest are counted in two groups in comparison, respectively, and the result is tabulated in a 3 × 3 contingency table which is tested by McCullagh's method to determine if the gene is dysregulated. In both simulated and real scRNA-seq data, RankCompV3 tightly controlled the FPR and demonstrated high accuracy, outperforming 11 other common single-cell DEG detection algorithms. Analysis with either regular single-cell or synthetic pseudo-bulk profiles produced highly concordant DEGs with the ground-truth. In addition, RankCompV3 demonstrates higher sensitivity to weak biological signals than other methods. The algorithm was implemented using Julia and can be called in R. The source code is available at https://github.com/pathint/RankCompV3.jl .

CONCLUSIONS

The REOs-based algorithm is a valuable tool for analyzing single-cell RNA profiles and identifying DEGs with high accuracy and sensitivity.

Revisiting Fold-Change Calculation: Preference for Median or Geometric Mean over Arithmetic Mean-Based Methods

Background: Fold change is a common metric in biomedical research for quantifying group differences in omics variables. However, inconsistent calculation methods and inadequate reporting lead to discrepancies in results. This study evaluated various fold-change calculation methods aiming at a recommendation of a preferred approach. Methods: The primary distinction in fold-change calculations lies in defining group expected values for log ratio computation. To challenge method interchangeability in a "stress test" scenario, we generated diverse artificial data sets with varying distributions (identity, uniform, normal, log-normal, and a mixture of these) and compared calculated fold-changes to known values. Additionally, we analyzed a multi-omics biomedical data set to estimate to what extent the findings apply to real-world data. Results: Using arithmetic means as expected values for treatment and reference groups yielded inaccurate fold-change values more frequently than other methods, particularly when subgroup distributions and/or standard deviations differed significantly. Conclusions: The arithmetic mean method, often perceived as standard or picked without considering alternatives, is inferior to other definitions of the group expected value. Methods using median, geometric mean, or paired fold-change combinations are more robust against violations of equal variances or dissimilar group distributions. Adhering to methods less sensitive to data distribution without trade-offs and accurately reporting calculation methods in scientific reports is a reasonable practice to ensure correct interpretation and reproducibility.

Leveraging neighborhood representations of single-cell data to achieve sensitive DE testing with miloDE

Single-cell RNA-sequencing enables testing for differential expression (DE) between conditions at a cell type level. While powerful, one of the limitations of such approaches is that the sensitivity of DE testing is dictated by the sensitivity of clustering, which is often suboptimal. To overcome this, we present miloDE-a cluster-free framework for DE testing (available as an open-source R package). We illustrate the performance of miloDE on both simulated and real data. Using miloDE, we identify a transient hemogenic endothelia-like state in mouse embryos lacking Tal1 and detect distinct programs during macrophage activation in idiopathic pulmonary fibrosis.

Directly selecting cell-type marker genes for single-cell clustering analyses

In single-cell RNA sequencing (scRNA-seq) studies, cell types and their marker genes are often identified by clustering and differentially expressed gene (DEG) analysis. A common practice is to select genes using surrogate criteria such as variance and deviance, then cluster them using selected genes and detect markers by DEG analysis assuming known cell types. The surrogate criteria can miss important genes or select unimportant genes, while DEG analysis has the selection-bias problem. We present Festem, a statistical method for the direct selection of cell-type markers for downstream clustering. Festem distinguishes marker genes with heterogeneous distribution across cells that are cluster informative. Simulation and scRNA-seq applications demonstrate that Festem can sensitively select markers with high precision and enables the identification of cell types often missed by other methods. In a large intrahepatic cholangiocarcinoma dataset, we identify diverse CD8+ T cell types and potential prognostic marker genes.

Traumatic brain injury promotes neurogenesis at the cost of astrogliogenesis in the adult hippocampus of male mice

Traumatic brain injury (TBI) can result in long-lasting changes in hippocampal function. The changes induced by TBI on the hippocampus contribute to cognitive deficits. The adult hippocampus harbors neural stem cells (NSCs) that generate neurons (neurogenesis), and astrocytes (astrogliogenesis). While deregulation of hippocampal NSCs and neurogenesis have been observed after TBI, it is not known how TBI may affect hippocampal astrogliogenesis. Using a controlled cortical impact model of TBI in male mice, single cell RNA sequencing and spatial transcriptomics, we assessed how TBI affected hippocampal NSCs and the neuronal and astroglial lineages derived from them. We observe an increase in NSC-derived neuronal cells and a concomitant decrease in NSC-derived astrocytic cells, together with changes in gene expression and cell dysplasia within the dentate gyrus. Here, we show that TBI modifies NSC fate to promote neurogenesis at the cost of astrogliogenesis and identify specific cell populations as possible targets to counteract TBI-induced cellular changes in the adult hippocampus.

Decoding frontotemporal and cell-type-specific vulnerabilities to neuropsychiatric disorders and psychoactive drugs

Frontotemporal lobe abnormalities are linked to neuropsychiatric disorders and cognition, but the role of cellular heterogeneity between temporal lobe (TL) and frontal lobe (FL) in the vulnerability to genetic risk factors remains to be elucidated. We integrated single-nucleus transcriptome analysis in 'fresh' human FL and TL with genetic susceptibility, gene dysregulation in neuropsychiatric disease and psychoactive drug response data. We show how intrinsic differences between TL and FL contribute to the vulnerability of specific cell types to both genetic risk factors and psychoactive drugs. Neuronal populations, specifically PVALB neurons, were most highly vulnerable to genetic risk factors for psychiatric disease. These psychiatric disease-associated genes were mostly upregulated in the TL, and dysregulated in the brain of patients with obsessive-compulsive disorder, bipolar disorder and schizophrenia. Among these genes, GRIN2A and SLC12A5, implicated in schizophrenia and bipolar disorder, were significantly upregulated in TL PVALB neurons and in psychiatric disease patients' brain. PVALB neurons from the TL were twofold more vulnerable to psychoactive drugs than to genetic risk factors, showing the influence and specificity of frontotemporal lobe differences on cell vulnerabilities. These studies provide a cell type resolved map of the impact of brain regional differences on cell type vulnerabilities in neuropsychiatric disorders.

Modeling type 1 diabetes progression using machine learning and single-cell transcriptomic measurements in human islets

Type 1 diabetes (T1D) is a chronic condition in which beta cells are destroyed by immune cells. Despite progress in immunotherapies that could delay T1D onset, early detection of autoimmunity remains challenging. Here, we evaluate the utility of machine learning for early prediction of T1D using single-cell analysis of islets. Using gradient-boosting algorithms, we model changes in gene expression of single cells from pancreatic tissues in T1D and non-diabetic organ donors. We assess if mathematical modeling could predict the likelihood of T1D development in non-diabetic autoantibody-positive donors. While most autoantibody-positive donors are predicted to be non-diabetic, select donors with unique gene signatures are classified as T1D. Our strategy also reveals a shared gene signature in distinct T1D-associated models across cell types, suggesting a common effect of the disease on transcriptional outputs of these cells. Our study establishes a precedent for using machine learning in early detection of T1D.

Identification and Pharmacological Targeting of Treatment-Resistant, Stem-like Breast Cancer Cells for Combination Therapy

Tumors frequently harbor isogenic yet epigenetically distinct subpopulations of multi-potent cells with high tumor-initiating potential-often called Cancer Stem-Like Cells (CSLCs). These can display preferential resistance to standard-of-care chemotherapy. Single-cell analyses can help elucidate Master Regulator (MR) proteins responsible for governing the transcriptional state of these cells, thus revealing complementary dependencies that may be leveraged via combination therapy. Interrogation of single-cell RNA sequencing profiles from seven metastatic breast cancer patients, using perturbational profiles of clinically relevant drugs, identified drugs predicted to invert the activity of MR proteins governing the transcriptional state of chemoresistant CSLCs, which were then validated by CROP-seq assays. The top drug, the anthelmintic albendazole, depleted this subpopulation in vivo without noticeable cytotoxicity. Moreover, sequential cycles of albendazole and paclitaxel-a commonly used chemotherapeutic -displayed significant synergy in a patient-derived xenograft (PDX) from a TNBC patient, suggesting that network-based approaches can help develop mechanism-based combinatorial therapies targeting complementary subpopulations.

Data normalization for addressing the challenges in the analysis of single-cell transcriptomic datasets

BACKGROUND

Normalization is a critical step in the analysis of single-cell RNA-sequencing (scRNA-seq) datasets. Its main goal is to make gene counts comparable within and between cells. To do so, normalization methods must account for technical and biological variability. Numerous normalization methods have been developed addressing different sources of dispersion and making specific assumptions about the count data.

MAIN BODY

The selection of a normalization method has a direct impact on downstream analysis, for example differential gene expression and cluster identification. Thus, the objective of this review is to guide the reader in making an informed decision on the most appropriate normalization method to use. To this aim, we first give an overview of the different single cell sequencing platforms and methods commonly used including isolation and library preparation protocols. Next, we discuss the inherent sources of variability of scRNA-seq datasets. We describe the categories of normalization methods and include examples of each. We also delineate imputation and batch-effect correction methods. Furthermore, we describe data-driven metrics commonly used to evaluate the performance of normalization methods. We also discuss common scRNA-seq methods and toolkits used for integrated data analysis.

CONCLUSIONS

According to the correction performed, normalization methods can be broadly classified as within and between-sample algorithms. Moreover, with respect to the mathematical model used, normalization methods can further be classified into: global scaling methods, generalized linear models, mixed methods, and machine learning-based methods. Each of these methods depict pros and cons and make different statistical assumptions. However, there is no better performing normalization method. Instead, metrics such as silhouette width, K-nearest neighbor batch-effect test, or Highly Variable Genes are recommended to assess the performance of normalization methods.

Recent advances in differential expression analysis for single-cell RNA-seq and spatially resolved transcriptomic studies

Differential expression (DE) analysis is a necessary step in the analysis of single-cell RNA sequencing (scRNA-seq) and spatially resolved transcriptomics (SRT) data. Unlike traditional bulk RNA-seq, DE analysis for scRNA-seq or SRT data has unique characteristics that may contribute to the difficulty of detecting DE genes. However, the plethora of DE tools that work with various assumptions makes it difficult to choose an appropriate one. Furthermore, a comprehensive review on detecting DE genes for scRNA-seq data or SRT data from multi-condition, multi-sample experimental designs is lacking. To bridge such a gap, here, we first focus on the challenges of DE detection, then highlight potential opportunities that facilitate further progress in scRNA-seq or SRT analysis, and finally provide insights and guidance in selecting appropriate DE tools or developing new computational DE methods.

Machine learning on alignment features for parent-of-origin classification of simulated hybrid RNA-seq

BACKGROUND

Parent-of-origin allele-specific gene expression (ASE) can be detected in interspecies hybrids by virtue of RNA sequence variants between the parental haplotypes. ASE is detectable by differential expression analysis (DEA) applied to the counts of RNA-seq read pairs aligned to parental references, but aligners do not always choose the correct parental reference.

RESULTS

We used public data for species that are known to hybridize. We measured our ability to assign RNA-seq read pairs to their proper transcriptome or genome references. We tested software packages that assign each read pair to a reference position and found that they often favored the incorrect species reference. To address this problem, we introduce a post process that extracts alignment features and trains a random forest classifier to choose the better alignment. On each simulated hybrid dataset tested, our machine-learning post-processor achieved higher accuracy than the aligner by itself at choosing the correct parent-of-origin per RNA-seq read pair.

CONCLUSIONS

For the parent-of-origin classification of RNA-seq, machine learning can improve the accuracy of alignment-based methods. This approach could be useful for enhancing ASE detection in interspecies hybrids, though RNA-seq from real hybrids may present challenges not captured by our simulations. We believe this is the first application of machine learning to this problem domain.

Drug target prediction through deep learning functional representation of gene signatures

Many machine learning applications in bioinformatics currently rely on matching gene identities when analyzing input gene signatures and fail to take advantage of preexisting knowledge about gene functions. To further enable comparative analysis of OMICS datasets, including target deconvolution and mechanism of action studies, we develop an approach that represents gene signatures projected onto their biological functions, instead of their identities, similar to how the word2vec technique works in natural language processing. We develop the Functional Representation of Gene Signatures (FRoGS) approach by training a deep learning model and demonstrate that its application to the Broad Institute's L1000 datasets results in more effective compound-target predictions than models based on gene identities alone. By integrating additional pharmacological activity data sources, FRoGS significantly increases the number of high-quality compound-target predictions relative to existing approaches, many of which are supported by in silico and/or experimental evidence. These results underscore the general utility of FRoGS in machine learning-based bioinformatics applications. Prediction networks pre-equipped with the knowledge of gene functions may help uncover new relationships among gene signatures acquired by large-scale OMICs studies on compounds, cell types, disease models, and patient cohorts.

A comparison of marker gene selection methods for single-cell RNA sequencing data

BACKGROUND

The development of single-cell RNA sequencing (scRNA-seq) has enabled scientists to catalog and probe the transcriptional heterogeneity of individual cells in unprecedented detail. A common step in the analysis of scRNA-seq data is the selection of so-called marker genes, most commonly to enable annotation of the biological cell types present in the sample. In this paper, we benchmark 59 computational methods for selecting marker genes in scRNA-seq data.

RESULTS

We compare the performance of the methods using 14 real scRNA-seq datasets and over 170 additional simulated datasets. Methods are compared on their ability to recover simulated and expert-annotated marker genes, the predictive performance and characteristics of the gene sets they select, their memory usage and speed, and their implementation quality. In addition, various case studies are used to scrutinize the most commonly used methods, highlighting issues and inconsistencies.

CONCLUSIONS

Overall, we present a comprehensive evaluation of methods for selecting marker genes in scRNA-seq data. Our results highlight the efficacy of simple methods, especially the Wilcoxon rank-sum test, Student's t-test, and logistic regression.

α1,3-fucosylation treatment improves cord blood CD34 negative hematopoietic stem cell navigation

For almost two decades, clinicians have overlooked the diagnostic potential of CD34neg hematopoietic stem cells because of their limited homing capacity relative to CD34posHSCs when injected intravenously. This has contributed to the lack of appeal of using umbilical cord blood in HSC transplantation because its stem cell count is lower than bone marrow. The present study reveals that the homing and engraftment of CD34negHSCs can be improved by adding the Sialyl Lewis X molecule via α1,3-fucosylation. This unlocks the potential for using this more primitive stem cell to treat blood disorders because our findings show CD34negHSCs have the capacity to regenerate cells in the bone marrow of mice for several months. Furthermore, our RNA sequencing analysis revealed that CD34negHSCs have unique adhesion pathways, downregulated in CD34posHSCs, that facilitate interaction with the bone marrow niche. Our findings suggest that CD34neg cells will best thrive when the HSC resides in its microenvironment.

scKWARN: Kernel-weighted-average robust normalization for single-cell RNA-seq data

MOTIVATION

Single-cell RNA-seq normalization is an essential step to correct unwanted biases caused by sequencing depth, capture efficiency, dropout, and other technical factors. Existing normalization methods primarily reduce biases arising from sequencing depth by modeling count-depth relationship and/or assuming a specific distribution for read counts. However, these methods may lead to over or under-correction due to presence of technical biases beyond sequencing depth and the restrictive assumption on models and distributions.

RESULTS

We present scKWARN, a Kernel Weighted Average Robust Normalization designed to correct known or hidden technical confounders without assuming specific data distributions or count-depth relationships. scKWARN generates a pseudo expression profile for each cell by borrowing information from its fuzzy technical neighbors through a kernel smoother. It then compares this profile against the reference derived from cells with the same bimodality patterns to determine the normalization factor. As demonstrated in both simulated and real datasets, scKWARN outperforms existing methods in removing a variety of technical biases while preserving true biological heterogeneity.

AVAILABILITY AND IMPLEMENTATION

scKWARN is freely available at https://github.com/cyhsuTN/scKWARN.

Exploring miRNA-target gene pair detection in disease with coRmiT

A wide range of approaches can be used to detect micro RNA (miRNA)-target gene pairs (mTPs) from expression data, differing in the ways the gene and miRNA expression profiles are calculated, combined and correlated. However, there is no clear consensus on which is the best approach across all datasets. Here, we have implemented multiple strategies and applied them to three distinct rare disease datasets that comprise smallRNA-Seq and RNA-Seq data obtained from the same samples, obtaining mTPs related to the disease pathology. All datasets were preprocessed using a standardized, freely available computational workflow, DEG_workflow. This workflow includes coRmiT, a method to compare multiple strategies for mTP detection. We used it to investigate the overlap of the detected mTPs with predicted and validated mTPs from 11 different databases. Results show that there is no clear best strategy for mTP detection applicable to all situations. We therefore propose the integration of the results of the different strategies by selecting the one with the highest odds ratio for each miRNA, as the optimal way to integrate the results. We applied this selection-integration method to the datasets and showed it to be robust to changes in the predicted and validated mTP databases. Our findings have important implications for miRNA analysis. coRmiT is implemented as part of the ExpHunterSuite Bioconductor package available from https://bioconductor.org/packages/ExpHunterSuite.

Mapping the transcriptome: Realizing the full potential of spatial data analysis

RNA sequencing in situ allows for whole-transcriptome characterization at high resolution, while retaining spatial information. These data present an analytical challenge for bioinformatics-how to leverage spatial information effectively? Properties of data with a spatial dimension require special handling, which necessitate a different set of statistical and inferential considerations when compared to non-spatial data. The geographical sciences primarily use spatial data and have developed methods to analye them. Here we discuss the challenges associated with spatial analysis and examine how we can take advantage of practice from the geographical sciences to realize the full potential of spatial information in transcriptomic datasets.

FEZ1 participates in human embryonic brain development by modulating neuronal progenitor subpopulation specification and migrations

Mutations in the human fasciculation and elongation protein zeta 1 (FEZ1) gene are found in schizophrenia and Jacobsen syndrome patients. Here, using human cerebral organoids (hCOs), we show that FEZ1 expression is turned on early during brain development and is detectable in both neuroprogenitor subtypes and immature neurons. FEZ1 deletion disrupts expression of neuronal and synaptic development genes. Using single-cell RNA sequencing, we detected abnormal expansion of homeodomain-only protein homeobox (HOPX)- outer radial glia (oRG), concurrent with a reduction of HOPX+ oRG, in FEZ1-null hCOs. HOPX- oRGs show higher cell mobility as compared to HOPX+ oRGs. Ectopic localization of neuroprogenitors to the outer layer is seen in FEZ1-null hCOs. Anomalous encroachment of TBR2+ intermediate progenitors into CTIP2+ deep layer neurons further indicated abnormalities in cortical layer formation these hCOs. Collectively, our findings highlight the involvement of FEZ1 in early cortical brain development and how it contributes to neurodevelopmental disorders.

FastCAR: fast correction for ambient RNA to facilitate differential gene expression analysis in single-cell RNA-sequencing datasets

Cell type-specific differential gene expression analyses based on single-cell transcriptome datasets are sensitive to the presence of cell-free mRNA in the droplets containing single cells. This so-called ambient RNA contamination may differ between samples obtained from patients and healthy controls. Current ambient RNA correction methods were not developed specifically for single-cell differential gene expression (sc-DGE) analyses and might therefore not sufficiently correct for ambient RNA-derived signals. Here, we show that ambient RNA levels are highly sample-specific. We found that without ambient RNA correction, sc-DGE analyses erroneously identify transcripts originating from ambient RNA as cell type-specific disease-associated genes. We therefore developed a computationally lean and intuitive correction method, Fast Correction for Ambient RNA (FastCAR), optimized for sc-DGE analysis of scRNA-Seq datasets generated by droplet-based methods including the 10XGenomics Chromium platform. FastCAR uses the profile of transcripts observed in libraries that likely represent empty droplets to determine the level of ambient RNA in each individual sample, and then corrects for these ambient RNA gene expression values. FastCAR can be applied as part of the data pre-processing and QC in sc-DGE workflows comparing scRNA-Seq data in a health versus disease experimental design. We compared FastCAR with two methods previously developed to remove ambient RNA, SoupX and CellBender. All three methods identified additional genes in sc-DGE analyses that were not identified in the absence of ambient RNA correction. However, we show that FastCAR performs better at correcting gene expression values attributed to ambient RNA, resulting in a lower frequency of false-positive observations. Moreover, the use of FastCAR in a sc-DGE workflow increases the cell-type specificity of sc-DGE analyses across disease conditions.

Alternative Transcripts Diversify Genome Function for Phenome Relevance to Health and Diseases

Manipulation using alternative exon splicing (AES), alternative transcription start (ATS), and alternative polyadenylation (APA) sites are key to transcript diversity underlying health and disease. All three are pervasive in organisms, present in at least 50% of human protein-coding genes. In fact, ATS and APA site use has the highest impact on protein identity, with their ability to alter which first and last exons are utilized as well as impacting stability and translation efficiency. These RNA variants have been shown to be highly specific, both in tissue type and stage, with demonstrated importance to cell proliferation, differentiation and the transition from fetal to adult cells. While alternative exon splicing has a limited effect on protein identity, its ubiquity highlights the importance of these minor alterations, which can alter other features such as localization. The three processes are also highly interwoven, with overlapping, complementary, and competing factors, RNA polymerase II and its CTD (C-terminal domain) chief among them. Their role in development means dysregulation leads to a wide variety of disorders and cancers, with some forms of disease disproportionately affected by specific mechanisms (AES, ATS, or APA). Challenges associated with the genome-wide profiling of RNA variants and their potential solutions are also discussed in this review.

GoM DE: interpreting structure in sequence count data with differential expression analysis allowing for grades of membership

Parts-based representations, such as non-negative matrix factorization and topic modeling, have been used to identify structure from single-cell sequencing data sets, in particular structure that is not as well captured by clustering or other dimensionality reduction methods. However, interpreting the individual parts remains a challenge. To address this challenge, we extend methods for differential expression analysis by allowing cells to have partial membership to multiple groups. We call this grade of membership differential expression (GoM DE). We illustrate the benefits of GoM DE for annotating topics identified in several single-cell RNA-seq and ATAC-seq data sets.

Identification of CD133+ intercellsomes in intercellular communication to offset intracellular signal deficit

CD133 (prominin 1) is widely viewed as a cancer stem cell marker in association with drug resistance and cancer recurrence. Herein, we report that with impaired RTK-Shp2-Ras-Erk signaling, heterogenous hepatocytes form clusters that manage to divide during mouse liver regeneration. These hepatocytes are characterized by upregulated CD133 while negative for other progenitor cell markers. Pharmaceutical inhibition of proliferative signaling also induced CD133 expression in various cancer cell types from multiple animal species, suggesting an inherent and common mechanism of stress response. Super-resolution and electron microscopy localize CD133 on intracellular vesicles that apparently migrate between cells, which we name 'intercellsome.' Isolated CD133+ intercellsomes are enriched with mRNAs rather than miRNAs. Single-cell RNA sequencing reveals lower intracellular diversity (entropy) of mitogenic mRNAs in Shp2-deficient cells, which may be remedied by intercellular mRNA exchanges between CD133+ cells. CD133-deficient cells are more sensitive to proliferative signal inhibition in livers and intestinal organoids. These data suggest a mechanism of intercellular communication to compensate for intracellular signal deficit in various cell types.

Single Cell Transcriptome of Stress Vulnerability Network in mouse Prefrontal Cortex

Increased vulnerability to stress is a major risk factor for the manifestation of several mood disorders, including major depressive disorder (MDD). Despite the status of MDD as a significant donor to global disability, the complex integration of genetic and environmental factors that contribute to the behavioral display of such disorders has made a thorough understanding of related etiology elusive. Recent developments suggest that a brain-wide network approach is needed, taking into account the complex interplay of cell types spanning multiple brain regions. Single cell RNA-sequencing technologies can provide transcriptomic profiling at the single-cell level across heterogenous samples. Furthermore, we have previously used local field potential oscillations and machine learning to identify an electrical brain network that is indicative of a predisposed vulnerability state. Thus, this study combined single cell RNA-sequencing (scRNA-Seq) with electrical brain network measures of the stress-vulnerable state, providing a unique opportunity to access the relationship between stress network activity and transcriptomic changes within individual cell types. We found especially high numbers of differentially expressed genes between animals with high and low stress vulnerability brain network activity in astrocytes and glutamatergic neurons but we estimated that vulnerability network activity depends most on GABAergic neurons. High vulnerability network activity included upregulation of microglia and mitochondrial and metabolic pathways, while lower vulnerability involved synaptic regulation. Genes that were differentially regulated with vulnerability network activity significantly overlapped with genes identified as having significant SNPs by human GWAS for depression. Taken together, these data provide the gene expression architecture of a previously uncharacterized stress vulnerability brain state, enabling new understanding and intervention of predisposition to stress susceptibility.

Deciphering the Biology of Circulating Tumor Cells through Single-Cell RNA Sequencing: Implications for Precision Medicine in Cancer

Circulating tumor cells (CTCs) hold unique biological characteristics that directly involve them in hematogenous dissemination. Studying CTCs systematically is technically challenging due to their extreme rarity and heterogeneity and the lack of specific markers to specify metastasis-initiating CTCs. With cutting-edge technology, single-cell RNA sequencing (scRNA-seq) provides insights into the biology of metastatic processes driven by CTCs. Transcriptomics analysis of single CTCs can decipher tumor heterogeneity and phenotypic plasticity for exploring promising novel therapeutic targets. The integrated approach provides a perspective on the mechanisms underlying tumor development and interrogates CTCs interactions with other blood cell types, particularly those of the immune system. This review aims to comprehensively describe the current study on CTC transcriptomic analysis through scRNA-seq technology. We emphasize the workflow for scRNA-seq analysis of CTCs, including enrichment, single cell isolation, and bioinformatic tools applied for this purpose. Furthermore, we elucidated the translational knowledge from the transcriptomic profile of individual CTCs and the biology of cancer metastasis for developing effective therapeutics through targeting key pathways in CTCs.

Best practices for single-cell analysis across modalities

Recent advances in single-cell technologies have enabled high-throughput molecular profiling of cells across modalities and locations. Single-cell transcriptomics data can now be complemented by chromatin accessibility, surface protein expression, adaptive immune receptor repertoire profiling and spatial information. The increasing availability of single-cell data across modalities has motivated the development of novel computational methods to help analysts derive biological insights. As the field grows, it becomes increasingly difficult to navigate the vast landscape of tools and analysis steps. Here, we summarize independent benchmarking studies of unimodal and multimodal single-cell analysis across modalities to suggest comprehensive best-practice workflows for the most common analysis steps. Where independent benchmarks are not available, we review and contrast popular methods. Our article serves as an entry point for novices in the field of single-cell (multi-)omic analysis and guides advanced users to the most recent best practices.

Mapping human adult hippocampal neurogenesis with single-cell transcriptomics: Reconciling controversy or fueling the debate?

The notion of exploiting the regenerative potential of the human brain in physiological aging or neurological diseases represents a particularly attractive alternative to conventional strategies for enhancing or restoring brain function. However, a major first question to address is whether the human brain does possess the ability to regenerate. The existence of human adult hippocampal neurogenesis (AHN) has been at the center of a fierce scientific debate for many years. The advent of single-cell transcriptomic technologies was initially viewed as a panacea to resolving this controversy. However, recent single-cell RNA sequencing studies in the human hippocampus yielded conflicting results. Here, we critically discuss and re-analyze previously published AHN-related single-cell transcriptomic datasets. We argue that, although promising, the single-cell transcriptomic profiling of AHN in the human brain can be confounded by methodological, conceptual, and biological factors that need to be consistently addressed across studies and openly discussed within the scientific community.

Transcriptome-Powered Pluripotent Stem Cell Differentiation for Regenerative Medicine

Pluripotent stem cells are endless sources for in vitro engineering human tissues for regenerative medicine. Extensive studies have demonstrated that transcription factors are the key to stem cell lineage commitment and differentiation efficacy. As the transcription factor profile varies depending on the cell type, global transcriptome analysis through RNA sequencing (RNAseq) has been a powerful tool for measuring and characterizing the success of stem cell differentiation. RNAseq has been utilized to comprehend how gene expression changes as cells differentiate and provide a guide to inducing cellular differentiation based on promoting the expression of specific genes. It has also been utilized to determine the specific cell type. This review highlights RNAseq techniques, tools for RNAseq data interpretation, RNAseq data analytic methods and their utilities, and transcriptomics-enabled human stem cell differentiation. In addition, the review outlines the potential benefits of the transcriptomics-aided discovery of intrinsic factors influencing stem cell lineage commitment, transcriptomics applied to disease physiology studies using patients' induced pluripotent stem cell (iPSC)-derived cells for regenerative medicine, and the future outlook on the technology and its implementation.

Methanol fixation is the method of choice for droplet-based single-cell transcriptomics of neural cells

The main critical step in single-cell transcriptomics is sample preparation. Several methods have been developed to preserve cells after dissociation to uncouple sample handling from library preparation. Yet, the suitability of these methods depends on the cell types to be processed. In this project, we perform a systematic comparison of preservation methods for droplet-based single-cell RNA-seq on neural and glial cells derived from induced pluripotent stem cells. Our results show that while DMSO provides the highest cell quality in terms of RNA molecules and genes detected per cell, it strongly affects the cellular composition and induces the expression of stress and apoptosis genes. In contrast, methanol fixed samples display a cellular composition similar to fresh samples and provide a good cell quality and little expression biases. Taken together, our results show that methanol fixation is the method of choice for performing droplet-based single-cell transcriptomics experiments on neural cell populations.

In vitro erythrocyte production using human-induced pluripotent stem cells: determining the best hematopoietic stem cell sources

BACKGROUND

Blood transfusion is an essential part of medicine. However, many countries have been facing a national blood crisis. To address this ongoing blood shortage issue, there have been efforts to generate red blood cells (RBCs) in vitro, especially from human-induced pluripotent stem cells (hiPSCs). However, the best source of hiPSCs for this purpose is yet to be determined.

METHODS

In this study, hiPSCs were established from three different hematopoietic stem cell sources-peripheral blood (PB), cord blood (CB) and bone marrow (BM) aspirates (n = 3 for each source)-using episomal reprogramming vectors and differentiated into functional RBCs. Various time-course studies including immunofluorescence assay, quantitative real-time PCR, flow cytometry, karyotyping, morphological analysis, oxygen binding capacity analysis, and RNA sequencing were performed to examine and compare the characteristics of hiPSCs and hiPSC-differentiated erythroid cells.

RESULTS

hiPSC lines were established from each of the three sources and were found to be pluripotent and have comparable characteristics. All hiPSCs differentiated into erythroid cells, but there were discrepancies in differentiation and maturation efficiencies: CB-derived hiPSCs matured into erythroid cells the fastest while PB-derived hiPSCs required a longer time for maturation but showed the highest degree of reproducibility. BM-derived hiPSCs gave rise to diverse types of cells and exhibited poor differentiation efficiency. Nonetheless, erythroid cells differentiated from all hiPSC lines mainly expressed fetal and/or embryonic hemoglobin, indicating that primitive erythropoiesis occurred. Their oxygen equilibrium curves were all left-shifted.

CONCLUSIONS

Collectively, both PB- and CB-derived hiPSCs were favorably reliable sources for the clinical production of RBCs in vitro, despite several challenges that need to be overcome. However, owing to the limited availability and the large amount of CB required to produce hiPSCs, and the results of this study, the advantages of using PB-derived hiPSCs for RBC production in vitro may outweigh those of using CB-derived hiPSCs. We believe that our findings will facilitate the selection of optimal hiPSC lines for RBC production in vitro in the near future.

Gene expression differentials driven by mass rearing and artificial selection in black soldier fly colonies

The micro-evolutionary forces that shape genetic diversity during domestication have been assessed in many plant and animal systems. However, the impact of these processes on gene expression, and consequent functional adaptation to artificial environments, remains under-investigated. In this study, whole-transcriptome dynamics associated with the early stages of domestication of the black soldier fly (BSF), Hermetia illucens, were assessed. Differential gene expression (DGE) was evaluated in relation to (i) generational time within the cultured environment (F2 vs. F3), and (ii) two selection strategies [no artificial selective pressure (NS); and selection for greater larval mass (SEL)]. RNA-seq was conducted on 5th instar BSF larvae (n = 36), representing equal proportions of the NS (F2 = 9; F3 = 9) and SEL (F2 = 9; F3 = 9) groups. A multidimensional scaling plot revealed greater gene expression variability within the NS and F2 subgroups, while the SEL group clustered separately with lower levels of variation. Comparisons between generations revealed 898 differentially expressed genes (DEGs; FDR-corrected p < 0.05), while between selection strategies, 213 DEGs were observed (FDR-corrected p < 0.05). Enrichment analyses revealed that metabolic, developmental, and defence response processes were over-expressed in the comparison between F2 and F3 larvae, while metabolic processes were the main differentiating factor between NS and SEL lines. This illustrates the functional adaptations that occur in BSF colonies across generations due to mass rearing; as well as highlighting genic dynamics associated with artificial selection for production traits that might inform future selective breeding strategies.

The shaky foundations of simulating single-cell RNA sequencing data

BACKGROUND

With the emergence of hundreds of single-cell RNA-sequencing (scRNA-seq) datasets, the number of computational tools to analyze aspects of the generated data has grown rapidly. As a result, there is a recurring need to demonstrate whether newly developed methods are truly performant-on their own as well as in comparison to existing tools. Benchmark studies aim to consolidate the space of available methods for a given task and often use simulated data that provide a ground truth for evaluations, thus demanding a high quality standard results credible and transferable to real data.

RESULTS

Here, we evaluated methods for synthetic scRNA-seq data generation in their ability to mimic experimental data. Besides comparing gene- and cell-level quality control summaries in both one- and two-dimensional settings, we further quantified these at the batch- and cluster-level. Secondly, we investigate the effect of simulators on clustering and batch correction method comparisons, and, thirdly, which and to what extent quality control summaries can capture reference-simulation similarity.

CONCLUSIONS

Our results suggest that most simulators are unable to accommodate complex designs without introducing artificial effects, they yield over-optimistic performance of integration and potentially unreliable ranking of clustering methods, and it is generally unknown which summaries are important to ensure effective simulation-based method comparisons.

Sensitive cluster-free differential expression testing

Comparing molecular features, including the identification of genes with differential expression (DE) between conditions, is a powerful approach for characterising disease-specific phenotypes. When testing for DE in single-cell RNA sequencing data, current pipelines first assign cells into discrete clusters (or cell types), followed by testing for differences within each cluster. Consequently, the sensitivity and specificity of DE testing are limited and ultimately dictated by the granularity of the cell type annotation, with discrete clustering being especially suboptimal for continuous trajectories. To overcome these limitations, we present miloDE - a cluster-free framework for differential expression testing. We build on the Milo approach, introduced for differential cell abundance testing, which leverages the graph representation of single-cell data to assign relatively homogenous, 'neighbouring' cells into overlapping neighbourhoods. We address key differences between differential abundance and expression testing at the level of neighbourhood assignment, statistical testing, and multiple testing correction. To illustrate the performance of miloDE we use both simulations and real data, in the latter case identifying a transient haemogenic endothelia-like state in chimeric mouse embryos lacking Tal1 as well as uncovering distinct transcriptional programs that characterise changes in macrophages in patients with Idiopathic Pulmonary Fibrosis. miloDE is available as an open-source R package at https://github.com/MarioniLab/miloDE.

Single-cell gene set enrichment analysis and transfer learning for functional annotation of scRNA-seq data

Although an essential step, cell functional annotation often proves particularly challenging from single-cell transcriptional data. Several methods have been developed to accomplish this task. However, in most cases, these rely on techniques initially developed for bulk RNA sequencing or simply make use of marker genes identified from cell clustering followed by supervised annotation. To overcome these limitations and automatize the process, we have developed two novel methods, the single-cell gene set enrichment analysis (scGSEA) and the single-cell mapper (scMAP). scGSEA combines latent data representations and gene set enrichment scores to detect coordinated gene activity at single-cell resolution. scMAP uses transfer learning techniques to re-purpose and contextualize new cells into a reference cell atlas. Using both simulated and real datasets, we show that scGSEA effectively recapitulates recurrent patterns of pathways' activity shared by cells from different experimental conditions. At the same time, we show that scMAP can reliably map and contextualize new single-cell profiles on a breast cancer atlas we recently released. Both tools are provided in an effective and straightforward workflow providing a framework to determine cell function and significantly improve annotation and interpretation of scRNA-seq data.

Molecular phenotyping of single pancreatic islet leader beta cells by "Flash-Seq"

AIMS

Spatially-organized increases in cytosolic Ca²⁺ within pancreatic beta cells in the pancreatic islet underlie the stimulation of insulin secretion by high glucose. Recent data have revealed the existence of subpopulations of beta cells including "leaders" which initiate Ca²⁺ waves. Whether leader cells possess unique molecular features, or localisation, is unknown.

MAIN METHODS

High speed confocal Ca²⁺ imaging was used to identify leader cells and connectivity analysis, running under MATLAB and Python, to identify highly connected "hub" cells. To explore transcriptomic differences between beta cell sub-groups, individual leaders or followers were labelled by photo-activation of the cryptic fluorescent protein PA-mCherry and subjected to single cell RNA sequencing ("Flash-Seq").

KEY FINDINGS

Distinct Ca²⁺ wave types were identified in individual islets, with leader cells present in 73 % (28 of 38 islets imaged). Scale-free, power law-adherent behaviour was also observed in 29 % of islets, though "hub" cells in these islets did not overlap with leaders. Transcripts differentially expressed (295; padj < 0.05) between leader and follower cells included genes involved in cilium biogenesis and transcriptional regulation. Providing some support for these findings, ADCY6 immunoreactivity tended to be higher in leader than follower cells, whereas cilia number and length tended to be lower in the former. Finally, leader cells were located significantly closer to delta, but not alpha, cells in Euclidian space than were follower cells.

SIGNIFICANCE

The existence of both a discrete transcriptome and unique localisation implies a role for these features in defining the specialized function of leaders. These data also raise the possibility that localised signalling between delta and leader cells contributes to the initiation and propagation of islet Ca²⁺ waves.

De-novo reconstruction and identification of transcriptional gene regulatory network modules differentiating single-cell clusters

Single-cell RNA sequencing (scRNA-seq) technology provides an unprecedented opportunity to understand gene functions and interactions at single-cell resolution. While computational tools for scRNA-seq data analysis to decipher differential gene expression profiles and differential pathway expression exist, we still lack methods to learn differential regulatory disease mechanisms directly from the single-cell data. Here, we provide a new methodology, named DiNiro, to unravel such mechanisms de novo and report them as small, easily interpretable transcriptional regulatory network modules. We demonstrate that DiNiro is able to uncover novel, relevant, and deep mechanistic models that not just predict but explain differential cellular gene expression programs. DiNiro is available at https://exbio.wzw.tum.de/diniro/.

The GRN concept as a guide for evolutionary developmental biology

Organismal phenotypes result largely from inherited developmental programs, usually executed during embryonic and juvenile life stages. These programs are not blank slates onto which natural selection can draw arbitrary forms. Rather, the mechanisms of development play an integral role in shaping phenotypic diversity and help determine the evolutionary trajectories of species. Modern evolutionary biology must, therefore, account for these mechanisms in both theory and in practice. The gene regulatory network (GRN) concept represents a potent tool for achieving this goal whose utility has grown in tandem with advances in "omic" technologies and experimental techniques. However, while the GRN concept is widely utilized, it is often less clear what practical implications it has for conducting research in evolutionary developmental biology. In this Perspective, we attempt to provide clarity by discussing how experiments and projects can be designed in light of the GRN concept. We first map familiar biological notions onto the more abstract components of GRN models. We then review how diverse functional genomic approaches can be directed toward the goal of constructing such models and discuss current methods for functionally testing evolutionary hypotheses that arise from them. Finally, we show how the major steps of GRN model construction and experimental validation suggest generalizable workflows that can serve as a scaffold for project design. Taken together, the practical implications that we draw from the GRN concept provide a set of guideposts for studies aiming at unraveling the molecular basis of phenotypic diversity.

TEAM: A MULTIPLE TESTING ALGORITHM ON THE AGGREGATION TREE FOR FLOW CYTOMETRY ANALYSIS

In immunology studies, flow cytometry is a commonly used multivariate single-cell assay. One key goal in flow cytometry analysis is to detect the immune cells responsive to certain stimuli. Statistically, this problem can be translated into comparing two protein expression probability density functions (pdfs) before and after the stimulus; the goal is to pinpoint the regions where these two pdfs differ. Further screening of these differential regions can be performed to identify enriched sets of responsive cells. In this paper, we model identifying differential density regions as a multiple testing problem. First, we partition the sample space into small bins. In each bin, we form a hypothesis to test the existence of differential pdfs. Second, we develop a novel multiple testing method, called TEAM (Testing on the Aggregation tree Method), to identify those bins that harbor differential pdfs while controlling the false discovery rate (FDR) under the desired level. TEAM embeds the testing procedure into an aggregation tree to test from fine- to coarse-resolution. The procedure achieves the statistical goal of pinpointing density differences to the smallest possible regions. TEAM is computationally efficient, capable of analyzing large flow cytometry data sets in much shorter time compared with competing methods. We applied TEAM and competing methods on a flow cytometry data set to identify T cells responsive to the cytomegalovirus (CMV)-pp65 antigen stimulation. With additional downstream screening, TEAM successfully identified enriched sets containing monofunctional, bifunctional, and polyfunctional T cells. Competing methods either did not finish in a reasonable time frame or provided less interpretable results. Numerical simulations and theoretical justifications demonstrate that TEAM has asymptotically valid, powerful, and robust performance. Overall, TEAM is a computationally efficient and statistically powerful algorithm that can yield meaningful biological insights in flow cytometry studies.

Distinct Mechanisms of Mismatch-Repair Deficiency Delineate Two Modes of Response to Anti-PD-1 Immunotherapy in Endometrial Carcinoma

Mismatch repair-deficient (MMRd) cancers have varied responses to immune-checkpoint blockade (ICB). We conducted a phase II clinical trial of the PD-1 inhibitor pembrolizumab in 24 patients with MMRd endometrial cancer (NCT02899793). Patients with mutational MMRd tumors (6 patients) had higher response rates and longer survival than those with epigenetic MMRd tumors (18 patients). Mutation burden was higher in tumors with mutational MMRd compared with epigenetic MMRd; however, within each category of MMRd, mutation burden was not correlated with ICB response. Pretreatment JAK1 mutations were not associated with primary resistance to pembrolizumab. Longitudinal single-cell RNA-seq of circulating immune cells revealed contrasting modes of antitumor immunity for mutational versus epigenetic MMRd cancers. Whereas effector CD8+ T cells correlated with regression of mutational MMRd tumors, activated CD16+ NK cells were associated with ICB-responsive epigenetic MMRd tumors. These data highlight the interplay between tumor-intrinsic and tumor-extrinsic factors that influence ICB response.

SIGNIFICANCE

The molecular mechanism of MMRd is associated with response to anti-PD-1 immunotherapy in endometrial carcinoma. Tumors with epigenetic MMRd or mutational MMRd are correlated with NK cell or CD8+ T cell-driven immunity, respectively. Classifying tumors by the mechanism of MMRd may inform clinical decision-making regarding cancer immunotherapy. This article is highlighted in the In This Issue feature, p. 247.

An Iterative Unsupervised Method for Gene Expression Differentiation

For several decades, intensive research for understanding gene activity and its role in organism's lives is the research focus of scientists in different areas. A part of these investigations is the analysis of gene expression data for selecting differentially expressed genes. Methods that identify the interested genes have been proposed on statistical data analysis. The problem is that there is no good agreement among them, as different results are produced by distinct methods. By taking the advantage of the unsupervised data analysis, an iterative clustering procedure that finds differentially expressed genes shows promising results. In the present paper, a comparative study of the clustering methods applied for gene expression analysis is presented to explicate the choice of the clustering algorithm implemented in the method. An investigation of different distance measures is provided to reveal those that increase the efficiency of the method in finding the real data structure. Further, the method is improved by incorporating an additional aggregation measure based on the standard deviation of the expression levels. Its usage increases the gene distinction as a new amount of differentially expressed genes is found. The method is summarized in a detailed procedure. The significance of the method is proved by an analysis of two mice strain data sets. The differentially expressed genes defined by the proposed method are compared with those selected by the well-known statistical methods applied to the same data set.

Organism-Wide Analysis of Sepsis Reveals Mechanisms of Systemic Inflammation

Sepsis is a systemic response to infection with life-threatening consequences. Our understanding of the impact of sepsis across organs of the body is rudimentary. Here, using mouse models of sepsis, we generate a dynamic, organism-wide map of the pathogenesis of the disease, revealing the spatiotemporal patterns of the effects of sepsis across tissues. These data revealed two interorgan mechanisms key in sepsis. First, we discover a simplifying principle in the systemic behavior of the cytokine network during sepsis, whereby a hierarchical cytokine circuit arising from the pairwise effects of TNF plus IL-18, IFN-γ, or IL-1β explains half of all the cellular effects of sepsis on 195 cell types across 9 organs. Second, we find that the secreted phospholipase PLA2G5 mediates hemolysis in blood, contributing to organ failure during sepsis. These results provide fundamental insights to help build a unifying mechanistic framework for the pathophysiological effects of sepsis on the body.

XCVATR: detection and characterization of variant impact on the Embeddings of single -cell and bulk RNA-sequencing samples

BACKGROUND

RNA-sequencing has become a standard tool for analyzing gene activity in bulk samples and at the single-cell level. By increasing sample sizes and cell counts, this technique can uncover substantial information about cellular transcriptional states. Beyond quantification of gene expression, RNA-seq can be used for detecting variants, including single nucleotide polymorphisms, small insertions/deletions, and larger variants, such as copy number variants. Notably, joint analysis of variants with cellular transcriptional states may provide insights into the impact of mutations, especially for complex and heterogeneous samples. However, this analysis is often challenging due to a prohibitively high number of variants and cells, which are difficult to summarize and visualize. Further, there is a dearth of methods that assess and summarize the association between detected variants and cellular transcriptional states.

RESULTS

Here, we introduce XCVATR (eXpressed Clusters of Variant Alleles in Transcriptome pRofiles), a method that identifies variants and detects local enrichment of expressed variants within embedding of samples and cells in single-cell and bulk RNA-seq datasets. XCVATR visualizes local "clumps" of small and large-scale variants and searches for patterns of association between each variant and cellular states, as described by the coordinates of cell embedding, which can be computed independently using any type of distance metrics, such as principal component analysis or t-distributed stochastic neighbor embedding. Through simulations and analysis of real datasets, we demonstrate that XCVATR can detect enrichment of expressed variants and provide insight into the transcriptional states of cells and samples. We next sequenced 2 new single cell RNA-seq tumor samples and applied XCVATR. XCVATR revealed subtle differences in CNV impact on tumors.

CONCLUSIONS

XCVATR is publicly available to download from https://github.com/harmancilab/XCVATR .

Dysfunctional Sars-CoV-2-M protein-specific cytotoxic T lymphocytes in patients recovering from severe COVID-19

Although the importance of virus-specific cytotoxic T lymphocytes (CTL) in virus clearance is evident in COVID-19, the characteristics of virus-specific CTLs related to disease severity have not been fully explored. Here we show that the phenotype of virus-specific CTLs against immunoprevalent epitopes in COVID-19 convalescents might differ according to the course of the disease. We establish a cellular screening method that uses artificial antigen presenting cells, expressing HLA-A^*24:02, the costimulatory molecule 4-1BBL, SARS-CoV-2 structural proteins S, M, and N and non-structural proteins ORF3a and nsp6/ORF1a. The screen implicates SARS-CoV-2 M protein as a frequent target of IFNγ secreting CD8⁺ T cells, and identifies M_198-206 as an immunoprevalent epitope in our cohort of HLA-A^*24:02 positive convalescent COVID-19 patients recovering from mild, moderate and severe disease. Further exploration of M_198-206-specific CD8⁺ T cells with single cell RNA sequencing reveals public TCRs in virus-specific CD8⁺ T cells, and shows an exhausted phenotype with less differentiated status in cells from the severe group compared to cells from the moderate group. In summary, this study describes a method to identify T cell epitopes, indicate that dysfunction of virus-specific CTLs might be an important determinant of clinical outcomes.

Single-cell profiling of healthy human kidney reveals features of sex-based transcriptional programs and tissue-specific immunity

Knowledge of the transcriptional programs underpinning the functions of human kidney cell populations at homeostasis is limited. We present a single-cell perspective of healthy human kidney from 19 living donors, with equal contribution from males and females, profiling the transcriptome of 27677 cells to map human kidney at high resolution. Sex-based differences in gene expression within proximal tubular cells were observed, specifically, increased anti-oxidant metallothionein genes in females and aerobic metabolism-related genes in males. Functional differences in metabolism were confirmed in proximal tubular cells, with male cells exhibiting higher oxidative phosphorylation and higher levels of energy precursor metabolites. We identified kidney-specific lymphocyte populations with unique transcriptional profiles indicative of kidney-adapted functions. Significant heterogeneity in myeloid cells was observed, with a MRC1+LYVE1+FOLR2+C1QC+ population representing a predominant population in healthy kidney. This study provides a detailed cellular map of healthy human kidney, and explores the complexity of parenchymal and kidney-resident immune cells.

Circuit-selective cell-autonomous regulation of inhibition in pyramidal neurons by Ste20-like kinase

Maintaining an appropriate balance between excitation and inhibition is critical for neuronal information processing. Cortical neurons can cell-autonomously adjust the inhibition they receive to individual levels of excitatory input, but the underlying mechanisms are unclear. We describe that Ste20-like kinase (SLK) mediates cell-autonomous regulation of excitation-inhibition balance in the thalamocortical feedforward circuit, but not in the feedback circuit. This effect is due to regulation of inhibition originating from parvalbumin-expressing interneurons, while inhibition via somatostatin-expressing interneurons is unaffected. Computational modeling shows that this mechanism promotes stable excitatory-inhibitory ratios across pyramidal cells and ensures robust and sparse coding. Patch-clamp RNA sequencing yields genes differentially regulated by SLK knockdown, as well as genes associated with excitation-inhibition balance participating in transsynaptic communication and cytoskeletal dynamics. These data identify a mechanism for cell-autonomous regulation of a specific inhibitory circuit that is critical to ensure that a majority of cortical pyramidal cells participate in information coding.