Cell types’ locations and relative proportions are consistent across morphologically normal gastruloids.

a. Representative gastruloid with the expression of T shown in magenta. Each spot is one transcript (count), nuclei are shown in gray. Thick dashed line denotes the AP axis, thin lines outline the anterior and posterior halves of the gastruloid. Scale bar is 200 μm. b. Gallery of representative gastruloids. Scale bar is the same as a. c. Proportion of each cell type across all gastruloids, ordered by median proportion. d. Coefficient of variation for distributions shown in c. e. Correlation between the ratio of anterior area to posterior area across gastruloids versus the proportion of specific cell types. Cell types shown had a statistically significant (p < 0.05) correlation; all cell types are shown in Figure S1.4e. f. Smoothed density of each cell type along the AP axis. Nuclei in each gastruloid were projected onto the AP axis, and their position was normalized to the total length (left) or symmetrically around a midpoint defined by T expression (right). Traces from individual gastruloids are shown by individual curves.

Pairwise cell type interactions quantify cell type mixing and interaction motifs.

a. Distribution of mixing indices across gastruloids. b. Exposure index, which represents how frequently the source cell type (row) is found next to the target cell type (column). The maximum possible value is 1, self-interactions excluded for readability. c. Significantly enriched and depleted triplet combinations of cell types, ordered by effect size from left to right. Bars are colored by the composition of the triplet they represent. Only triplets found in over half of the gastruloids are shown.

Progressive NMP differentiation is revealed by the single cell L-score.

a. Illustration of NMP differentiation with predictions for where along the trajectory gene expression is shared or unique. b. Per cell expression of all NMP exclusive genes summed versus all other gene categories summed, one plot per type. n=25514 cells typed as NMP, presomitic mesoderm, or spinal cord in n=26 gastruloids. Pearson correlation is shown with colored lines. c. Expression plot showing the relative amount of Eogt vs Pax6 (left) and Eogt vs Rfx4 (right) per cell in a representative gastruloid. Pax6 and Rfx4 are both annotated as spinal cord-associated genes. d. Per cell expression of all spinal cord exclusive genes summed versus expression of all presomitic mesoderm genes summed. n=25514 cells typed as NMP, presomitic mesoderm, or spinal cord in n=26 gastruloids. e. Per cell expression scatterplots of the two pairs of genes shown in b). The y axis of each is the per-cell expression of Eogt. The x-axis is the per-cell expression of Pax6 (left) or Rfx4 (right). R is Pearson’s r, scL is scL-score. Count data is shown in black; smoothed 2D densities are shown in orange. f. Distribution of scL-score values for pairs of spinal cord and presomitic mesoderm genes (left), NMP and presomitic mesoderm genes (center), and NMP and spinal cord genes (right). g. Hierarchical clustering heatmap of scL-score vectors for NMP, presomitic mesoderm, spinal cord, and combined category genes. scL-score difference vectors were used for clustering; scL-score values are shown in the heatmap. Blocks highlighted are discussed in the text. h. Tree of the hierarchical relationships between genes resulting from the clustering shown in e). Color indicating the gene type is shown at the bottom, legend is the same as in g. The density plots shown are the summed, averaged gene expression of all genes in the leaves up until that node, smoothed with a 2D density kernel estimate (see Methods for details about smoothing).

Clustering L-score vectors clearly resolves cell types and reveals novel genetic interactions.

a. Heatmap of all scL-score values for all genes in the panel (excluding poorly detected and cell cycle genes, see Methods, n=166). Colored bars on the top and right hand sides indicate if a gene is associated with a particular cell type. Heatmap was hierarchically clustered by row. b. Expansion of a presomitic mesoderm cluster (i.). Clustering relationships are indicated with the dendrogram, and summed densities for all genes in an example gastruloid show where the genes in each cluster are expressed spatially. Clustering relationships determined from the full gene set shown in a. c. Expansion of the posterior cell type cluster (ii.). Clustering relationships are indicated with the dendrogram, and summed densities for all genes in an example gastruloid show where the genes in each cluster are expressed spatially. Clustering relationships determined from the full gene set shown in a. d. Expansion of the endothelial cluster (iii.). Clustering relationships are indicated with the dendrogram, and summed densities for all genes in an example gastruloid show where the genes in each cluster are expressed spatially. Clustering relationships determined from the full gene set shown in a. e. Expression of two example pairs of genes from the endothelial cluster: Cldn5 and Tgfb1, and Cldn5 and Sox17.

Spatial L-score reveals tissue-level patterns of gene expression.

a. Illustration of how the density-based L-score (spatial L-score) is calculated, and how the value changes as a function of how much the density estimate is smoothed. b. Clustered heatmap of the spatial L-score for a representative gastruloid. Purple box indicates a mixed cluster of endothelial and endoderm genes. c. An image of the gastruloid used to generate the heatmap in b. Cell types are indicated by color, legend is the same as in b.

Endothelial precursors show unique organization and distinct, spatially-dependent cell states.

a. Representative gastruloid showing two distinct morphologies of endothelial cells: somite-associated (anterior) and endoderm-associated (posterior). b. Exposure index for endothelial cells compared to all other cell types. Dataset-wide average (n=26 gastruloids) shown in the colored bars, grey lines indicate the values for the gastruloid shown in a. c. Exposure index for endoderm cells compared to all other cell types. Dataset-wide average (n=26 gastruloids) shown in the colored bars, grey lines indicate the values for the gastruloid shown in a. d. Genes differentially expressed in somite-associated or endoderm-associated endothelial cells. Bar color represents the cell type associated with the gene (if any); genes that were not previously known to be associated with a cell type in gastruloids are shown in gray. Bars are ordered by significance (adjusted p value) from greatest (Pecam1, adjusted p=8.49e-22) to least (Top2a, adjusted p=7.65e-3). e. Spatial distribution of expression for example endoderm-associated genes shown in d). Each plot shows one gene, each dot is a single transcript, and cells typed as endothelial are outlined in black. f. Spatial distribution of expression for example somite-associated genes shown in d). Each plot shows one gene, each dot is a single transcript, and cells typed as endothelial are outlined in black.

Gastruloid shape analysis and representative raw seqFISH images, spot detection, and gene deconvolution.

a. Proportion of each plate of gastruloids that, at 120 hours, were scored as ‘correct’. The correct phenotype was elongated with one axis and a clear anterior and posterior domain. The experiments used to generate the samples used in this paper are shown with colored dots, with colors indicating the day. One of the datasets used multiple plates pooled together, but the plates are shown individually for clarity. b. PCA projection of the morphological characteristics of all gastruloids generated in the 3 separate experiments analyzed in this study. n = 529 total gastruloids, 84-90 gastruloids per plate. Highlighted points are the extreme values in each dimension; images of the relevant gastruloid are shown in insets. c. Area and elongation of all gastruloids generated for this study (n=529), colored by date. d. Representative images showing nuclear staining, spot detection by hybridization (hyb) and deconvolution of Dll1.

Comparison of AP-axis gene expression centers of mass between individual gastruloids reported in this work and previously-reported tomoSeq performed on gastruloids.

a. Center of mass of gene expression in all gastruloids compared to the center of mass of genes expression in one gastruloid analyzed by tomo-seq as reported in (van den Brink et al., 2020). The number of genes in each plot is shown in the title; numbers vary by gastruloid as some of the genes in the seqFISH panel were not consistently detected across gastruloids (see Methods). Orange points are the Hox genes included in the panel. b. Summary of the Pearson correlation of all samples shown in a) for all genes (blue) and Hox genes in the seqFISH panel (orange).

Gallery of individual gastruloids colored by cell type.

a. Gastruloids from the experiment conducted on 05/07/2025; color of each nucleus indicates cell type (legend on bottom left). b. Gastruloids from the experiment conducted on 09/24/2024; color of each nucleus indicates cell type (legend on bottom left). c. Gastruloids from the experiment conducted on 02/09/2024; color of each nucleus indicates cell type (legend on bottom left).

Cell type entropy scores, cell type proportion covariation, and correlation between individual cell type proportions and morphological characteristics.

a. Cell type score entropy distributions for each cell type. Each violin shows all the values for that cell type across all gastruloids. n=79607 nuclei across 26 gastruloids. b. Per-gastruloid cell type proportions (summarized in Main Figure b). c. Per-gastruloid cell type proportions including untyped cells (’None’). d. Co-variation in proportion (normalized by Centered Log-Ratio (CLR) transformation) between cell types. Green boxes indicate significant hits (adjusted p value < 0.05 at a 5% FDR). e. Correlation between the per-gastruloid mixing index with the proportion of the labeled cell type in that gastruloid. Pearson’s r and the significance (p) for each relationship is shown in the black box. ** = p < 0.01, * = p < 0.05.

Analysis of the AP-axis location and rank order of cell types.

a. Bootstrapped null distributions for the difference in peak location to the overall mean for each cell type along the AP axis length normalized (top row) or normalized to length with the midpoint determined by T expression (bottom row). Red lines show the average value for the observed peak differences in the 26 samples shown in Figure S1.3a,b,c. The fold change from the mean is shown in the inset, as is the p value for the difference. b. Comparison between the fold changes reported in the top and bottom row of a). c. Rank heatmap for the peak location of each cell type for each gastruloid (length normalized). d. Rate of rank swapping between cell types; 0=conserved, 0.5=equal likelihood of occurring at two adjacent ranks.

Exposure indices triplet motifs.

a. Exposure indices for all cell types shown individually as bar plots. b. Significantly varying (red) and conserved (blue) neighbor interactions across gastruloids (n=26). c. Exposure index matrix with self-interactions included.

Cell type scoring robustness.

a. Entropy cutoff used to generate a subset of high-confidence cells. At this cutoff (entropy <= 1.5) 99.7% of ‘none’ type cells are removed, as well as ∼50% of cardiac mesoderm and paraxial mesoderm cells. b. Correlation between the overall mixing index per-gastruloid, calculated before and after entropy filtering (n=26 gastruloids). c. Correlation between the proportion of each cell type and the change in mixing index before and after entropy filtering. Pearson’s r and the associated p value are shown in the inset. d. Full exposure index matrix recalculated after entropy filtering.

Single cell L-score of NMP, presomitic mesoderm, and spinal cord genes.

a. Left: cell type scores for all nuclei in the gastruloid used in this figure. Right: cell type scores for all nuclei typed as NMP, presomitic mesoderm, or spinal cord in each of those categories. b. Spatial distribution of the expression of Nkx1-2 and Rfx4 in an example gastruloid. c. Nkx1-2 transcripts in the same gastruloid as in c. d. Distribution of scL-score values between pairs of terminal cell type genes and genes annotated as belonging to that cell type and NMP genes.

Illustration of L-score calculation.

a. Illustration of the calculation of the scL-score.

Comparison between scL-score, Exclusive Expression Index (EEI), and coefficient of coexpression (COEX) calculated from simulated data representing common expression scenarios.

a. Simulations of several common gene expression scenarios. The top row shows the simulated per-cell expression distributions for each scenario, and the bottom row shows the expression values rank-ordered by the expression of gene 1. This rank-ordering is used to calculate the L-score. b. scL-score, EEI, and COEX (from the COexpression Table ANalysis (COTAN) framework) for the scenarios shown in a. We compared the scL-score to two existing measures of exclusivity. The Exclusively Expressed Index (EEI) (Nakajima et al., 2021) is bounded below by 0 and increases with mutual exclusivity. EEI is computed from binary zero/non-zero quantification and is designed specifically to measure exclusivity but not coexpression. The coefficient of coexpression (COEX) from the COexpression Table ANalysis (COTAN) framework (Galfrè et al., 2021) can also be used to quantify relationships between genes: positive values indicate coexpression, negative values indicate mutual exclusivity, and values near 0 indicate little structured relationship. In the two mutually exclusive simulations, all three methods detected exclusivity. In the three simulations of genes independently expressed in all cells (but with varying relative expression levels), EEI and COEX were 0, while the scL-score remained close to 0 (0.102, -0.225, and -0.025, respectively), consistent with little structured relationship. In the weak coexpression simulation, the scL-score was positive (0.770), EEI remained close to 0, and COEX was positive (0.340), indicating detectable but modest coexpression. In the perfect coexpression simulation, the scL-score reached 1.000, EEI was 0, and COEX was strongly positive (1.000). Together, these simulations show that all three methods detect strong mutual exclusivity, and both scL-score and COEX distinguish positive coexpression from exclusivity and from unstructured expression.

Comparison between scL-score, Exclusive Expression Index (EEI), and coefficient of coexpression (COEX) for 4 example gene pairs.

a. Expression of 4 example gene pairs in a single gastruloid (top) and their calculated scL-score, EEI, and COEX values (bottom). We calculated the scL-score, EEI, and COEX for four gene pairs from one gastruloid sample (2025-05-07_roi2). The scL-score consistently delineated gene pairs possessing opposing expression profiles, while EEI was not always able to measure those exclusivity patterns (Pax6-Eogt: scL-score=-0.572 and -0.526, EEI=0; Rfx4-Eogt: scL-score=-0.970 and -0.955, EEI=0.171; Nkx1-2-Rfx4: scL-score=-0.537 and -0.615, EEI∼0). Only the scL-score was able to detect the coexpression pattern present between the positively associated expression profiles of Cdx4 and Cdx2 (scL-score=0.413, EEI=0, COEX -0.331) (shown visually in Figure S3.4a). Thus, while EEI was informative for measuring gene expression relationships characterized by mutual exclusivity, the scL-score more clearly separated positive, random, and mutually exclusive relationships on a single signed bounded scale. The COEX value trended in the opposite direction than expected, but this may be due to the fact that it cannot be calculated on single gene pairs and necessarily uses information from the entire count table, which here only consisted of 6 genes. These comparisons combined with the simulations in Figure S3.3, clarify a conceptual advantage of the L-score over the EEI and COTAN frameworks. In contrast, by leveraging the ranked structure of transcript counts across cells, the L-score framework does not binarize expression and does not require fitting a parametric distribution. It can be calculated on single gene pairs, and is more sensitive to mutual exclusivity.

Refining and assessing marker gene panels with scL-score values.

a. A representative gastruloid with cell types as scored by the full marker gene panel (left) and entropy values of the score distributions of each cell (right). b. The same gastruloid scored with a reduced panel of marker genes filtered by the properties of their scL-score values with all other genes. c. Heatmap showing the average scL-score values between sets of marker genes from the hand-curated panel used in this paper. Panels were filtered by average per-cell expression, see Methods for details.

Using the scL-score to assess unsupervised clustering and identify marker genes.

a. Left: UMAP of all gastruloid nuclei from the dataset taken on 05/07/2025 with > 40 counts projected onto a UMAP and colored according to Leiden clusters. Right: Cluster assignment for a representative gastruloid projected into spatial coordinate. b. Heatmap showing the average scL-score values for the top 10 genes associated with each cluster shown in a. We reasoned that since the L-score quantifies mutually-exclusive gene expression, a property of good marker genes, it might be possible to use a gene’s scL-score values to refine our marker gene panel. We calculated the scL-score vectors, and then filtered for only genes which had at least 10% of their L-score values < -0.8, and at least 12% > -0.3. This did shift the balance of differentiation front and presomitic mesoderm genes, however it slightly increased the entropy scores (Figure S3.5a,b), indicating that for small, hand-chosen marker gene panels, removing genes decreases scoring confidence, albeit only slightly. However, we also used the average scL-score of genes associated with one type compared to another type to quantify type similarity (in terms of per-cell gene expression), and found marker pairs that distinguished cell types or were concordantly expressed in the same cell types (Figure S3.5c). This method can also be used to post-facto analyze similarity of clusters produced by Leiden or other clustering methods. To demonstrate this we pseudo-bulked all the nuclei for our entire dataset and clustered using a standard scanpy workflow (Figure S3.6a). We then took the top genes for each cluster, and calculated the average scL-score between the groups (Figure S3.6b). This analysis demonstrated that the clusters produced by unsupervised clustering, while they could be qualitatively mapped onto the cell types we expected to see (Figure S3.5a), were much more similar in terms of gene expression. Just using the top genes from this analysis would not be sufficient to produce marker genes. However, we used scL-score values to find genes that were divergently expressed, even among similar clusters (Figure S3.6b).

Heatmap of L-score values including cell cycle genes and NMP subcluster.

a. Heatmap of all scL-score values for all genes in the panel (excluding poorly detected but including cell cycle genes n=202 genes). Colored bars on the top and right hand sides indicate if a gene is associated with a particular cell type. Heatmap was hierarchically clustered by row. b. Highlight of the NMP cluster from Figure 4a. Clustering relationships are indicated with the dendrogram, and summed densities for all genes in an example gastruloid show where the genes in each cluster are expressed spatially. Clustering relationships determined from the full gene set shown in Figure 4a.

Quantification of cell type clustering by scL-score analysis and robustness of clustering to number of genes included.

a. Average intra-cell type mean pairwise branch length (cell type dispersion) for the tree produced by scL-score clustering on 166 well-detected genes (excluding cell cycle genes) for real (red line) and permuted (gray distribution) leaf identities. b. Cell type dispersion of the tree produced by scL-score clustering on 202 well-detected genes (including cell cycle genes) for real (red line) and permuted (gray distribution) leaf identities. c. Cell type dispersion of observed (red) and permuted (gray) leaf identities as a function of the number of genes used to produce the underlying tree. d. Gap between the mean permuted cell type dispersion and the real cell type dispersion as a function of the number of genes used to produce the underlying tree.

Clustering on transformed scL-scores.

a. scL score values for all pairs of well-detected non-cell cycle genes (n=166). Row order is determined by hierarchical clustering on transformed scL-score values: (1-scL)/2. Values in the heatmap on the left are the untransformed values, the dendrogram on the right shows the clustering relationships. Column order was set to be the same as the row order. b. scL score values for all pairs of well-detected genes (n=202). Row order is determined by hierarchical clustering on transformed scL-score values: (1-scL)/2. Values in the heatmap on the left are the untransformed values, the dendrogram on the right shows the clustering relationships. Column order was set to be the same as the row order.

Comparison of scL-score clusters to cNMF clusters.

a. cNMF stability analysis for varying numbers of components. cNMF was run with standard parameters on all nuclei from the dataset taken on 05/07/2025. b. Top 24 genes in each cluster averaged in space and summed together for visualization. 24 genes were chosen as this was the average number of genes per cluster when the L-score tree was truncated to produce 7 clusters for comparison (see c) c. Truncation of the scL-score tree to produce 7 clusters. The genes in each cluster were then visualized as in b. d. Qualitative clusters on the scL-score tree that have many genes associated with a single cell type. e. Quantification of the overlap between scL-score clusters and cNMF clusters. Overlap between each scL-score cluster and each cNMF cluster was quantified by Jaccard Index or Adjusted Rand Index. The scL-score clusters were then shuffled and the same values were quantified. The value of the overlap between each scL-score cluster (blue) or permuted scL-score cluster (gray) and the most similar cNMF cluster is shown. We ran cNMF (Kotliar et al., 2019) on our data to identify gene programs. Stability analysis indicated that 7 was the optimal number of components. To compare how the genes identified by cNMF compared to those which are related by scL-score similarity, we truncated the L-score tree to produce 7 clusters. We visualized the genes for each method in space. There was clear correlation between the clusters produced by both methods. The cNMF clusters looked ’cleaner’ (because the top 24 genes in each program were plotted, and thus many genes that did not meet this threshold in any program weren’t visualized at all), while the L-score clusters had a hierarchy of relatedness which was absent from the cNMF programs. Furthermore, we could select tree nodes qualitatively as being enriched for cell type marker genes, and found that these clusters also closely matched the cNMF clusters.

scL-score analysis of scRNA-seq gastruloid dataset.

a. Dendrogram produced by hierarchical clustering the scl-scores of the 207 genes in the panel in the work as measured by sc-RNAseq in (van den Brink et al., 2020): 4 samples of 120 hour gastruloids were pooled; after filtering for read quality 14304 cells remained. scL-score difference vectors were used as the distance measure for clustering as in Figure 4. b. Cell type dispersion (average pairwise branch length per cell type) of the real (red line) and permuted (gray distribution) tree shown in a. c. scL score heatmap of values calculated from (van den Brink et al., 2020); the same quality control, filtering, and clustering method was performed as is shown in a), but all well-detected genes (19075) were used to perform clustering. The scL-score value, rather than the distance, is shown for each gene pair. d. Example cluster found on the diagonal of the heatmap shown in c. e. Example cluster found on the diagonal of the heatmap shown in c. f. Example cluster found on the diagonal of the heatmap shown in c.

Spatial L-score heatmap of values averaged across all samples and direct comparison of scL-score and spatial L-score values for all gene pairs.

a. Heatmap of clustered spatial L-score values for all genes averaged across all gastruloids. n=202 genes. b. Comparison of scL-score and spatial L-score values for cell type markers. Colored dots represent intra-type pairs, with the color specifying the cell type. The gray trace is the smoothed density estimate for all pairs, including inter-type pairs and pairs including cell cycle or non-marker genes.

Clustered heatmap of the scL-score for all genes in an example gastruloid.

a. Heatmap of clustered scL-score values for all genes for the representative gastruloid shown in the main figure. Purple rectangle highlights endoderm and endothelial gene clusters.

Estimation of the effect of spot mis-assignment.

a. Cell type entropy score distributions as a function of nuclear dilation (higher = more pixels added to the original segmentation). Nuclei from the gastruloid in Figure 6a are shown. b. Expanded images from Figure 6e, f.

Expression of differentially expressed endothelial genes in a gastruloid with only one endothelial population.

a. Expanded images showing the expression of the same genes highlighted in Figure 6e, f in the gastruloid shown in Figure S1.3a iii.

Average nearest neighbor distances between cell centroids in each gastruloid in our dataset from 05/07/2025 (n=18 gastruloids).

These values multiplied by 2 were used as bandwidths for individual two-dimensional Gaussian kernels fitted to each RNA spot to produce a smoothed spatial profile of expression for a given gene on a particular gastruloid.

Marker genes for all scored cell types.

Genes that mark more than one cell type are shown with both.