Multicondition and multimodal temporal profile inference during mouse embryonic development

Abstract

The emergence of single-cell time-series data sets enables modeling of changes in various types of cellular profiles over time. However, because of the disruptive nature of single-cell measurements, it is impossible to capture the full temporal trajectory of a particular cell. Furthermore, single-cell profiles can be collected at mismatched time points across different conditions (e.g., sex, batch, disease) and data modalities (e.g., scRNA-seq, scATAC-seq), which makes modeling challenging. Here, we propose a joint modeling framework, Sunbear, for integrating multicondition and multimodal single-cell profiles across time. Sunbear can be used to impute single-cell temporal profile changes, align multi-data set and multimodal profiles across time, and extrapolate single-cell profiles in a missing modality. We apply Sunbear to reveal sex-biased transcription during mouse embryonic development and predict dynamic relationships between epigenetic priming and transcription for cells in which multimodal profiles are unavailable. Sunbear thus enables the projection of single-cell time-series snapshots to multimodal and multicondition views of cellular trajectories.

Figure 1.

Sunbear framework. (A) Sunbear takes as input a collection of measurements of single cells at multiple time points in two or more biological conditions (top) or data modalities (bottom). (B) During the training phase, Sunbear learns to decompose the original time-series profiles into four components: cell embedding, time point, batch, and condition. Batch and condition factors are represented by one-hot encodings. The time factor is represented by a sinusoidal encoding. The cell embedding is learned from the original profile and is conditionally independent of the other factors. In the multimodal setting, cell identities are aligned between data modalities. (C) In the prediction phase, Sunbear concatenates the query cell's identity factor while varying other factors to impute the cell's profile across time and conditions. By sharing cell embeddings across modalities, Sunbear allows joint temporal modeling of multimodal profiles.

Figure 2.

Single-cell profile inference across time and conditions. (A) Sunbear is trained on scRNA-seq profiles of whole-mouse embryos collected at alternating sexes along developmental time points (Qiu et al. 2024). (B) Sunbear is validated in three scenarios. In each scenario, one data block is held out from the training, and Sunbear is used to predict the profile of the missing block based on cells in the query block (outlined in turquoise). Sunbear's prediction is compared against the baselines using the held-out block's nearest existing measurements with the desired sex factors. The pseudobulk Pearson's correlation per major cell trajectory is calculated between the held-out profile and predictions/baselines. (C) Cross-time evaluation: query and baseline are selected either from the closest previous time point (left) or from the closest subsequent time point (right). The pseudobulk Pearson's correlation between the original held-out profile and predicted (y-axis) and baselines (x-axis) is plotted for each major cell trajectory in each held-out time point. Each dot represents a cell trajectory per held-out time point, and numbers indicate the number of dots above and below the diagonal line. P-values are calculated by a one-sided Wilcoxon rank-sum test. (D) Cross-sex prediction: similar to that in C, except query cells are selected from the opposite sex to the held-out data. (E) Cross-sex prediction: similar to that in D, except we enforce a strict baseline model by taking the mean of the previous and subsequent time point per cell trajectory.

Figure 3.

Sex differences in mouse embryonic development. (A) Pairwise comparison of Sunbear prediction and the nearest neighbor baseline in recapitulating differential expression patterns in sex-matched time points. Each dot indicates the AUROC score of recapitulating female-/male-biased patterns in each sex-matched time point and cell type. (B) Similar to in A, pairwise comparison of Sunbear prediction and the nearest neighbor baseline in ranking escape genes to be more female-biased than all other genes on the X Chromosome. (C,D) Predicted temporal sex-biased log fold change in glutamatergic neurons and border-associated macrophages. Each line represents a gene that is predicted to be consistently higher in females than males and is colored by whether the gene is a known constitutive escape gene or not. (E) Distribution of predicted sex-biased scores of genes (0 = extremely male-biased, 1 = extremely female-biased), grouped and colored by whether the gene is upregulated (pink) or downregulated (blue) in Kdm6a KO versus WT samples in CD4⁺ cells. P-values are calculated by one-sided Wilcoxon rank-sum tests. (F) Gene Ontology biological processes enriched in consistently female- and male-biased genes in border-associated macrophages. Nonredundant terms with the smallest FDR are selected for visualization. No enrichment of biological processes is found in glutamatergic neurons.

Figure 4.

Multimodal temporal inference. (A) A UMAP embedding suggests that scRNA-seq and scATAC-seq profiles are well aligned across time and batch. Only time points with both scRNA-seq and scATAC-seq available are shown. (B) AUROC of the predicted differential accessibility pattern relative to those derived from the original data sets. AUROC is calculated per cell type, and differential accessibility is calculated between each held-out time point and each query time point (shown as “query time point → held-out time point”). (C) Peak-wise AUROC of scATAC-seq profiles predicted based on scRNA-seq relative to the original scATAC-seq profile in each held-out time point. AUROCs are calculated across all cells. (D) Workflow for calculating the dynamic association between peaks and genes. A query cell's scRNA-seq profile is fed into Sunbear to predict temporal patterns of gene expression and chromatin accessibility. For each pair of chromatin region and its proximal gene, we calculate the correlation coefficient between them with incremental time shifts, which results in a TLCC vector (column). (E) Visualization of predicted peak region–gene relationships. Heatmap of TLCC matrices on randomly selected 5000 peak regions with accessibility changes ahead of (“before”) or subsequent to (“after”) nearby gene expression. Peak regions are sorted based on the time shift with the maximum TLCC.