we computed pseudobulk samples and carried out the de analysis using the same procedure and cut - offs as in the protocol - based de analysis. to estimate the transcriptomic similarity between neurons in hnoca and the human developing brain atlas, we first summarized the average expression of each neural cell type in the primary reference, as well as in each dataset of hnoca. for each hnoca dataset, only neural cell types with at least 20 cells were considered. highly variable genes were identified across the neural cell types in the primary reference using a chi - squared test - based variance ratio test on the generalized linear model with gamma distribution ( identity link ), given coefficient of variance of transcript counts across neural cell types as the response and the reciprocal of average transcript count across neural cell types as the independent variable. genes with benjamini – hochberg adjusted p values less than 0. 01 were considered as highly variable genes. similarity between one neural cell type in the primary atlas and its counterpart in each hnoca dataset was then calculated as the spearman correlation coefficient across the identified highly variable genes. to estimate the similarity of the core transcriptomic identity, which is defined by the coexpression of transcription factors, the highly variable genes were subset to only transcription factorsfor calculating spearman correlations. the list of transcription factors was retrieved from the animaltfdb v. 4. 0 database. to identify metabolically stressed cells in the datasets, we used the scanpy score _ genes function with default parameters to score the ‘ canonical glyco lysis ’ gene set obtained from the enrichr go _ biological _ process _ 2021 database across all neuronal cells from hnoca and refs.. to estimate the significance of the difference between the correlation of glyco lysis scores and whole transcriptomic similarities, and the correlation of glyco lysis scores and core transcriptomic identity similarities, we generated 100 subsets of highly variable genes, each with the same size as the highly variable transcription factor. transcriptomic similarities were calculated on the basis of those subsets, and then correlated with the glyco lysis scores. to characterize heterogeneity of telencephal ic npcs and neurons in hnoca, we first transferred the cell type labels ( as indicated as the ‘ type ’ label in the given metadata ) from the human neocortical development atlas to the hnoca telencephal ic npcs, intermediate progenitor cells and neurons, on the basis of transcriptomic correlation. in brief, each primary atlas cluster we obtained as mentioned above was assigned to a cell type as the most abundant cell type among cells in the cluster. the label of the best - correlated primary cluster was then transferred to every query cell. given the transferred label, together with the level 2 cell type annotation shown in fig. 1c, as the annotation label, scpoli was applied to the telencephal ic subset of hnoca for data integration. to benchmark how well different integration strategies recover the neuron subcell type heterogeneity, we generated four different clustering labels: ( 1 ) louvain clustering ( resolution, 2 ) with the original scpoli latent representation; ( 2 ) louvain clustering ( resolution, 2 ) with the updated scpoli representation; ( 3 ) louvain clustering ( resolution, 2 ) with pca of hnoca telencephal ic subset ( based on scaled expression of 3, 000 highly variable genes of the telencephal ic subset with flavor = ‘ seurat ’ ) and ( 4 ) louvain clustering ( resolution, 1 ) for each sample separately ( each with 3, 000 highly variable genes identified with flavor = ‘ seurat ’, followed by data scale and pca ). next, for each sample with at least 500 dorsal telencephal ic neurons, the adjusted mutual information scores were calculated between each of those four clustering labels with the transferred cell type label mentioned above as the gold standard, across the dorsal telencephal ic neurons as annotate d as the level 2 annotation. to create a comprehensive primary atlas of dorsal telencephal ic neurons for de analysis between neural organoids and primary tissues, we subset dorsal telencephal ic neurons or neocortical neurons from four different primary atlases. cells in five author - defined clusters ( 60, 57, 79, 45, 65 ) with high expression of map2, dcx and neurod6 were selected. cells with the following ‘ clusterv2 - final ’ labels were selected: ‘ neuron _ 28 ’, ‘ neuron _ 34 ’, ‘ gw19 _ 2 _ 29neuronneuron ’, ‘ neuron _ 30 ’, ‘ neuron _ 66neuron ’, ‘ gw18 _ 2 _ 42neuronneuron ’, ‘ neuron _ 33 ’, ‘ neuron _ 39neuron ’, ‘ neuron _ 35 ’, ‘ neuron _ 63neuron ’, ‘ neuron _ 9 ’, ‘ neuron _ 11 ’, ‘ neuron _ 20 ’, ‘ neuron _ 22 ’, ‘ neuron _ 5neuron ’, ‘ neuron _ 21 ’, ‘ neuron _ 18 ’, ‘ neuron _ 101neuron ’, ‘ neuron _ 17 ’, ‘ neuron _ 19 ’, ‘ neuron _ 16 ’, ‘ neuron _ 50neuron ’, ‘ neuron _ 12 ’, ‘ neuron _ 13 ’, ‘ neuron _ 68neuron ’, ‘ neuron _ 100neuron ’, ‘ neuron _ 25 ’, ‘ neuron _ 27 ’, ‘ neuron _ 53neuron ’, ‘ neuron _ 23 ’, ‘ neuron _ 26 ’, ‘ neuron _ 24 ’, ‘ neuron _ 102neuron ’, ‘ neuron _ 72neuron ’, ‘ neuron _ 15 ’, ‘ neuron _ 29 ’ and ‘ neuron _ 35neuron ’ on the basis of their high expression of neurod6 and foxg1. cells dissected from dorsal telencephalon that were annotate d as neurons with and only with the vglut ntt label were selected. cells annotate d as excitatory neurons were selected. the curate d clusters of the wang et al. primary atlas, as described earlier, were also subset to those with excitatory neuron labels. the selected dorsal telencephal ic neuron subsets of the atlases were merged into the joint neocortical neuron atlas. next, cells in the joint neocortical neuron atlas were correlated with the average expression profile of each excitatory neuron cluster of the wang et al. atlas. the cluster label of the best - correlated cluster was assigned to each cell in the joined neocortical neuron atlas, so that cell cluster labels were harmonized for all cells in the atlas. label - aware data integration was then performed using scpoli. on the basis of the scpoli latent representation, louvain clustering was performed on the joint neocortical neuron atlas ( resolution, 1 ). this cluster label was transferred to the dorsal telencephal ic neurons in hnoca with max - correlation manner across highly variable genes defined on average transcriptomic profiles of clusters in the joint neocortical neuron atlas. we used scarches to map scrna - seq data from the neural organoid morphogen screen to both the scanvi model of the human developing brain atlas and the scpoli model of the hnoca. in both cases, the ‘ dataset ’ field of the screen data was used as the batch covariate, which indicates belonging to one of the three categories: ‘ organoid screen ’, ‘ secondary organoid screen ’ or ‘ fetal striatum 21pcw ’. for mapping to the primary reference, we used the scvi - tools implementation of scarches without the use of cell type annotations and trained the model for 500 epochs with weight _ decay of 0 and otherwise default parameters. for mapping to hnoca we used scarches through scpoli and trained the model for 500 epochs without unlabelled prototype training. we included 11 scrna - seq datasets of neural organoids, which were designed to model 10 different neural diseases including microcephaly, amyot rophic lateral sclerosis, alzheimer ’ s disease, autism, fxs, schizophrenia, neuronal heterotopia, pitt – hopkins syndrome, myotonic dystrophy and glioblastoma. count matrices and metadata were directly downloaded for the ten datasets with processed data provided in the gene expression omnibus or arrayexpress.