Selecting gene features for unsupervised analysis of single-cell gene expression data.

Motivation: Cell populations are never truly homogeneous; individual cells exist in biochemical states that define functional differences between them. New technology based on microfluidic arrays combined with multiplexed quantitative polymerase chain reactions now enables high-throughput single-cell gene expression measurement, allowing assessment of cellular heterogeneity. However, few analytic tools have been developed specifically for the statistical and analytical challenges of single-cell quantitative polymerase chain reactions data.Results: We present a statistical framework for the exploration, quality control and analysis of single-cell gene expression data from microfluidic arrays. We assess accuracy and within-sample heterogeneity of single-cell expression and develop quality control criteria to filter unreliable cell measurements. We propose a statistical model accounting for the fact that genes at the single-cell level can be on (and a continuous expression measure is recorded) or dichotomously off (and the recorded expression is zero). Based on this model, we derive a combined likelihood ratio test for differential expression that incorporates both the discrete and continuous components. Using an experiment that examines treatment-specific changes in expression, we show that this combined test is more powerful than either the continuous or dichotomous component in isolation, or a t-test on the zero-inflated data. Although developed for measurements from a specific platform (Fluidigm), these tools are generalizable to other multi-parametric measures over large numbers of events.Availability: All results presented here were obtained using the SingleCellAssay R package available on GitHub (http://github.com/RGLab/SingleCellAssay).Contact:rgottard@fhcrc.orgSupplementary information:Supplementary data are available at Bioinformatics online.

Analysis of single-cell gene expression data and genome assemblies of five diverse metazoan species shows a tight link between conserved gene order and its relationship with cell type-specific gene expression.

We present single-cell clustering using bifurcation analysis (SCUBA), a novel computational method for extracting lineage relationships from single-cell gene expression data and modeling the dynamic changes associated with cell differentiation. SCUBA draws techniques from nonlinear dynamics and stochastic differential equation theories, providing a systematic framework for modeling complex processes involving multilineage specifications. By applying SCUBA to analyze two complementary, publicly available datasets we successfully reconstructed the cellular hierarchy during early development of mouse embryos, modeled the dynamic changes in gene expression patterns, and predicted the effects of perturbing key transcriptional regulators on inducing lineage biases. The results were robust with respect to experimental platform differences between RT-PCR and RNA sequencing. We selectively tested our predictions in Nanog mutants and found good agreement between SCUBA predictions and the experimental data. We further extended the utility of SCUBA by developing a method to reconstruct missing temporal-order information from a typical single-cell dataset. Analysis of a hematopoietic dataset suggests that our method is effective for reconstructing gene expression dynamics during human B-cell development. In summary, SCUBA provides a useful single-cell data analysis tool that is well-suited for the investigation of developmental processes.

ABSTRACTFoundation Models (FMs) have made significant strides in both industrial and scientific domains. In this paper, we evaluate the performance of FMs for single‐cell sequencing data analysis through comprehensive experiments across eight downstream tasks pertinent to single‐cell data. Overall, the top FMs include scGPT, Geneformer, and CellFM by considering model performances and user accessibility among ten single‐cell FMs. However, by comparing these FMs with task‐specific methods, we found that single‐cell FMs may not consistently excel than task‐specific methods in all tasks, which challenges the necessity of developing foundation models for single‐cell analysis. In addition, we evaluated the effects of hyperparameters, initial settings, and stability for training single‐cell FMs based on a proposed scEval framework, and provide guidelines for pre‐training and fine‐tuning to enhance the performances of single‐cell FMs. Our work summarizes the current state of single‐cell FMs, points to their constraints and avenues for future development, and offers a freely available evaluation pipeline to benchmark new models and improve method development.

BackgroundSerial Analysis of Gene Expression (SAGE) is a DNA sequencing-based method for large-scale gene expression profiling that provides an alternative to microarray analysis. Most analyses of SAGE data aimed at identifying co-expressed genes have been accomplished using various versions of clustering approaches that often result in a number of false positives.Principal FindingsHere we explore the use of seriation, a statistical approach for ordering sets of objects based on their similarity, for large-scale expression pattern discovery in SAGE data. For this specific task we implement a seriation heuristic we term ‘progressive construction of contigs’ that constructs local chains of related elements by sequentially rearranging margins of the correlation matrix. We apply the heuristic to the analysis of simulated and experimental SAGE data and compare our results to those obtained with a clustering algorithm developed specifically for SAGE data. We show using simulations that the performance of seriation compares favorably to that of the clustering algorithm on noisy SAGE data.ConclusionsWe explore the use of a seriation approach for visualization-based pattern discovery in SAGE data. Using both simulations and experimental data, we demonstrate that seriation is able to identify groups of co-expressed genes more accurately than a clustering algorithm developed specifically for SAGE data. Our results suggest that seriation is a useful method for the analysis of gene expression data whose applicability should be further pursued.

Background - In the NeoTRIPaPDL1 phase III trial, triple-negative breast cancer (TNBC) patients were randomized to receive nab-paclitaxel/carboplatin for 8 cycles (CT arm) with or without atezolizumab (CT/A arm). We previously reported that the majority of patients (65%) with PD-L1- at baseline converted to PD-L1+ (the majority with IC2/IC3) after the first cycle of treatment in the atezolizumab arm. Here we studied treatment associated changes in the molecular tumour features and immune microenvironment by PD-L1 expression. Methods – A total of 158 (56.4%) of patients enrolled in the NeoTRIPaPDL1 trial (76 and 82 from CT/A and CT arm, respectively) were included in this analysis by satisfying the following criteria: i) availability of a paired core biopsy at baseline and at day 1 of the second treatment cycle (D1C2); ii) evaluation of stromal TILs (sTILs) and staining for PD-L1 (Ventana SP142); iii) successful RNA-seq gene expression profiling. PD-L1 groups were defined as IC0/IC1 (<5% PD-L1low) vs IC2/IC3 (>=5% PD-L1high). Presence of different immune cell populations was quantified by gene expression profile deconvolution using ConsensusTME R package. Cancer hallmark gene set collection and custom signature activation status were estimated in each sample using singscore R package. Differences in score distribution were evaluated by 2-sided t-test. Results At baseline, both sTILs and immune cell signatures were upregulated in PD-L1high compared to PD-L1low tumours (p<0.001). No significant differences were found between the two treatment arms in the PD-L1high subpopulation. In the PD-L1low cases, both sTILs and immune related signatures were slightly downregulated in the CT/A arm (p<0.05), suggesting a modest unbalance among treatment arms with CT arm being slightly more inflamed than CT/A arm.At D1C2, PD-L1high tumours in the CT arm systematically had high sTILs (median=70%, range=30-90%), while PD-L1low tumours receiving CT/A had low sTILs in a significant proportion of cases (median=30%, range=0-90%; p<0.001). Similarly, at D1C2 several gene expression-estimated immune cell populations and immune-related signatures were upregulated in the CT arm compared to CT/A arm in PD-L1high tumours, with the weakest association observed for M2 macrophages (p=0.058). No tumour-related signatures were differentially expressed among the two treatment arms within groups with PD-L1high, suggesting that different treatment. modulate PD-L1 by engaging a different immune mileau instead of modulating tumor related features.Considering PD-L1low groups at D1C2, in the CT arm 21 tumors had >30% sTILs (n=21/69, 30.4%), while in the CT/A arm only 2 had sTILs >30% (2/27, 7.4%) (p<0.001). Analysis of gene expression data identified IFN-related signatures as the most upregulated in the CT compared to CT/A arm in PD-L1low cases at D1C2 (p<0.05). Conclusions Integrated dynamic analysis of PD-L1 expression and gene expression data highlighted significant treatment-specific changes of the immune landscape according to PD-L1 expression, when this biomarker is assessed during treatment. This indicates that the immune milieu associated with PD-L1 status is strongly dependent on the time of assessment in relationship to treatment received. Such observation may explain why in the NeoTRIP trial, baseline PD-L1 but not on-treatment PD-L1 was predictive of pCR in CT/A arm. In addition, our findings could have implications related to the use of PD-L1 as a predictive biomarker in pre-treated patients, especially when assessed early on during treatment. Citation Format: Maurizio Callari, Chiun-Sheng Huang, Daniel Egle, Begoña Bermejo, Claudio Zamagni, Matteo Dugo, Marc Thill, Antonio Anton, Marco Barreca, Stefania Russo, Eva Maria Ciruelos, Richard Greil, Stefania Zambelli, Balázs Gyorffy, Chanel Smart, Olivia Biasi, Pinuccia Valagussa, Giuseppe Viale, Luca Gianni, Giampaolo Bianchini. Immune milieu associated with PD-L1 status in TNBC is dependent on time of biomarker assessment and treatment received: A secondary analysis of the NeoTRIPaPDL1 trial [abstract]. In: Proceedings of the 2021 San Antonio Breast Cancer Symposium; 2021 Dec 7-10; San Antonio, TX. Philadelphia (PA): AACR; Cancer Res 2022;82(4 Suppl):Abstract nr P1-04-02.

Tumor-associated macrophages (TAMs) constitute up to 50% of tumor bulk in glioblastoma (GBM) and play an important role in tumor maintenance and progression. The recently discovered differences between invading tumour periphery and hypoxic tumor core implies that macrophage biology is also distinct by location. This may provide further insight into the observed treatment resistance to immune modulation. We hypothesize that macrophage activation occurs through processes that are distinct in tumor periphery versus core. We therefore investigated regional differences in TAM recruitment and evolution in GBM by combining open source single cell and bulk gene expression data. We used single cell gene expression data from 4 glioblastomas (total of 3589 cells) and 122 total bulk samples obtained from 10 different patients. Cell identity, ontogeny (bone-marrow derived macrophages-BMDM vs microglia), and macrophage activation state were inferred using verified gene expression signatures. We captured the spectrum of immune states using cell trajectory analysis with pseudotime ordering. In keeping with previous studies, TAMs carrying BMDM identity were more abundant in tumor bulk while microglia-derived TAMs dominated the tumor periphery across all macrophage activation states including pre-activation. We note that core TAMs evolve towards a pro-inflammatory state and identify a subpopulation of cells based on a gene program exhibiting strong, opposing correlation with Programmed cell Death-1 (PD-1) signaling, which may correlate to their response to PD-1 inhibition. By contrast, peripheral TAMs evolve towards anti-inflammatory phenotype and contains a population of cells strongly associated with NFkB signaling. Our preliminary analysis suggests important regional differences in TAMs with regard to recruitment and evolution. We identify regionally distinct and potentially actionable cell subpopulations and advocate the need for a multi-targeted approach to GBM therapeutics.

BackgroundMicroarrays have become extremely useful for analysing genetic phenomena, but establishing a relation between microarray analysis results (typically a list of genes) and their biological significance is often difficult. Currently, the standard approach is to map a posteriori the results onto gene networks in order to elucidate the functions perturbed at the level of pathways. However, integrating a priori knowledge of the gene networks could help in the statistical analysis of gene expression data and in their biological interpretation.ResultsWe propose a method to integrate a priori the knowledge of a gene network in the analysis of gene expression data. The approach is based on the spectral decomposition of gene expression profiles with respect to the eigenfunctions of the graph, resulting in an attenuation of the high-frequency components of the expression profiles with respect to the topology of the graph. We show how to derive unsupervised and supervised classification algorithms of expression profiles, resulting in classifiers with biological relevance. We illustrate the method with the analysis of a set of expression profiles from irradiated and non-irradiated yeast strains.ConclusionIncluding a priori knowledge of a gene network for the analysis of gene expression data leads to good classification performance and improved interpretability of the results.

Recent research has focused on incorporating biological function and pathway information into the analysis of gene expression data, partly as a means of compensating for insufficient experimental replications, low signal to noise, lack of reproducibility and/or multiple testing confounds. A Bayesian approach seems to be ideal for incorporating functional information into gene expression data analysis. In this study, we tested the feasibility of using literature derived gene relationships in a Bayesian model to analyze gene expression data. Prior distributions were constructed based on gene associations derived from the biomedical literature using Latent Semantic Indexing (LSI). The LSI model was built using more than 1 million Medline abstracts corresponding to 22,000 human and mouse genes. A key advantage of LSI is that both explicit and implicit gene relationships can be derived from the literature. Gene neighborhoods were determined using latent Gaussian Markov random fields and logistic transformation of the latent variables. We tested the procedure on a microarray dataset for interferon-stimulated genes in mouse embryonic fibroblasts. By integrating functional information from literature, Bayesian approach identified relevant genes that previously did not meet the 0.05 significance level. In comparison to a standard mixture model, spatial mixture model has more power for identifying direct and indirect interferon regulated genes. The spatial model enhanced the ranks of some genes which are known to be affected by interferon treatment, such as Nmi (NMI N-myc and STAT interactor) and ifi35 (interferon-induced protein 35). It also identified some genes that previously were ignored because of the marginal p-values, such as dpysl2, map2k1, msn, Psck5, and Il6st. Interestingly, these genes appear to be indirectly related to interferon treatment. In summary, we show that our procedure increases statistical power and produces more biologically meaningful gene lists. These results suggest that Bayesian methods which incorporate functional information from the literature may improve analysis of gene expression data.

The transition between epithelial and mesenchymal states has fundamental importance for embryonic development, stem cell reprogramming, and cancer progression. Here, we construct a topographic map underlying epithelial-mesenchymal transitions using a combination of numerical simulations of a Boolean network model and the analysis of bulk and single-cell gene expression data. The map reveals a multitude of metastable hybrid phenotypic states, separating stable epithelial and mesenchymal states, and is reminiscent of the free energy measured in glassy materials and disordered solids. Our work not only elucidates the nature of hybrid mesenchymal/epithelial states but also provides a general strategy to construct a topographic representation of phenotypic plasticity from gene expression data using statistical physics methods.

【In DNA microarray studies, the number of genes far exceeds the number of samples and the gene expression measures are highly correlated. Partial least squares regression (PLSR) is one of the popular methods for dimensional reduction and known to be useful for the classifications of microarray data by several studies. In this study, we suggest a modified version of the partial least squares regression to analyze gene expression data with survival information. The method is designed as a new gene selection method using PLSR with an iterative procedure of imputing censored survival time. Mean square error of prediction criterion is used to determine the dimension of the model. To visualize the data, plot for variables superimposed with samples are used. The method is applied to two microarray data sets, both containing survival time. The results show that the proposed method works well for interpreting gene expression microarray data.】

The analysis of gene expression data has made significant contributions to understanding disease mechanisms and developing new drugs and therapies. In such analysis, gene selection is often required for identifying informative and relevant genes and removing redundant and irrelevant ones. However, this is not an easy task as gene expression data have inherent challenges such as ultra-high dimensionality, biological noise, and measurement errors. This study focuses on the measurement errors in gene selection problems. Typically, high-throughput experiments have their own intrinsic measurement errors, which can result in an increase of falsely discovered genes. To alleviate this problem, this study proposes a gene selection method that takes into account measurement errors using generalized liner measurement error models. The method consists of iterative filtering and selection steps until convergence, leading to fewer false positives and providing stable results under measurement errors. The performance of the proposed method is demonstrated through simulation studies and applied to a lung cancer data set.

Cytoscape is a bioinformatic data analysis and visualization platform that is well-suited to the analysis of gene expression data. To facilitate the analysis of yeast microarray data using Cytoscape, we constructed an interaction network (interactome) using the curated interaction data available from the Saccharomyces Genome Database (www.yeastgenome.org) and the database of yeast transcription factors at YEASTRACT (www.yeastract.com). These data were formatted and imported into Cytoscape using semi-automated methods, including Linux-based scripts, that simplified the process while minimizing the introduction of processing errors. The methods described for the construction of this yeast interactome are generally applicable to the construction of any interactome. Using Cytoscape, we illustrate the use of this interactome through the analysis of expression data from a recent yeast diauxic shift experiment. We also report and briefly describe the complex associations among transcription factors that result in the regulation of thousands of genes through coordinated changes in expression of dozens of transcription factors. These cells are thus able to sensitively regulate cellular metabolism in response to changes in genetic or environmental conditions through relatively small changes in the expression of large numbers of genes, affecting the entire yeast metabolome.

As single-cell gene expression data analysis continues to grow, the need for reliable clustering techniques has become increasingly important. This growth has also led to a rise in heuristic means for technique choice, which could lead to inaccurate reports if a comprehensive evaluation of the resulting clusters is omitted. Previous work in the field of microarray data provided measures of stability to evaluate clustering techniques for highly correlated data. Additional work on aggregation in the same era presented a way to combine ranked lists using several performance measures. In this paper, we adapt the stability measures from the microarray era to function on single-cell data and employ them alongside several existing measures to evaluate characteristics of clustering techniques at varying parameter choices. We then compare individual measures alongside an aggregation method (AM) on all measures to rank these techniques for five datasets with known groups and one without known groups. Our findings show that techniques ranked highest by a single validation measure tend to have mediocre performance for other measures for multiple datasets. Alternatively, our method of aggregating techniques based on their ranked performance on all measures ensures that techniques with a strong overall performance are ranked highest. The techniques selected in this fashion had above average performance in most of the measured clustering characteristics and perform well in other ways, such as consistently estimating cluster counts that closely match the true number of biological groups.

Identifying common and distinctive processes underlying multiset data