Cellular stress in brain organoids is limited to a distinct and bioinformatically removable subpopulation

Organoids enable disease modeling in complex and structured human tissue, in vitro. Like most 3D models, they lack sufficient oxygen supply, leading to cellular stress. These negative effects are particularly prominent in complex models, like brain organoids, where they can prevent proper lineage commitment. Here, we analyze brain organoid and fetal single cell RNA sequencing (scRNAseq) data from published and new datasets totaling over 190,000 cells. We describe a unique stress signature found in all organoid samples, but not in fetal samples. We demonstrate that cell stress is limited to a defined organoid cell population, and present Gruffi, an algorithm that uses granular functional filtering to identify and remove stressed cells from any organoid scRNAseq dataset in an unbiased manner. Our data show that adverse effects of cell stress can be corrected by bioinformatic analysis, improving developmental trajectories and resemblance to fetal data.


Introduction
Organoids are 3D stem cell cultures that enable human tissue modeling with unprecedented structure and complexity (Eiraku et al., 2008;Kadoshima et al., 2013;Lancaster et al., 2013;Paşca et al., 2015;Qian et al., 2016). At the same time, single cell transcriptomics has become widely used for their characterization. Alongside these recent technological breakthroughs, it has become clear that 3D tissue culture is affected by limited oxygen and nutrient supply to the center of the tissue.
As most models lack functional vascularization (Garreta et al., 2021), and therefore rely on limited passive transport across the tissue, diffusion-limited hypoxia is an intrinsic problem in organoids.
Brain organoids are among the most complex and physically largest organoids and are therefore most affected by the limited nutrient supply of the center . Nevertheless, this problem has only been recently addressed in brain organoids Giandomenico et al., 2019;Mansour et al., 2018;Qian et al., 2020) and its extent is still debated.
It therefore remains an open question if stress is a global or a local issue, thus, how far it affects the 3D tissue culture model. A recent paper claimed that in vitro conditions lead to a pervasive stress across the whole organoid, causing immaturity, misspecification, and dissimilarity to fetal tissue . These observations are in contrast to the previous understanding of spatially limited stress . This raises the question: how should we handle the data affected by an artificial stress signature? It is unclear if stress is an 'acute' signature on top of a cell's original state, or if stress is leading to a completely different cell fate. While the same stress pathways are also active in the fetal brain, reports disagree on whether it is equivalent to those observed in vitro Gordon et al., 2021).
These approaches aim to increase nutrient supply, but neither is currently as scalable as the standard organoid culture which therefore remains the mainstay of organoid research. Until a widely applicable and scalable experimental solution emerges, tissue health and cellular stress persist as a problem for the field.
While most large single-cell RNA-seq studies on diverse brain organoid systems reported glycolytic or ER stressed clusters Tanaka et al., 2020;, there is no consensus on how to identify them, what is happening in these cells, and what are the consequences of stressed cells on the organoid? To measure the prevalence and consequences of stress in brain organoids, we analyzed differentiation, maturation, and identity of ~160,000 single cells from newly presented and published cortical and cerebral organoid (together: brain organoids) datasets.
We found stressed cells in all organoid samples, forming a distinct subpopulation. Stressed cells showed widespread transcriptional changes beyond stress pathway activity, we therefore call it 'stressed-state'. We do not find this stress-state in vivo, therefore it is likely an artifact. Eliminating artificial cell populations is essential to truly recapitulate in vivo conditions. As stressed cells are currently unavoidable, we developed granular functional filtering (Gruffi), an unbiased computational algorithm to isolate stressed cells. Gruffi can clarify developmental trajectories and increase similarity of in vitro fetal datasets.

Results
A distinct population of ER stressed-and glycolytic-cells exist in all analyzed organoids We reanalyzed recent, landmark single cell transcriptomics studies and performed new experiments Khan et al., 2020; 2019) (Fig 1A) to answer 3 questions: What stress pathways are active in organoids? Does stress occur in all or only certain organoid protocols? Is cellular stress limited to a group of cells or is it pervasive?
We included all those mature, wild type samples that were prepared on the 10X Chromium singlecell platform starting from .fastq files using the same pipeline (methods). To focus on stress in the neural lineage, we removed all cells that are not part of brain development and are a result of mispatterning, sometimes observed in organoids. For a proportional representation of datasets, we sub-sampled ~160K from a total of 300K cells.
Cellular stress can lead to a perturbation of essential processes, thus affecting cell quality in scRNAseq. Therefore, we applied a minimal filtering, keeping all cells with >500 genes, less than 20% mitochondrial-and 30% ribosomal-reads (Ilicic et al., 2016;Luecken and Theis, 2019). This resulted in a median depth of 3651 UMI/cell (methods). We integrated and analyzed the resulting datasets in Seurat (v4), and found the previously reported cell types (Fig 1B). The UMAP separated dividing cells and glia cells from neurons (horizontally) and excitatory-from interneurons (vertically).
Besides, there were multiple clusters in the center of the UMAP, which were less well defined by marker gene expression.
To better understand the nature of stress in these cells, we analyzed all significantly enriched genes in the stress-clusters. We found that, in either cluster, more than half of the top 50 coding genes were part of 'response to stress' (GO:0006950), and apoptosis related terms were among the strongest enriched terms (Table S1). The biggest enriched stress pathways were 'regulation of cell death' (GO:0010941), 'response to hypoxia' (GO:0001666), and 'response to endoplasmic reticulum stress' (GO:0034976). Surprisingly, metabolic terms were both among the strongest and largest enriched terms, highlighting that metabolic shift is a hallmark of stressed cells in organoids.
We then calculated the GO-term enrichment within the 150 strongest enriched coding genes of both stress clusters together, and visualized these on the protein-protein interaction (PPI) map (Fig 1E,methods). We highlighted enriched GO-terms (FDR < 5e-7) forming connected PPI clusters, revealing the interplay of glycolysis, hypoxia, unfolded protein response, and translation with the general stress response (Fig 1E, Table S1). To identify all genes co-regulated with stress, we applied scWGCNA (single cell weighted gene co-expression network analysis, (Morabito et al., 2021)) and found 12 gene modules (Appendix Figure S1A), one of which was specific to stressed cells (Stress module Fig 1F). Gene set enrichment analysis (GSEA) on the stress module identified the strongest enrichment for "response to hypoxia", "cellular response to ER stress", and GO-terms of glycolytic processes (Fig 1G, Appendix Figure S1B). To test whether stress occurred in all samples, we quantified the contribution of each dataset to the stress clusters. We found that all datasets contained stressed cells, to a generally high (median 13%), but highly variable fraction (50% CV, Fig 1F). Thus, stressed cells are a general feature of organoids regardless of conditions, lab of origin or protocol used.
While the initial clustering-based approach identified cells with stress signatures, it had three major limitations. First, while some clusters are too large and comprise mixed populations (Appendix Figure   S1C), others can be too small to find marker genes by differential gene expression analysis (DGEA).
Second, cluster boundaries are often not well-defined, especially when dealing with developmental trajectories. Third, the resulting limitations in DGEA obstruct the identification of stress genes so that results vary by the dataset and parameters used. Together, these limitations affecting DE could explain why previous studies identified disparate gene sets, like 'Glycolytic cells' in Nowakowski et al., 2017) vs. 'ER stressed cells' in Tanaka et al., 2020).
To overcome these issues, we tested different clustering resolutions (Appendix Figure S1D). As none of these could separate the distinct populations of cells within 'Stressed Neurons' (Fig 1B), we concluded that a new approach is needed to identify and exclude stressed cells.

Granular functional filtering identifies stressed cells unbiasedly
Functional scoring highlights cellular stress regardless of cluster boundaries To universally identify stressed cells, we established a sample-and data-independent definition for stressed genes. Using gene lists from well-characterized pathways defined as GO-terms (methods), we aggregated information from all genes per pathway by an expression-scoring method widely used for cell cycle scoring . Therein, we downloaded gene lists per GO-term from Ensembl, calculated their average expression and normalized it to randomly sampled control genes of matching expression level (methods). Finally, we evaluated if functional scoring helps to characterize stressed cells. We found that high scores for 'glycolytic process' (GO:0006096) and 'response to endoplasmic reticulum stress' (GO:0034976), were the strongest signatures of stressclusters (Fig 2A) and provided clearer separation between stressed and non-stressed cells than cluster boundaries. Importantly, high scores marked mostly overlapping cell populations. The coactivation of additional scores, such as 'response to starvation' (GO:0042594, Fig EV1F) and 'cellular response to hypoxia' (GO:0071456, Fig EV1G) corroborated a complex stress-identity.
Comparing neuronal and glial cell types revealed that all non-dividing glial cells showed higher ER stress scores (Fig 2B). We therefore designed our algorithm to accommodate for cell-type specific background when identifying stressed cells.

Change of cellular identity in stressed cells
Besides increased stress-gene expression, stress clusters were characterized by low expression of pan-neural markers (NEUROD6, DCX, MAP2, NCAM1, ELAVL4; Appendix Figure S1E). To test if this marker depletion is also reflected by a general change in glial or neural fates, we extended our scoring approach. We calculated scores for the two cardinal cell states in neural development, neurogenesis (GO:0022008), and gliogenesis (GO:0042063). Both terms were depleted in stressed cells (Fig 2A). Compared to both glial and neural clusters, stressed cells also showed remarkably low scores for 'cell differentiation' (GO:0030154), and 'forebrain development' (GO:0030900), suggesting that stressed cells are in a more 'basal' state compared to differentiating neurons (Appendix Figure S2 AB).Thus, chronic stress in organoids comes at the expense of neurogenic cell differentiation and leads to a switch in cell fate. We therefore refer to the reduced expression of marker genes as the 'stress identity'. As stressed cells lose both glial and neuronal identity and stop differentiating, it explains why these cells occupy the middle of the UMAP (Fig 2A).

Granular evaluation overcomes noise inherent to single-cell data
Single-cell gene expression measurements are inherently noisy. While GO-scores are computed across multiple genes per cell, these may still suffer from high variability and noise. Indeed, some cells within the stress-clusters showed low stress-scores, even if clustering together (Fig 2A). At the same time, sporadic cells in well-defined cell types showed high stress-scores. These cells expressed stress genes inconsistently, but expressed respective cell type markers, which are otherwise absent in stress-identity.
To overcome variability in single cell measurements, one can either denoise the data, e.g., by modelbased imputations, or group cells and evaluate them together. Many different imputation methods have been developed recently, however, imputed values often vary (Hou et al., 2020), and they can induce false signals (Andrews and Hemberg, 2018). This is probably due to the complexity of the imputation problem. We therefore took a grouping approach where we partitioned cells into groups of 100-200 cells by ultra-high-resolution SNN-clustering in PCA-space (Methods) resulting in small groups of cells, that we term granules.
The ultra-high-resolution clustering approach can overcome the problems of boundaries by breaking down the data into minute groups of cells. To get sufficient coverage for robust gene scoring, and because clustering creates some very small granules, we added a reclassification step, where cells in granules with <30 cells are reassigned to the closest granule above threshold (Methods).
To test the granular approach, we compared stressed cells identified by Gruffi's granular method (gSC) and stressed cells identified on single-cell scores (scSC). We compared cells only identified by either, both or neither of the approaches. First of all, scSC were evenly scattered across all clusters, while most gSC were close to stress-clusters and the cells identified by both methods (Appendix Figure S2C). By definition, single-cell selection on stress-scores identified the cells with highest stress-scores. However, scSC showed less of all other features defining stress-identity: lack of cell differentiation, lower mitochondrial and higher ribosomal mRNA content. Granular identification found cells that shared these features more with stress-identity (Appendix Figure S2D-I). We implemented both methods in Gruffi, but we concluded that the granular approach is more suitable if one aims to exclude stress-identity, whereas the single cell approach is more suitable if one aims to simply find cells with the highest stress gene expression, but otherwise properly specified cells.

Granular functional filtering (Gruffi) isolates and removes stressed cells
As clustering-based identification failed to detect stressed cells specifically and robustly, we built on the concepts above and developed Granular Functional Filtering or Gruffi (Fig 2B). Gruffi takes a number of Gene Ontology Pathways (1) to obtain corresponding gene sets (2), and computes cellwise GO scores (2). At the same time, it identifies a suitable resolution (I), clusters cells in granules (II), and reassigns cells of too small granules (III). Merging these, it then computes the multiple granule scores (4), estimates a threshold separating stressed and non-stressed cells (5) and assigns a 'stress' label integrating multiple scores (6).
To uniformly determine the prevalence of stressed cells across organoids and protocols, we applied Gruffi to the integrated organoid dataset. After pathway scores calculation, we estimated the optimal clustering resolution, which was 1009 granules with a median of 154 cells (Fig 2C). Stress identification must be robust across all datasets, therefore we incorporated 3 scores: the two most specific pathways: 'glycolytic process' (GO:0006096) and 'response to endoplasmic reticulum stress' (GO:0034976), and a negative filter on 'gliogenesis' (GO:0042043) score accounting for the higher native expression of ER genes in glia. Addition of further negative filters did not improve identification; however, we implemented this option in our algorithm.
Gruffi then calculated the average and cell number normalized functional scores per granule (Fig   2D) resulting in a 3-dimensional functional score for each granule. This combinatorial approach is flexible to the type and number of scores used, which may be useful for applications beyond its original scope. Next, Gruffi estimated the thresholds separating stressed from normal cells, accounting for score variability among non-stressed cells (Methods). At this point the retrieved thresholds can and are advised to be inspected and refined via the implemented interactive Shiny App interface. Throughout our analysis presented here, we did not further adjust the by Gruffi initially estimated thresholds. Finally, combining all thresholds, Gruffi classified stressed and non-stressed cells (Fig 2F), which largely overlapped with the expression of key stress markers, and with high stress scores (Fig 1C). While this suggested a correct identification, we continued the in-depth analysis of stressed cells.
Stress has a profound, yet limited transcriptional effect and stress-identity is not present in vivo.

Stressed cells show a profound transcriptomic change
As all samples had stress-identity cells, we searched for their defining features and their consequences on the whole organoid. Stressed cells fell in two distinct clusters with either a more glial or a more neural character (Fig 1). We therefore divided Gruffi's stressed cells into these two categories to investigate the heterogeneity of stress response in organoids (methods). We quantified the total expression of mitochondrially encoded genes and of ribosomal mRNAs, which correspond respectively to ~2% and ~10% of the total transcriptome, respectively. These are widely used to assess quality and cell state in scRNAseq (Luecken and Theis, 2019). We hypothesized that chronic hypoxia and glycolysis diminish the need for oxidative phosphorylation, which may translate to fewer mitochondrial UMIs. Indeed, stressed neurons showed 52%, whereas stressed progenitors showed a 25% reduction in mitochondrial read fractions, as compared to their non-stressed counterparts ( Fig   3A). In contrast, ribosomal mRNA fractions were 40% and 23% higher, respectively (Table S2), perhaps to compensate for ER-dysfunction (MWW p.val < 2e-16 in all cases)

Increased catabolism in stressed cells
Because we saw the decrease in mitochondrial mRNAs, we looked for transcriptomic signatures for the active degradation of mitochondria. To that end, we applied Gruffi's scoring method for relevant GO terms, and found opposite signs: both stressed groups showed positive scores for 'autophagy of the mitochondrion' (GO:0000422) on average, while normal cells do not (Fig 3A). At the same time the groups scored similarly for 'autophagy' (GO:0006914) (Appendix Figure S3A, Table S2).

Translation in ER stressed cells might lead to protein degradation via the ubiquitin-dependent ERAD
pathway. Therefore, we calculated the corresponding score (GO:0030433), and as for mitophagy, found that stressed cells scored positively, while normal cells scored negatively (Fig 3A, Table S2).
Altogether, these changes show that stress induces major changes to the cell's physiology and metabolism that go beyond acute stress response.
Finally, we analyzed cell types by clustering PROGENy's activity-score across signaling pathways and clusters (Schubert et al., 2018). We found that stressed cells form an outgroup and are marked by the upregulation of Hypoxia and VEGF pathways, and the downregulation of the PI3K pathway, highlighting oxygen deficiency and quiescence (Fig 3B, methods). To ask how stressed cells identified by Gruffi differ from those identified by the naïve approach of identifying stress by clustering (Fig 1), we subcategorized stressed cells identified in only one of the two classifications. The complete lack of separation of cells found by the naïve approach contrasted the salient stress features found by Gruffi (Fig 3B).

The presence of stressed cells does not affect specification and maturation of non-stressed neurons 1
A previous study reported that stress in organoids leads to impaired cell-type fidelity, and incomplete maturation as a global phenomenon in organoids . To test those effects, we compared cell types in the organoid datasets to fetal cortical data of comparable age (de la Torre -Ubieta et al., 2018; Polioudakis et al., 2019) (methods). We defined the fetal brain as the reference data and constructed modules from co-expressed genes ( Table S3, methods) as described before (Trapnell et al., 2014). The resulting 65 aggregated gene modules were then used to calculate Pearson correlation across clusters and then visualized in a heatmap (Fig 3C). Major cell types (excitatory neurons, interneurons, and progenitors) formed the largest clusters across the dataset.
Most organoid cell types pairwise best matched the corresponding fetal cell type, indicating that organoids undergo proper cell type specification unlike suggested previously .
Stressed neurons, however, were uncorrelated to all cell types, except stressed progenitors (Fig 3C). This dissimilarity to all fetal and organoid cell types, along with the analyses above, indicated that stress neurons lost most of their identity and formed a new transcriptional cell-state.
Interestingly, while stressed progenitors showed a similarity to stressed neurons, they also showed a strong progenitor identity, suggesting either an increased robustness or more distinct transcriptome of the progenitor state. Altogether we found no evidence of imparied cell-type fidelity in organoids, as cell types in organoids match respective cell types in vivo, and that stressed cells show little resemblance to cell types found in the fetal cortex.

Stressed cells in organoids have no fetal counterpart
Stressed cells might also exist in vivo, even if they did not form a recognized cluster in published studies. While some previous reports have argued that stressed cells similar to organoids exist in vivo (Gordon et al., 2021;Tanaka et al., 2020), others claim that it is an artifactual population specific to organoids . Therefore, we integrated the fetal brain dataset with a matching number of randomly downsampled cells from the organoid dataset. Using the published fetal cell type annotation, we found that the organoid dataset was generally well recapitulating the fetal data (Fig 3D). At the same time, CGE and MGE (caudal and medial ganglionic eminences) interneuron differences and deep layer excitatory neuron differentiation were clearer in fetal data.
Interestingly, there were two populations entirely of organoid origin: a population near the neural trajectory (I) and a progenitor population (II, Fig 3D). To test for stress identity, we calculated stress scores as before, and found that precisely these populations score high for ER stress and glycolysis (Fig 3E). To quantify stressed cells and validate our method on both datasets, we ran Gruffi on the integrated object. This identified stressed cells almost exclusively in the two aforementioned populations (Fig 3G) and almost exclusively in organoid samples (Fig 3H). Altogether, our analysis suggests that stress-identity is an in vitro artifact and there is minimal to no stress-identity in vivo.

Stressed cells do not affect the maturation of other cells, but their removal
improves data quality and interpretability.
The removal of stressed cells reveals clear trajectories that recapitulate fetal neurodevelopment Removing stressed cells might create a better model of the fetal brain development. Stress genes majorly contribute to variable genes, which determine both clustering and visualization (Luecken and Theis, 2019). We therefore removed all stressed and low-quality cells, reidentified variable genes, and recalculated all dependent representations (PCA, UMAP, snn-graph) and downstream analysis with identical parameters.
Starting from the progenitors, the resulting UMAP revealed 3 trajectories ("E", "I" & "MB" in  (Fig 1B). Instead, midbrain cells were linked to their progenitors only by a small, separated population of the 'yellow' cluster. After stressed cell removal, these trajectories now lead to distinct populations of mature neurons, as opposed to the continuum of connected clusters before Gruffi (before stress removal). Notably, this lineage separation recapitulates fetal neurodevelopment.
To ensure the robustness of our approach we repeated the analysis on a single dataset (Fig EV2A to E), as well as downsampling the full dataset to ~8200 cells, which is the typical output of a single 10X experiment (Fig 4C-D). Both analyses revealed that expected lineage trajectories are missing or broken in the UMAPs before Gruffi (midbrain in Fig 4A, and excitatory neurons in Fig 4C), but they are correctly recovered after running Gruffi (B, D). Correct and continuous trajectories in low dimensional representations are essential for most pseudotime methods that use these as basis for pseudotemporal cell assignment.
We then tested if our method can be applied to independent datasets by re-analyzing the data from a recent publication reporting a large "undefined" cluster (Samarasinghe et al., 2021). We ran Gruffi on the precomputed R data object obtained from the authors and marked stressed cells (Fig EV3A   and B). This dataset is particularly well suited to demonstrate the versatility of Gruffi, as it contained secretory choroid plexus cells, which are undergoing high physiological ER-stress due to secretion.
We therefore derived a choroid plexus score as an additional negative filter score. As no GO-term exists for choroid plexus, we turned to the recent publication of choroid plexus organoids  identified marker genes, and derived the choroid plexus score (methods). Gruffi labeled 74% of the "Undefined" cluster as stressed and an average of 1% of cells from other clusters ( Table   S4). We then performed DGEA on stressed vs. non-stressed cell and GO-term enrichment on all enriched coding genes (methods). Visualization of enrichments on the protein interaction network using STRING (Fig EV3C) showed that apoptosis, stress, unfolded protein response, and hypoxia dominate cells identified as stressed in this dataset as well.
Organoids show proper cell type specification along the excitatory lineage A previous study reported that stress in organoids not only leads to impaired cell-type fidelity but also incomplete maturation . To investigate whether stressed-identity cells affect specification along the excitatory lineage, we compared progenitor and neuron signatures (methods , Table S5). This led to the expected bimodal separation of the fetal samples (Fig EV4A).
As this dataset was used to derive the signatures, we confirmed the separation of neurons vs progenitors using an independent fetal dataset with multiple samples around mid-gestation (Fig 4C,   Fig EV4B). Next, we asked how organoid samples separate using those signatures. We calculated signature scores on Bhaduri et al. organoid datasets and reproduced the previously reported lack of specification (Fig 4D). In contrast, when testing the other datasets analyzed in this study, we detected a fetal-like specification (Fig 4E, Fig EV4C-F), suggesting proper specification in most organoid datasets.
Nevertheless, organoid cells still separated less than fetal cells, and more cells were scoring low on both progenitor and neural axes. As stressed cells are characterized by the lack of both glial and neural signatures (Fig 2A), we hypothesized that stressed cells may populate the area between neurons and progenitors. After annotating stressed cells, we found two populations in between progenitors and neurons: stressed cells and intermediate progenitors, or IPCs. After stress removal, most of those remaining cells are IPCs (Fig 4F, in red), which are indeed a transitory stage between glia and neurons. Altogether, we find no evidence for general misspecification in organoids. Instead, progenitors and excitatory neurons properly separate, while two specific populations, IPCs, and stressed cells, lack specific neuronal or progenitor signatures.

The presence of stressed cells does not affect the maturation of other cell types
Besides a lack of cell type specification, incomplete maturation due to stress was previously suggested . As a positive control for maturation, we took two media conditions that affected maturation . We grew pairwise matched samples in two different media conditions, then analyzed together (Fig 4G). We calculated the pseudo-temporal trajectory of the excitatory lineage and graded each dataset for maturation along this trajectory (Fig   EV4G, methods).
Our results showed that the fraction of stressed cells does not explain maturation differences (Fig   4H), but the choice of media does: low-nutrient media improved neural maturation. To understand the impact on maturation on single cell level, we plotted individual cells along the maturation trajectory, split by media condition (Fig 4I). Additionally, to assess expression changes associated with mature states we generate scores for mature neurons (Fig 4I, Fig EV4H to J, Table S6). This revealed that while organoids grown in either condition had abundant excitatory neurons, those in HN media remained mostly immature, while organoids grown in LN media contained more mature neurons (Fig 4I, Fig EV4J to L). In sum, the presence and abundance of stressed cells in a sample have negligible effects on neural maturation, while measurable differences arise by the choice of media.

Discussion
Brain organoids generate complex, structured tissue in vitro (Eiraku et al., 2008;Kadoshima et al., 2013;Lancaster et al., 2013;Paşca et al., 2015;Qian et al., 2016). Besides their tremendous potential for modeling human diseases (Sidhaye and Knoblich, 2021), it is critical to understand and account for their limitations. Here, we showed that a population of stressed cells exists in all analyzed organoid samples, and that this is a biologically distinct population, which is not found in vivo. We provide an in-depth analysis of these cells that hopefully will help to decipher the needs of 3D tissue in culture.
While an experimental solution is the end goal, stress is present in published and current experiments. To tackle this issue, we developed Gruffi, an in silico tool to bioinformatically identify and remove these cells, based on stress pathway activity scoring. Gruffi comes with a graphical and interactive interface. It integrates into a typical single-cell analysis workflow using Seurat, but can be used in other pipelines as well. The resulting stress-decontaminated samples displayed a clearer representation of the fetal neural development and showed higher similarity to in vivo samples. Even if future organoid protocols may resolve cellular stress, earlier published data still suffers from stress, which negatively impacts data integration. Gruffi, however, can recover such data for comparison, reducing the need for performing new experiments.
We observed diverse stress pathway activity, and it is important to understand how they are connected on a cell biological level. Our results are compatible with earlier observations that the organoid core, but not surface, is hypoxic , explaining why stress characterizes only in a defined set of cells. The central role of hypoxia can explain the other transcriptional shifts.
The lack of oxygen triggers a metabolic shift, from oxidative phosphorylation to anaerobic glycolysis.
Hypoxia also triggers ER stress, in two ways. First, glycolysis is much less efficient in energy production, leading to energy depletion, and consequently stronger metabolite transport is needed.
These transporters are secreted via the ER-pathway (Loike et al., 1992), triggering unfolded protein response (UPR) (Lee et al., 2020). At the same time, the depletion of energy leads to a pH imbalance, affecting organelles that rely on ATP-dependent transporters for ion homeostasis (Chiche et al., 2010).
Our results are consistent with a previous observation that acute hypoxia in cortical spheroids triggers a strong ER-stress response (Pașca et al., 2019). However, a simplistic, one-dimensional distance-to-surface model of nutrient availability cannot explain the heterogeneity of stress marker expression. It is an interesting future direction to determine different cellular niches, e.g. by local variation in oxygen and nutrient levels. Similarly, an interesting question for future studies is, how cellular heterogeneity leads to the differential expression of stress markers in close neighboring cells.
Importantly, stress identification is just the first application and proof of principle for granular functional filtering. This flexible framework can be extended to many other applications in single-cell analysis. As long as a group of cells form an identity, so that they group together in any low dimensional representation, and coexpress a defined geneset (GO-terms, KEGG-pathways, etc.), the cells can be identified, studied and removed. Here, we applied Gruffi to remove stressed cells from brain organoid datasets, but we think that there are many other applications possible, such as selecting cells from a lineage or cells responding to a treatment. Currently, brain organoids are the largest and longest-cultured 3D organoid systems and are therefore particularly affected by stress.
As 3D tissue models and investigations become ever more sophisticated, cell culture induced artifacts are more important to account for. Therefore, we expect that our approach will find many applications beyond its original scope.

Data and code availability
The single-cell RNA-sequencing data have been uploaded to Gene Expression Omnibus (GEO) under reference number GSE1XXXXXXXXX (ncbi.nlm.nih.gov/geo/query/acc.cgi?acc= GSE1XXXXXXXXX). We used publicly available raw sequencing data from the following publications  and obtained wild type patient data from the authors of , complying with ethical and data safety requirements (Khan et al., 2020;Samarasinghe et al., 2021). The Gruffi package will be available under github.com/jn-goe/gruffi. The code for analysis will be accessible on     midbrain neurons, and cortical excitatory neurons. Clustering and colors are the same as in (Fig 1A).
(C-D) We repeated the analysis on a smaller subset of the data downsampled to ~8200 cells, a typical outcome of a single 10X experiment (this subset does not contain midbrain neurons). Note that almost any lineage relationship could be inferred from the representation before Gruffi (  used for this analysis (as in Fig 1E). Enriched genes are marked red. DGEA and all enrichment terms are in (Table S6).

Fig EV4 -Proper specification and maturation in non-stressed cells (A)
Progenitor (x-axis) and excitatory neuron scores (y-axis) calculated in fetal brain datasets across brain development (materials and methods, Bhaduri et al 2020). Early datasets from Carnegie stage 13 (CS13) are enriched in progenitors, while late datasets from gestational week (GW) 18 and 22 show both neurons and progenitors. (B) Progenitor (x-axis) and excitatory neuron scores (y-axis) applied to multiple datasets of mid-gestation (Fig 4E) show reliable separation of neurons and progenitors, while only intermediate progenitors (IPCs, red) cluster in between the two populations.  (Table S7) is specifically enriched in the UL-EN cluster (Cluster 8 in Fig 4I). cells identified by Gruffi's granular method (gSC) and stressed cells identified on single-cell scores (scSC). Cells are colored by clusters as in (Fig. 1), and separated by stress identification classes  Table Descriptions   Table S1 A small, clearly distinct population of ribo-high cells was found to come from a single dataset, and constituted mis-patterned, non-neural tissue. This table quantifies the distribution of these ribo-high cells (>30% ribosomal reads) across datasets and clusters.

Table S2
Comparison of ribosomal, mitochondrial reads and stress scores in stressed and unstressed clusters. Table S3 We calculated modules of co-regulated genes on the fetal reference dataset , that were subsequently aggregated to correlate fetal and organoid data.

Table S4
Samarasinghe DGEA results and stressed-cell enriched coding genes. STRING DB's GO -Biological Process enrichments on the above gene lists. Terms highlighted in Fig EV3C and permalink.

Table S5
We calculated differentially expressed genes enriched in excitatory neurons (first tab) and progenitors (second tab) in the fetal reference dataset . The top 30 genes were selected for progenitor and excitatory neuron scores, respectively (third tab). carrying the Yamanaka reprogramming factors OCT3/4, SOX2, c-MYC and KLF4 factors. All cell lines were used within 10 passages from last STR profiling and tested regularly for mycoplasma contamination. We additionally used the above cell lines (177 and 178) for the media comparison experiments (Fig. 4). The cell lines were evaluated and cultured as the HipSci lines. Briefly, cells were seeded on vitronectin (Stemcell Technologies, cat#07180) coated plates and fed every day with E8 essential media. Cells were passaged as single cells using Accutase (Sigma) with Revitacell cell supplement (1/100, Invitrogen, cat#A2644501), and grew until 70% confluency, then we replated.
Cultured cell lines were routinely tested for mycoplasma contamination by PCR (Janetzko et al. 2014).
Organoid culture Organoids were generated as described in (Esk et al. 2020). Briefly, 150 uL/well of Essential 8 media supplemented with RevitaCell (1/100) containing the corresponding cell suspension for 9000 cells/well were plated for each cell line using low attachment 96-well plates (Sigma CLS7007) . Briefly, the protocol entailed the following steps: On day 3, media was replaced to Essential 8 media and from day 6 on, embryoid bodies were transferred to neural induction media (NI) and 200uL/well was exchanged every day. On day 10, when embryoid bodies are about 500-600um in thickness and neuroepithelium is evident, the aggregates were transferred to 10 cm dishes and embedded in matrigel (MG) droplets. On day 13 and 14, the media was changed to Improved Differentiation Media without ascorbic acid (Imp-A) containing 3 uM CHIR. After that, the media was replaced every 3-4 days. On day 19, the dishes were transferred to an orbital shaker. On day 25, the organoids were fed with Improved Differentiation media with ascorbic acid (Imp+A) and the media was replaced every 3-4 days. On day 40, the two different culture methods (Brainphys & Imp+A) diverged. The "Imp+A" organoids were further cultured in Imp+A supplemented with 1%MG, BDNF (20ng/mL), GDNF (20ng/mL) db-cAMP (1mM). The "Brainphys" (BP) organoids were gradually transitioned to BP media in 3 feeding steps: 75%-25%: 50%-50%, 25%-75% (Imp+A & BP). From that point they were cultured in BP supplemented with 1%MG, BDNF (20ng/mL), GDNF (20ng/mL) db-cAMP (1mM). Single-cell analysis

Media composition
We first aligned reads to the GRCh38 human reference genome with Cell Ranger 3.1 (10x Genomics) using pre-mRNA gene models and default parameters to produce the cell-by-gene, Unique Molecular Identifier (UMI) count matrix. UMI counts were then analyzed in R, using the Seurat v4. We filtered for high quality cells based on the number of genes detected (>500). Thereafter, expression matrices of high quality cells were normalized ("LogNormalize") and scaled to a total expression of 10K UMI for each cell. Regression of variables at this step did not improve clustering results, hence no variables were regressed nor removed.

Non-neuronal cell exclusion
Before integration datasets were checked for quality, as certain IPS lines are prone to misdifferentiation. As a consequence, multiple datasets included non-neuronal cells that would interfere with the CCA integration, henceforth we applied initial filtering for CNS cells. To exclude non-neuronal cells we used a recently published fetal organ atlas (Cao et al. 2019). Processed data was downloaded and cell type annotation was modified to reflect major cell types for a basic classification. All CNS cell types were grouped together under one annotation to determine properly specified clusters. Next, an xgboost classifier was trained to distinguish major cell types on the RNA assay data using the top variable genes of the fetal dataset with parameters determined by crossvalidation. This classifier was applied to each dataset: 1. Datasets were pre-processed individually and clustered in UMAP space; 2. The expression of the RNA assay was used to classify each cell according to the cell groups of the training dataset; 3. Classification was summarized per cluster and all clusters that were not classified as CNS cells were filtered out. The cleaned datasets were used for CCA integration.
To establish the maximal mitochondrial-, and ribosomal RNA fractions, we plotted these against feature counts and each other, and set thresholds to remove extreme outliers. A group of cells showed a distinctly high ribosomal fraction (>30%). We found that these cells correspond to one cluster coming from one dataset (Velasco organoid 21, Table S1) and are non neural in gene expression. We used the same threshold value for maximal mitochondrial read fraction for simplicity.
Next, we aligned and merged sequencing libraries by Seurat's canonical correlation analysis or CCA (dimensions: 50) (Butler et al. 2018) using the intersection variable genes across datasets.
Next, principal components were calculated on the variable genes, and the first 50 components were then used to calculate UMAP coordinates. For clustering we used Seurat's implementation of snn/Louvain clustering. Therein, we first calculate the k-nearest-neighbor (knn) graph of cells in PCAspace (dimensions:50). Based on Jaccard similarity scores on the knn graph the shared nearest neighbor (snn) graph is computed. Louvain clustering on the snn graph identified clusters of cells.
Differentially expressed genes were identified by Wilcoxon-test, and filtered for p-values below 0.001, and fold change larger than 2.
We found a group of 1304 interneurons that formed a separate cluster on the very top of the UMAP (Fig 1). These constituted 7.61% of all interneurons and were 94% originating from the Kanton S3 dataset. Both interneuron clusters showed similar expression of classic interneuron markers, and pairwise differential gene expression analysis showed no meaningful differences. Therefore, we lumped these 2 clusters together Integration of organoid and fetal data We obtained raw data for fetal cortical single cell datasets covering age comparable to organoid Individual analysis of Velasco 7 dataset (Fig EV2) The dataset was filtered for high quality cells with a higher gene count than 1000 and analyzed using

Single-cell scoring
We defined specific GO-terms relevant for functional processes in stress and differentiation. Gene lists for each GO-term were downloaded from Ensembl via BiomaRt (Durinck et al. 2009) and intersected with detected genes. Alternative database access is also implemented (see R-package documentation). We then generalized a widely used cell-cycle scoring method based on aggregated gene set activity , and used its implementation in Seurat (AddModuleScore).
Briefly, in the AddModuleScore function the following steps have been implemented as in the original: (A) Take a target set of genes; (B) Calculate their average expression; (C) Create control sets of genes. The control gene-sets are used to control for the cell-to-cell variability in quality and depth.
To create the control sets, first all genes are binned by expression levels (25 bins), then for each gene in the target set, randomly select 100 genes from its expression bin, and finally (D) subtract the average of control from each target gene. The expression binning is important, because expression levels affect the variability of gene expression . Gene lists of 'glycolytic process' (GO:0006096) and 'response to endoplasmic reticulum stress' (GO:0034976) were downloaded and intersected with detected genes, then used to evaluate stress state.

Data partitioning and reclassification
Gene detection in single cells is noisy. To overcome this noise, we grouped cells into small aggregates by high-resolution snn-clustering (as in the manuscript, using the algorithm of Seurat).
Gruffi's aut.res.clustering(), finds an optimal granule resolution with a median of 100-200 cells per granule (cluster). Next, clustering is performed at the determined resolution, resulting in cells separated into 100's of granules, depending on the size of the dataset. Finally, all cells in granules with <30 cells are reassigned to the nearest cluster center in the 3D UMAP space (Euclidean distance) using reassign.small.clusters().

Thresholding and stress annotation
Finally, the average GO-scores for each granule were calculated, and stress level per granule was evaluated. We propose two possible methods to estimate an upper threshold for the assignment of stressed granules for one GO term. (a) Determining manually, based on the expression of stress genes and stress score values on UMAP. Based on these, an empirical quantile (90% if observing 10% stressed cells) as threshold can be assigned. (b) Automatic stress threshold estimation by Shiny. GO.thresh(). For this, we refer to cell number normalized, mean GO scores per granule. In the following this will be referred to as granule score. We assume that: b1) The granules consist of a statistically sufficient number of cells. b2) GO scores of non-stressed cells independently follow the same unknown distribution. b3) GO scores of stressed cells are significantly higher than GO scores of non-stressed cells. b4) The dataset consists of more non-stressed than stressed cells.
Assumptions b1) and b2) together with cell number normalization allow us to use the central limit theorem and hence fit a normal distribution to the granule scores of non-stressed granules. Based on b3) and b4) we conclude that respective non-stressed GO-scores are small, and the mean of the normal distribution can be estimated by the median of granule scores. The standard variation of the normal distribution is now estimated only w.r.t. to GO scores smaller than the median of granule scores. Now the theoretical 99% quantile of the fitted normal distribution can be computed and used as an upper threshold.
Assumption b1) is fulfilled by the automatic clustering resolution search and the reassignment of cell granules with less than 30 cells. Although we based our analysis on thresholds retrieved by this method, since we cannot assure that assumptions b2) to b4) hold true, we highly recommend further inspections and refinements of suggested thresholds in any case. To do so, we propose to visually monitor further manual adjustments via the implemented Shiny App interface.
When considering a combination of GO terms, e.g., response to Endoplasmic Reticulum Stress and Glycolytic Process, we combine the respective thresholds such that a cell is assigned as stressed if either upper threshold is crossed. In case one additionally wants to include a GO term for nonstressed cells, e.g., gliogenesis, the above threshold method can be applied, too, but in this case granules with a score higher than the threshold are assigned as non-stressed. Finally, based on the thresholds on each score, stressed cells are annotated and can be excluded from the dataset.

Other analyses Protein-protein interaction maps
We selected all genes enriched in either stress clusters and jointly ranked them by descending log2FC. We selected the top 150 coding genes, and visualized the 'high-confidence' connected component of the protein protein interaction network using the STRING database (v11.5) (Szklarczyk et al. 2019), links denoting the confidence of connection (permalink: bxso1NJafq8R).

Pathway visualization using ShinyGO and KEGG
The top 150 coding genes (as above) were provided for ShinyGO v0.741 ( (Ge, Jung, and Yao 2020) with default parameters and the background gene list of all 26439 detected genes (from the RNA assay). ShinyGO's visualized enriched KEGG pathways using Pathview and relevant pathways were selected.

Comparison of granular and single-cell scoring
For granular scoring, we used the annotation and approach from (Fig 2). For single-cell scoring, we used the exact same approach, but skipped the granule average calculation of stress scores, and instead we calculated the stress threshold on single-cell scores. For fair comparison, we adjusted the quantile cutoff parameter in the single-cell scoring so that it results in a similar number of stressed cells, as in the granular approach. We then took the symmetric difference of these to find cells only flagged by either, but not both methods. Group median values for were plotted for the four categories (Both, gSC, scSC, Non-stressed) The separation of stressed neurons Stressed cells clearly separated into two major groups, as also seen in (Fig 1). Therefore, we separated Gruffi's classification into two categories. Low resolution clustering (res.0.1.ordered) separated glia (cl.1-2) from neurons (remaining clusters) both better (less clustering artifacts) and simpler than higher resolutions (0.3, 0.5). Intersecting this binary annotation with Gruffi's stress annotation (T, F), separated cells into the four clusters: Neurons, Stressed Neurons, Stressed Progenitors, Progenitors, visible in (Fig S3B).

Progeny pathway activity scoring
Progeny pathways scoring was performed as in vignette, with the following parameters top 200 genes. To visualize the differences between stressed cells identified by typical clustering (Fig. 1B) and Gruffi (Fig. 2F), we separated "Stressed Neurons" and "Stressed Prog." clusters into subsets identified, or not identified by Gruffi, yielding four groups: "Stressed Prog. (Clustering)", "Stressed Prog. (Gruffi)", "Stressed Neurons (Clustering)", "Stressed Neurons (Gruffi)" (Fig S3C). As progeny failed to run on the full object, we randomly downsampled the full dataset to 33.3% of the cells (>50K cells). We visualized the scores using pheatmap with ward.D2 hierarchical clustering and separated the 3 most distinct clusters. For (Fig 3B) we displayed all clusters >2% of all cells (Fig S3D), and all clusters are displayed in (Fig S3E).
Choroid Plexus scoring and stress identification in Samarasinghe et al.
We obtained the Seurat R object of ) from cells.ucsc.edu and performed DGEA by Wilcoxon test in Seurat. We used the clustering presented in the paper, and contrasted "mature choroid plexus" to all other clusters. We calculated a choroid plexus score from the resulting 192 genes (log2fc>1, p.adj<0.01, pct.expr>33%, Supplementary Table S5) and provided this to Gruffi as a negative score (like gliogenesis). We then calculated differential gene expression on stressed cells vs. non-stressed cells, as identified by Gruffi. The resulting 16 genes (log2fc>1, p.adj<0.01) were then analyzed in STRINGdb as before (permalink: bwEKXY7CP0p8).
Neural and glial identity scores for cell type specification As previously , we grouped all neural or progenitor classes to define the respective tow signatures in fetal samples, based on DGE analysis. After intersection with genes that are detected in organoid datasets, the top 30 genes were used for EN, or progenitor signatures (Supplementary Table S6). From these, per-cell subtype scores were calculated using the AddModuleScore() function of Seurat. XY-scatter density plots were drawn by plotting a progenitor score and a neuronal score on the X and Y axis respectively. The cells were colored based on cell type reflecting progenitors, neuronal cells and intermediate progenitors (IPCs), that would physiologically be an intermediate state.
Pseudotime analysis of maturation