Cellograph: a semi-supervised approach to analyzing multi-condition single-cell RNA-sequencing data using graph neural networks

Abstract

With the growing number of single-cell datasets collected under more complex experimental conditions, there is an opportunity to leverage single-cell variability to reveal deeper insights into how cells respond to perturbations. Many existing approaches rely on discretizing the data into clusters for differential gene expression (DGE), effectively ironing out any information unveiled by the single-cell variability across cell-types. In addition, DGE often assumes a statistical distribution that, if erroneous, can lead to false positive differentially expressed genes. Here, we present Cellograph: a semi-supervised framework that uses graph neural networks to quantify the effects of perturbations at single-cell granularity. Cellograph not only measures how prototypical cells are of each condition but also learns a latent space that is amenable to interpretable data visualization and clustering. The learned gene weight matrix from training reveals pertinent genes driving the differences between conditions. We demonstrate the utility of our approach on publicly-available datasets including cancer drug therapy, stem cell reprogramming, and organoid differentiation. Cellograph outperforms existing methods for quantifying the effects of experimental perturbations and offers a novel framework to analyze single-cell data using deep learning.

Keywords: Graph neural networks; Semi-supervised learning; Single-cell genomics.

Fig. 1

Illustrative overview of Cellograph algorithm. Single-cell data collected from multiple sample drug treatments (A, B) is converted to a kNN graph (C), where cells are nodes, and edges denote connections between transcriptionally similar cells. The colored rectangles (B) correspond to the different samples represented by the drugs in A. This kNN is fed in as input to a two-layer GCN (D) that quantitatively and visually learns how prototypical each cell is of its experimental label through the learned latent embedding. E A mathematical schematic of the first layer, where each cell’s gene expression and its neighbors’s gene expression is aggregated to produce a lower-dimensional representation of the cell in a latent space. F A mathematical schematic of the second layer respectively, where the output embedding of the first layer is mapped to softmax probabilities of cells belonging to each of the drug treatments

Fig. 2

Cellograph identifies treatment groups and distinguishes genes defining these groups on a human organoid dataset. A PHATE projection of learned latent space, with cells colored by treatment labels, probabilities of belonging to control or KPT-treated cells, clusters obtained by k-means clustering of the learned latent embedding with $k=3$, and gene expression of GDF15 and KLK7. B Heatmap of top 25 weighted genes from parameterized gene weight matrix. C Heatmap of differentially expressed genes between clusters derived from Cellograph. D Compositional plot of predicted treatment groups from the softmax probabilities ($z_{\mathit{ij}}>0.5$) (left) and cell types annotated by the original study (right) partitioned by clusters

Fig. 3

Cellograph defines genetic signatures of distinct drug responses in the drug holiday dataset. A PHATE embeddings of the learned latent space colored according to the treatment labels, clusters, and treatment probabilities (day 0 not shown). B Heatmap of top 25 weighted genes from learned parameterized gene weight matrix. C The distribution of treatment probabilities for Day 11 cells partitioned by treatment groups. D The distribution of gene expression between clusters 0, 5, 3, and 2 of select differentially expressed genes (INHBA, TUBA1B)

Fig. 4

Cellograph distinguishes the molecular mechanisms of transdifferentiation and dedifferentiation in myogenesis. A PHATE embeddings of learned latent space annotated according to treatment conditions, clusters, and softmax probabilities of all conditions except for MEFs, defining the in-group variation. B Heatmap of top weighted genes from parameterized gene weight matrix, identifying pertinent genes such as cyclin D1 and CRABP1. C Violin plot of softmax probabilities of cells belonging to the MyoD/day 4 treatment group, showing similarities to the MyoD/day 2 population. D Violin plots of top 20 differentially expressed genes between clusters 1 and 8 and clusters 3 and 9, which define the ${\text{Pax7}}^{+}$ cells and MyoD+FRC/day 8 treated cells, respectively. E Compositional plot of predicted cell types partitioned by cluster

Fig. 5

Results of running Milo and CNA on the datasets evaluated. A Output of running Milo and CNA on the human organoid dataset. B Output of running Milo and CNA on the drug holiday dataset. C Output of running Milo and CNA on the myogenesis dataset

Fig. 6

Runtime of Cellograph’s performance versus MELD’s on optimal parameter settings. Cellograph consistently outperforms MELD on each dataset, while using fewer computing resources (y-axis is log-scaled)

Fig. 7

Boxplots of NMI values per clustering algorithm. Distributions of 100 independent NMI calculations for each clustering algorithm for all three datasets evaluated, quantifying concordance between the cluster assignments and ground truth labels