a random - walk - with - restart procedure was then applied to smooth the raw scores across the knn graph of the primary atlas. in brief, we first represented the primary atlas knn graph as its adjacency matrix ( a ), followed by row normalization to convert it into a transition probability matrix ( p ). with the raw scores represented as a vector s 0, in each iteration t, we generated s t as [ 12pt ] $ $ _ = } } _ + ( 1 - ) ^ _ $ $ s t = α s 0 + ( 1 − α ) p t s t − 1 source paper: pmc11578878 this procedure was performed 100 times to get the smooth presence scores that were subsequently log transformed. scores lower than the 5th percentile or higher than the 95th percentile were trimmed. the trimmed scores were normalized into the range of as the final presence scores in the hnoca dataset. given the final presence scores in each of the hnoca datasets, the max presence scores in the whole hnoca data were then easily calculated as the maximum of all the presence scores for each cell in the primary atlas. a large ( close to one ) max presence score indicates a high frequency of appearance for the cell type or state in at least one hnoca dataset whereas a small ( close to zero ) max presence score suggests under - representation in all the hnoca datasets. to test the cell type compositional changes on admission of certain morphogens from different organoid differentiation protocols, we used the pertpy implementation of the sccoda algorithm. sccoda is a bayesian model for detecting compositional changes in scrna - seq data. for this, we have extracted the information about the added morphogens from each differentiation protocol and grouped them into 15 broad molecule groups on the basis of their role in neural differentiation ( supplementary table 1 ). these molecule groups were used as a covariate in the model. the region labels transferred from the primary atlas were used as labels in the analysis ( cell _ type _ identifier ). for cell types without regional identity, the cell type labels presented in fig. 1c were used. pluripotent stem cells and neuroepithelium cells were removed from the analysis because they are mainly present in the early organoid stages. we used bio _ sample as the sample _ identifier. we ran sccoda sequentially with default parameters, using no - u - turn sampling ( run _ nuts function ) and selecting each cell type once as a reference. we used a majority vote - based system to find the cell types that were credibly changing in more than half of the iterations. to complement the composition analysis conducted with sccoda, we devised an alternative approach to test for differential composition using regularized linear regression. we fit a generalized linear model with the region composition matrix as the response y and molecule usage as independent variables x: [ 12pt ] $ $ y x } $ $ y ~ x β source paper: pmc11578878 the model was fit with lasso regularization ( alpha = 1 ) using gaussian noise and an identity link function. the regularization parameter lambda was automatically determined through cross - validation as implemented in the function cv. glmnet ( ) from the glmnet r package. all non - zero coefficients β were considered as indications of enrichment and depletion. to study the transcriptomic differences between organoid and primary cells, we subset hnoca using the final level 1 annotation to cells labelled ‘ neuron ’. we furthermore subset the human developing brain atlas to cells that had been assigned a valid label in the neuron _ ntt _ label annotation column. we added an extra two datasets of fetal cortical cells from ref. for the data from ref., we subset the data to cells labelled ‘ fetal ’ and estimated transcripts per million reads for each gene in each cell using rsem given the star mapping results. we then computed a pca, a knn graph, umap and leiden clustering ( resolution 0. 2 ) using scanpy. we then selected the cluster with the highest stmn2 and neurod6 expression as the cortical neuron cluster and used only those cells. for the data from ref. we subset the datasets to cells annotate d as ‘ neuronal ’ in supplementary table 5 ( ‘ cortex annotations ’ ) of their publication and computed a pca, neighbourhood graph and umap to visualize the dataset. we found that only samples from the individuals cs14 _ 3, cs20, cs22 and cs20 contained detectable expression of stmn2 and neurod6 so we subset the dataset further to only cells from those individuals. to compute de between hnoca cells and their primary counterparts, we first aggregated cells of the same regional neural cell type into pseudobulk samples by summing the counts for every sample ( annotation columns, ‘ batch ’ for hnoca; ‘ sampleid ’ for the human developing brain atlas; ‘ sample ’ for ref. and ‘ individual ’ for ref. ) using the python implementation of decoupler ( v. 1. 4. 0 ) while discarding any samples with fewer than ten cells or 1, 000 total counts. we then subsetted the feature space to the intersection of features of all datasets and removed any cells with fewer than 200 genes expressed. we further removed any genes expressed in less than 1 % of neurons in hnoca and any genes located on the x and y chromosomes. out of the remaining 11, 636 genes, on average, 99 % were reported in each of the constituent hnoca datasets. for each regional neural cell type, we removed any sample from the pseudobulk data that was associated with an organoid differentiation assay with fewer than two total samples or fewer than 100 total cells. we next used edger to iteratively compute de genes between each organoid differentiation protocol and primary cells of the matching regional neural cell types for every regional neural cell type while correcting for organoid age in days, number of cells per pseudobulk sample, median and standard deviation of the number of detected genes per pseudobulk sample. we used the data from ref. ( the human developing brain atlas mentioned above ), ref. as primary data for the de comparison in the cell type ‘ dorsal telencephal ic neuron nt - vglut ’, whereas for all other cell types we used the human developing brain atlas as the fetal dataset. we used the edger genewise negative binomial generalized linear model with quasi - likelihood f - tests. we deemed a gene significantly de if its false - discovery rate ( benjamini – hochberg ) corrected p value was smaller than 0. 05 and it had an absolute log 2 - fold change above 0. 5. we used the gseapy python package to carry out functional enrichment analysis in our de results using the ‘ go _ biological _ process _ 2021 ’ gene set. to evaluate the effect of different primary datasets on the de results, we computed the de between dorsal telencephal ic neuron nt - vglut from the hnoca subset generated with the protocol from ref. and the matching cell type from the braun et al. primary dataset as well as the data from ref.. to prevent technology effects to affect this analysis, we only used cells generated with the 10x genomics 3 ′ v. 2 protocol in this comparison. we generate pseudobulk samples as described above and corrected organoid age in days and number of cells per pseudobulk sample in the de comparison. we used the same edger - based procedure and cut - offs as described above. we used the scipy fcluster method to cluster genes on the basis of their log - fold changes in the two primary datasets. we grouped clusters to represent consistently upregulated, consistently downregulated and three different inconsistently regulated groups of genes. we computed functional enrichment of each gene group as described above. to evaluate the effect of different organoid datasets on the protocol - based de analysis, we computed de between dorsal telencephal ic neuron nt - vglut of every organoid publication ( further split by protocol, where more than one protocol was used in a publication ) and the matching cell type in the dataset from ref..