a, Schematic of normal brain developmental trajectories queried in this study and their comparison to organoid models. Normal cortical development requires the emergence of a diversity of progenitor cell types from a seemingly uniform neuroepithelium. Through a sequence of cell-type specification and maturation, progenitor cells undergo neurogenesis and gliogenesis to generate the cellular diversity of the cortex. Areal identities are specified during this process and comprise a core property of developing neurons.

a, Cortical organoid protocols using different levels of directed differentiation were evaluated using scRNA-seq and immunohistochemistry. Stem cells were expanded on matrigel, dissociated to single cells, and re-aggregated in v-bottom low-adhesion plates. Small molecules were used to promote forebrain induction and after 18 days cells were moved to 6-well plates on an orbital shaker. Organoids were maintained in culture and collected from weeks 3 to 24. b, Protocol schematics for other methods used to differentiate whole brain and cortical organoids, which have published single-cell data. Publicly available data were used for comparative analyses with our collected data.

a, Organoids derived from the 13234 induced PSC line underwent the least and directed differentiation protocols, were collected at weeks 5 and 10, and were processed for immunohistochemistry. Organoids from both protocols were stained with SOX2 to mark progenitors, HOPX to identify outer radial glia and TBR2 to label intermediate progenitor cells. Cultures were also stained with CTIP2 to mark deep layer neurons and SATB2 to identify upper layer neurons. At week 5 all progenitor subtypes were present. and by week 10 both deep and upper layer neurons were detected. b, Organoids from the H28126 induced PSC line were differentiated using the least, directed and most directed protocols. All progenitor types marked by SOX2, HOPX and TBR2 and CTIP2⁺ and SATB2⁺ neuronal populations were present by week 5. Expression of all markers decreased by week 10. Organoid staining validation of broad cell types was repeated independently three times.

a, t-SNE plots depicting the single-cell analysis of primary cortical cells as coloured by cluster, age of sample and cortical area. Stacked histograms showing composition of each cluster for these metadata properties are also included. b, t-SNE plots depicting the single-cell analysis of cortical organoid cells as coloured by cluster, protocol, pluripotent stem cell line (induced PSC or human embryonic stem cell) and age of sample. Stacked histograms showing composition of each cluster for these metadata properties are also included. Source data

a, Re-analysis of published single-cell sequencing in organoid samples. t-SNE plot is coloured by cell-type designation, and the feature plots depict the same cell populations as presented in Fig. 1. b, t-SNE plots depicting the single-cell analysis of published organoid cells as coloured by cluster, protocol (including paper of origin) and FOXG1 expression. c, Recapitulation of the heat map in Fig. 2, using published organoid clusters from above and comparing to primary reference dataset from this paper. Quantification of correspondence shows the quantitative correlation from the best match in the heat map for each category of class, state, type and subtype, averaged across all clusters (primary: n = 189,409 cells from five individuals collected independently; published organoid data: n = 109,813 cells from 7 datasets collected independently by different scientific groups; two-sided Welch’s t-test evaluating mean + s.d.; subtype versus type, * P = 0.0193; subtype versus state, ***P = 0.00017). Source data

a, Re-analysis of published single-cell sequencing⁵, in which a reproducible cortical organoid protocol was presented. The t-SNE plot is coloured by cell-type designation, and the feature plots depict the same cell populations as presented in Fig. 1. b, Recapitulation of the heat map in Fig. 2, using published organoid clusters⁵ and comparing to the primary reference dataset from this paper. Quantification of correspondence shows the quantitative correlation from the best match in the heat map for each category of class, state, type and subtype, averaged across all clusters (primary: n = 189,409 cells from five individuals collected independently; organoid data⁵: n = 166,241 cells from an independently collected dataset; two-sided Welch’s test was used to evaluate mean + s.d.; subtype versus state ****P = 1.8 × 10⁻⁷). c, Pseudoage analysis of published organoids⁵ mirrors the organoids in this study with low correspondence between pseudoage and actual age. Pseudoage calculation is indicated by the graph line and shading represents the geometric density standard error of the regression. d, Area identity was assigned for all excitatory neurons from ref. ⁵ and each organoid consisted of heterogeneous areal identities, consistent with the observations in the organoids from this study. Source data

a, Composition of each organoid by cell-type designation. FOXG1 expression across all organoid samples is plotted by feature on the right. b, Comparison of organoid subtype as determined by this study versus three control analyses. Graphically, the column indicates subtype correspondence; error bar, s.d. The first analysis was performed by halving the primary dataset randomly and without overlap and then comparing the subclusters from the two datasets. This age- and method-matched analysis shows that primary variation is significantly lower than the variation between organoids and primary cells, as indicated by the significantly higher subtype correlation between primary datasets (organoids: n = 242,349 cells collected from 37 organoids from 4 biologically independent samples from 4 independent experiments; primary data: 189,409 cells from 5 biologically independent samples from 5 experiments; ****P = 2.0 × 10⁻²⁴, two-sided Welch’s t-test). A similar analysis was performed comparing the primary data from this study to data collected by microfluidic approaches¹⁹. Although the ages, capture method and number of cells varied greatly, subtype correlation between the published primary data and the data in this study is significantly higher than the subtype similarity between organoids and primary samples¹⁹ (n = 4,261 cells from 48 biologically independent samples across more than 35 independent experiments, ****P = 2.0 × 10⁻⁵). We additionally performed this analysis between two published datasets for cells from adult humans, comparing middle temporal gyrus⁴² (MTG, n = 15,928 cells) from an older adult with distinct brain regions from young adults in the control samples of a study on autism spectrum disorder⁴³ (ASD, n = 104,559). Despite differences across ages and individuals, who could be expected to have unique cortical gene-expression profiles based upon sensory experience, the distinct cortical regions isolated and the different capture methods, the subtype correlation between these two primary datasets is significantly higher than the correlation between organoid cells and primary cells (**P = 0.0076). c, Subtype correlation as calculated and shown in Fig. 2, broken down by protocol and pluripotent line, in which bars indicate subtype correlation and error bars show s.d. The least directed protocol was significantly better at recapitulating cell subtype than the most directed protocol (*P = 0.0483, two-sided Welch’s t-test), consistent with recent findings⁵. We also observed that the induced PSC line 1323_4 generated significantly more similar subtypes to primary samples than WTC10 or H1 (**P = 0.0013 and 0.0089, respectively). d, Clustering and subtype analysis was performed between all organoids and primary PFC samples and primary V1 individually. Subtype correlation did not change regardless of the area to which organoids were compared. ‘Overall’ refers to the subtype correlation observed when comparing all organoids cells to all primary cells and is shown for comparison. Histogram bars show subtype correlation and error bars show s.d. (n = 242,349 cells from 37 organoids across 4 independent experiments). e, Subtype correlation analysis was performed across all organoid stages (n = week 3: 38,417 cells, week 5: 26,787 cells, week 8: 11,023 cells, week 10: 50,550 cells, week 15: 2,722 cells, week 24: 4,506 cells from 4 independent experiments) and all primary ages (n = GW6: 5,970 cells, GW10: 7,194 cells, GW14: 14,435 cells, GW18: 78,157 cells, GW22: 83,653 cells from 5 independent experiments). Histogram bars show subtype correlation and error bars show s.d. Week-3 organoids are more similar to younger primary stages, and week-15 organoids are most similar to older primary ages. Other ages correspond similarly well to the primary stages of peak neurogenesis (GW10–24), and altogether the organoids are most significantly similar to GW14 (**P = 0.0015, two-sided Welch’s t-test). ‘Overall’ refers to the subtype correlation observed when comparing all organoids cells to all primary cells and is shown for comparison. The last histogram shows the average gene score of each sample and error bars show s.d. Younger primary samples and organoids have a relatively lower gene score related to their marker specificity; this specificity increases substantially over time in primary cells but less so in organoid cells. Source data

a, Canonical correlation analysis from Seurat v3 was performed using reference-based integration. For this analysis, 20,000 cells were randomly subsetted from both the primary and organoid datasets and their counts matrices were merged. The primary samples were designated as the reference, and using CCA the organoid cells were projected into that reference space. A UMAP plot of the intersection is shown. The stacked histogram shows the relative contributions of each sample to each cluster. Most clusters were primarily one dataset or the other, validating the observations of limited primary subtype recapitulation in organoids. b, For the clusters with at least 20% contribution from both primary and organoid cells, differential expression was performed across all of these clusters jointly using a two-sided Wilcoxon rank-sum test. The full differential expression is presented in Supplementary Table 5, but genes upregulated in organoid cells were examined with Enrichr pathway analysis, and a summary of the top Gene Ontology terms are presented (organoid: n = 20,000 cells from 37 organoids across 4 independent experiments; primary: n = 20,000 cells from 5 individuals across 5 independent experiments). c, Canonical correlation analysis from Seurat v3 was performed using the integration-based method. For this analysis, 20,000 cells were randomly subsetted from both the primary and organoid datasets and their counts matrices were merged. A UMAP plot of the intersection is shown. The stacked histogram shows the relative contributions of each sample to each cluster. Most clusters were primarily one dataset or the other, validating the observations of limited primary subtype recapitulation in organoids. d, For the clusters with at least 20% contribution from both primary and organoid cells, differential expression was performed across all of these clusters jointly using a two-sided Wilcoxon rank-sum test. The full differential expression is presented in Supplementary Table 5, but genes upregulated in organoid cells were examined with Enrichr pathway analysis, and a summary of the top Gene Ontology terms is presented (organoid: n = 20,000 cells from 37 organoids across 4 independent experiments; primary: n = 20,000 cells from 5 individuals across 5 independent experiments). e, scAlign was performed for integration of datasets. For this analysis, 20,000 cells were randomly subsetted from both the primary and organoid datasets and their counts matrices were merged. A UMAP plot of the intersection is shown. The stacked histogram shows the relative contributions of each sample to each cluster. Many clusters were primarily one dataset or the other, validating the observations of limited primary subtype recapitulation in organoids. f, For the clusters with at least 20% contribution from both primary and organoid cells, differential expression was performed across all of these clusters jointly using a two-sided Wilcoxon rank-sum test. The full differential expression is presented in Supplementary Table 5, but genes upregulated in organoid cells were examined with Enrichr pathway analysis, and a summary of the top Gene Ontology terms is presented (organoid: n = 20,000 cells from 37 organoids across 4 independent experiments; primary: n = 20,000 cells from 5 individuals across 5 independent experiments). g, LIGER was performed for integration of datasets. For this analysis, 20,000 cells were randomly subsetted from both the primary and organoid datasets and their counts matrices were merged. A UMAP plot of the intersection is shown. The stacked histogram shows the relative contributions of each sample to each cluster. Although the clusters were well mixed, they had very diffuse marker gene expression suggesting key that biological drivers of variation were obscured by the analysis. h, For the clusters with at least 20% contribution from both primary and organoid cells, differential expression was performed across all of these clusters jointly using a two-sided Wilcoxon rank-sum test. The full differential expression is presented in Supplementary Table 5, but genes upregulated in organoid cells were examined with Enrichr pathway analysis, and a summary of the top Gene Ontology terms is presented (organoid: n = 20,000 cells from 37 organoids across 4 independent experiments; primary: n = 20,000 cells from 5 individuals across 5 independent experiments). i, MetaNeighbour was performed using unsupervised analysis to compare the clusters from primary and organoid samples. MetaNeighbour uses cell–cell similarity scores based upon neighbour voting and AUROC calculations to quantify the similarities between cells. These pairwise values were used as an input to hierarchical clustering, and almost entirely segregated primary clusters from organoid clusters. Box-and-whiskers plot shows quantification of the similarities within organoid and primary datasets versus the comparison of the two showed that the primary alone comparisons were significantly higher (organoid to organoid: ***P = 0.00078; primary to organoid: ***P = 0.00036, two-sided Welch’s t-test) (organoid: n = 20,000 cells from 37 organoids across 4 independent experiments; primary: n = 20,000 cells from 5 individuals across 5 independent experiments). The bars show range of subtype correlation with middle line indicating the mean and error bars the maximum and minimum. These results further validate our observations that there are important distinctions between the organoid and primary subtypes. j, The gene score for each of the four integration methods is presented, and all are significantly lower than for primary clustering alone (organoid subtype: ****P = 5.3 × 10⁻³⁸; CCA v3 Projected: ****P = 5.5 × 10⁻⁹⁴; CCA v3 Integrated: ****P = 2.8 × 10⁻²⁴; scAlign: ****P = 2.1 × 10⁻²³; LIGER: ****P = 2.9 × 10⁻⁹⁴, two-sided Welch’s t-test). The one method that substantially integrated the samples (LIGER) had the lowest gene score. Box-and-whisker plot shows mean score and error bars show maximum and minimum (n = 242,349 cells from 37 organoids across 4 independent experiments). The differentially expressed genes that were upregulated in primary samples from all four analyses are intersected. A substantial number of these genes were found in all four datasets, and these genes included examples that we identified from other methods in this study, including PTPRZ1, MEF2C and SATB2, validating the accuracy of our analytical methods and our main findings. Source data

a, Variance partition was run on both primary and organoid datasets across the metadata properties shown. Each dot represents a gene and the amount of variance of that gene explained by the relevant metadata property. b, ChEA analysis of type genes identified in primary cortical samples. The x-axis shows the −log₁₀(adjusted P) of the transcription factors indicated; results obtained from Enrichr datasets included a variety of experimental systems but have been shortened for ease of reading to the relevant transcription factor (n = 189,409 cells from 5 biologically independent samples; two-sided Wilcoxon rank-sum test). Type genes in organoid samples were not unified for significant transcription factor regulation. c, Violin plots of radial glia and neuron markers in primary (orange) and organoid (blue) radial glia and neurons in which width of coloured section indicates distribution of expression of each data point within a sample. In some cases, organoids show expression of multiple markers, lower expression of key markers, or similar expression to that seen in primary samples (organoids: n = 242,349 cells from 37 organoids across 4 independent experiments; primary: n = 189,409 cells from 5 biologically independent samples from 5 independent experiments). d, Dot plots from Fig. 2 shown with one colour only to avoid dot overlap. e, Lower-magnification images of PTPRZ1 and HOPX overlap as shown in Fig. 2c show domains of overlapping expression in the primary oSVZ and distinct domains of expression in the organoid ventricular zone. Validation stains were repeated independently three times. Source data

a, WGCNA networks generated from annotated primary radial glia (Methods) were applied to both primary and organoid radial glia cells. Module eigengenes shown in the heat map indicate overall higher expression in primary compared to organoid radial glia. b, Pseudoage (x-axis) versus actual age (y-axis) in PFC and V1 radial glia showing that PFC neurons are more mature than V1 radial glia. c, Box-and-whisker plot (minimum to maximum, bar at mean) across all cells within a single organoid from all organoids within this study show heterogeneity of maturation is within a single organoid and not between individuals (n = 242,349 cells from 37 organoids across 4 independent experiments). d, The parallel pseudoage analysis to the analysis in Fig. 3c is shown, but starting with organoid networks for the pseudoage calculation. Graph line shows mean pseduoage score against actual age, shading represents the geometric density standard error of the regression. The same pattern is observable, with organoids failing to recapitulate the molecular maturation of primary radial glia, though genes related to the switch from neurogenesis to gliogenesis are preserved and may account for some of the limited correlation. Source data

a, Organoid areal assignments by age, line and protocol indicate heterogeneous areal identity. b, Heat maps showing normalized module eigengene signature of each area in primary samples (with known area on the right) and in organoid samples. c, Summary of assigned area in primary samples compared to actual area. In many cases, they correspond strongly, and in others there is evidence of lack of distinction. For example, parietal cells still strongly express temporal signatures, suggesting that they have not yet been distinctly specified in primary samples, although this specification does exist in organoids. d, Box-and-whisker plot (minimum to maximum, bar at mean, error bars show s.d.) is the same comparison as shown in Fig. 4c, but across all areas (primary: n = 122,958 excitatory neurons from 5 individuals from 5 independent experiments; organoids: n = 97,531 excitatory neurons from 37 organoids from 4 biologically independent stem cell lines. In some cases there is no significant difference between strength of area signal in primary cells and organoid cells (PFC, NS (not significant), P = 0.5373), in other cases either the primary or organoid sample is significantly stronger (motor: *P = 0.0148; all other areas: ****P < 0.0001; Welch’s two-sided t-test). Source data

a, Markers of metabolic stress are expressed across cortical organoid protocols. Violin plots show both data from our experiments (1–3) and published datasets from other protocols (4–12), which have significantly increased expression of the glycolysis gene PGK1 and the ER stress genes ARCN1 and GORASP2 compared to primary samples (n = 5 individual replicates, GW14 shown). Width of the colored area indicates mean gene-expression level of each dataset and overlaid dots show each individual data point. All protocols have significantly higher expression of these three markers compared to primary samples (****P = < 0.0001, two-sided Student’s t-test). b, Single-cell sequencing identified increased expression of genes in organoids, which was validated across all stages of organoid differentiation evaluated (weeks 3–14). Validation staining experiments were repeated independently three times. Representative images from week-14 organoids differentiated using the least directed differentiation protocol. Colonies of induced PSCs also express the ER stress markers ARCN1 and GORASP2 (n = 3 biologically independent samples across 3 experiments). Scale bar, 50 μm. c, Primary cortical tissue express glycolysis and ER stress genes at undetectable levels (n = 3 biologically independent samples across 3 experiments). When tissue was cultured for one week, there was no significant increase in cellular stress (n = 3 biologically independent samples across three experiments). Scale bar, 50 μm. Source data

a, Metabolic stress network module eigengene expression across all cells is shown in box-and-whisker plots (minimum to maximum, bar at average, error bars show s.d.) across 11 datasets generated either in this manuscript or from publicly available datasets. Data are shown for expressed genes from KEGG pathway glycolysis and ER stress networks. This study: n = 242,349 cells from 37 organoids across 4 independent experiments; published datasets as annotated. b, The same box-and-whisker plots are shown for organoids (n = week 3: 38,417 cells, week 5: 26,787 cells, week 8: 11,023 cells, week 10: 50,550 cells, week 15: 2,722 cells, week 24: 4,506 cells from 4 independent experiments) and all primary ages (n = GW6: 5,970 cells, GW10: 7,194 cells, GW14: 14,435 cells, GW18: 78,157 cells, GW22: 83,653 cells from 5 independent experiments). ER stress and glycolysis networks decrease over time in primary samples but decrease less in organoids and are significantly higher in most organoid stages than in primary samples. Significance was calculated for each organoid sample with respect to each primary sample, and a one-sided Welch’s t-test was performed (to evaluate whether organoid expression was higher than primary). All comparisons were either not significant (ns) or significant with ****P < 0.0001. c, Cellular stress genes are expressed at low levels during human cortical development. GW13 and 17 samples were stained for the glycolysis gene, PGK1, and showed little expression at either age. The ER stress gene ARCN1 had little expression at either age, but there was modest expression of the ER stress gene GORASP2 at GW13 that decreased by later neurogenesis. Staining validation studies were performed independently four times. d, Dissociated primary cells were cultured for one week. Across five independent studies, there was no detectable expression of the glycolysis gene PGK1, but the ER stress genes ARCN1 and GORASP2 showed significantly increased expression. e, Immunostaining of primary aggregates (n = 5 biologically independent samples), which express markers of oRG cells (HOPX and SOX2), IPCs (TBR2) and neurons (CTIP2). Aggregates also had increased cellular stress indicated by PGK1, ARCN1 and GORASP2 staining. Violin plots show expression level and data distribution for each marker in primary cells, primary cells after organoid transplantation and primary cells after being aggregated together. The expression of PGK1 and GORASP2 are increased in post-transplanted primary cells from the organoid as well as in primary cell aggregates. Cell types and physical distribution in the primary aggregate are shown. Scale bar, 50 mm. Representative image shown (n = 3 replicates). Source data

a, FACS plots showing dummy infection (left) and transplanted organoids (right) in terms of the their GFP signal (x-axis) versus sidesscatter (y-axis). Cells in the gated region were collected (% of parent written on plot) and sequenced for transplantation 2.5 weeks after incubating in the organoid, representative plot shown on right, n = 5. b, Immunohistochemical validation that cells infected with GFP virus were all SOX2-labelled progenitor in cells dissociated from primary cortical tissue GW14–20. Scale bar, 50 mm, representative image shown (n = 5 replicates). c, An additional example of primary cell integration into organoids after transplant, in which the primary cells integrate into organoid rosettes (n = 7 primary samples into 21 organoids across 2 independent studies). d, t-SNE of pre- and post-transplant primary cells, as well as the cluster designations. Many cell types represented in pre-transplanted cells are not present in the post-transplant population. e, Subtype similarity correlation between pre-transplant, post-transplant, and primary aggregate samples. Includes plot (bar is average subtype correlation, error bars are s.e.) as a replicate of the experiment in Fig. 5b, validating that at older organoid ages (week 12) the post-transplanted cells are still significantly impaired in their subtype specification (****P = 1.46 × 10⁻¹¹, n = 2 primary biologically independent samples into 2 organoids in addition to n = 5 biologically independent samples into 10 organoids in Fig. 5, two-sided Welch’s t-test). Primary aggregates are significantly impaired in their subtype specification (**P = 0.0016), but are significantly better than post-transplanted primary cells (**P = 0.0037). This may be related to non-neural populations in the aggregates. f, Transplanted organoid cells were visualized in the mouse cortex after 2 and 5 weeks post-transplant (n = 13 independent mice transplanted with 14 organoids derived from 2 induced PSC lines across 2 independent experiments). Human cells were visualized by GFP and human nuclear antigen (HNA) expression. Organoid-derived cells expressed markers of progenitors (SOX2 and PAX6), neurons (CTIP2, SATB2 and NEUN) and astrocytes (GFAP and HOPX). Mouse-derived vascular cells (laminin and CD31) innervate the organoid transplant. g, After 2 weeks post-transplantation, organoid cells showed reduced expression of the glycolysis gene PGK1 and ER stress genes ARCN1 and GORASP2 (n = 6 transplanted mice stained with each marker independently from 2 induced PSC lines across 2 independent experiments). h, Subtype correlation analysis of pre- and post- transplanted organoid cells shows an increase in oRG subtype identity (similarity to primary cluster 26) and in newborn neurons (similarity to primary cluster 22). Source data