We identified the alpha and beta ER stress clusters in Figure 5 by performing downstream analysis, specified in this section, on the integrated joint embedding produced by Harmony. After Harmony integration, we performed clustering analysis to find novel subtypes. Clustering was done on the trimmed shared nearest neighbor graph with the Louvain algorithm45 , as implemented in the Seurat package BuildSNN and RunModularityClustering functions. We used parameters resolution=0.8, k=30, and nn.eps=0. We identified several clusters within the alpha, beta, and ductal cell populations. For each cluster, we performed differential expression analysis within the defined cell type. That is, we compared alpha clusters to all other alpha cells. For differential expression, we used the R Limma package18 (link) on the normalized data. We included technology and library complexity (log number of unique genes) as covariates in the linear models. We used the top 100 over-expressed genes for each cluster, weighted by the t-statistic, to perform pathway enrichment with the enrichR46 ,47 R package, using the three Gene Ontology genesets48 ,49 . The ductal subpopulation was enriched for ribosomal genes; we did not follow up on this cluster.