comprehensive integrated reference atlases across the human lifespan, in health and disease, together with diverse organoid models in normal and perturbed conditions, will help to clarify the full potential of the human genome. altogether, the heoca will serve as a valuable resource for the organoid research community and a foundation to expand the ability to model human biology. the experiments conducted in this study did not require approval from a specific ethics board. to integrate the atlas, all available datasets were included, with no sample exclusion. for integration method comparisons and sample variance effect analyses, random samples were selected from the full dataset. reproducibility codes for the analyses are available online via github, as detailed in the ‘ code availability ’ section. all statistical methods are described in the corresponding sections of the paper. human intestinal tissue samples were obtained and experimental procedures performed within the framework of the nonprofit foundation htcr ( munich, germany ) including informed patient consent. ileal organoids were derived and maintained according to previously published culture conditions. for the cytokine treatments, organoids were dissociated into five - to ten - cell fragments using tryple ( invitrogen ) and reseeded in matrigel. after 6 days, organoids were treated for 6 days by supplementing the culture medium with 50 ng ml tnf and 200 ng ml rankl ( acro biosystems ) for tnf treatment or with 400 ng ml il - 4 and 40 ng ml il - 13 ( acro biosystems ) for il - 13 and il - 4 treatment. to model host cell responses to viral infection, organoids were treated for 1 day with 1 ng ml ifnα, 1 ng ml ifnγ ( acro biosystems ) and 5 ng ml ifnβ ( peprotech ). to model acute pathogenic inflammation, organoids were treated for 1 day with 10 ng ml tnf, 10 ng ml il - 6, 500 ng ml il - 17a, 1 ng ml ifnγ ( acro biosystems ), 50 ng ml il - 18, 100 ng ml osm ( biolegend ) and 100 ng ml scf ( medchemexpress ). after the indicated treatment durations, organoids were dissociated for scrna - seq using the neural tissue dissociation kit ( p ) ( miltenyi biotec ) as described previously. first, culture medium was removed and organoids were incubated in cell recovery solution ( corning ) for 40 min at 4 °c. next, organoids were transferred to 1 % bovine serum albumin ( bsa ) - coated tubes using hbss – 1 % bsa buffer while pipetting thoroughly to fragmentize the organoids. organoid fragments were centrifuged at 500 g, 5 min, 4 °c. each cell pellet was resuspended in prewarmed buffer x mixed with 25 μl of enzyme p. cells were incubated for 15 min at 37 °c combined with mechanical dissociation by pipetting every 5 min. next, 5 μl of enzyme a in 10 μl of buffer y was added to the digest and incubated for a further 10 min combined with pipetting every 5 min. cells were subsequently washed twice with hbss – 1 % bsa buffer and filtered through a 40 - μm filter coated with 1 % bsa. single cells were counted using a countess 3 fl automate d cell counter ( invitrogen ) and kept on ice. dilutions of ~ 1, 000 cells per μl in 50 – 60 μl of hbss – 1 % bsa buffer were prepared and immediately processed using the 10x chromium next gem single cell 3 ′ reagent kit ( v. 3. 1 ) according to the manufacturer ’ s instructions. libraries were sequenced on illumina ’ s novaseq6000. the scrna - seq data used in this study were obtained from the original papers ( supplementary table 1 ). if the raw fastq files were available, they were downloaded. the seq2science ( v. 1. 2. 2 ) method was used to download the raw fastq files from the gene expression omnibus database ( https: / / www. ncbi. nlm. nih. gov / geo / ) or biostudies database ( https: / / www. ebi. ac. uk / biostudies / ). the reads were aligned to the grch38 genome and ensembl 98 gene annotation using starsolo ( star v. 2. 7. 10b ). in cases in which the raw fastq files were not available, the raw counts were downloaded instead. the downloaded counts and the counts obtained from the realigned reads were merged for subsequent analysis. to integrate the data, we combined the count data from all the samples into a unified dataset. for subsequent analysis, we retained only the genes classified as protein - coding genes and long noncoding rna genes. the low - quality cells in each sample were filed. the raw counts were then normalized to a total count of 10, 000 and log - transformed. given these normalized counts, the top 3, 000 highly variable genes were identified using the default settings in scanpy. these highly variable genes were selected for further downstream analysis. cell - type annotation was performed using the snapseed method ( https: / / github. com / devsystemslab / snapseed ). for each sample, the raw counts were normalized to a total count of 10, 000 and then log - transformed. from these normalized counts, the top 3, 000 highly variable genes were identified using the default settings in scanpy. these highly variable genes were selected as the subset for further downstream analysis. principal component analysis ( pca ) was performed on the normalized data, and the top 30 principal components were chosen for calculating the k nearest neighbors ( k nn ). using the k nn, a umap was generated to visualize the data in a lower - dimensional space. to cluster the data, the leiden clustering method with a resolution of 2 was applied. this clustering approach helped to identify distinc t groups of cells based on their gene expression patterns. previously defined marker genes associated with specific cell types were used to guide the annotation process. to annotate cell types in each cluster, the snapseed method was used. this method calculates the area under the receiver operating characteristic ( roc ) curve ( auc ) and fold change values for each marker gene in relation to the cluster. if multiple markers were available for a particular cell type, the maximum auc and fold change values were selected. the average auc and fold change values were used to represent the specific cell type, and the most specific cell type was annotate d for each cluster based on these criteria. for gene expression level, we merged all counts in each organoid sample by genes using the adpbulk method and applied a natural logarithm transformation to one plus the counts. we then selected the top 500 highly variable genes and calculated pca based on their expression. from principal components 1 and 2 ( pc1 and pc2 ), we selected the top 200 and bottom 200 loading genes for go enrichment analysis using the gseapy method. for the scpoli embedding level, we calculated the mean scpoli embedding for each organoid sample using the adpbulk method, followed by pca based on the mean embedding. a linear model was then used to calculate the covariance between principal components and sample counts, stem cell source, scrna - seq method, tissue type and publication. to benchmark and compare different integration methods, we selected ten random samples from the dataset for validation, repeating this process ten times. twelve integration methods, including pca, seurat ( v. 3, v. 4 and v. 5 ), scvi, scanvi, scpoli, bbknn, harmony, combat, css ( pearson ) and css ( spearman ) were applied to the data to assess their performance in integrating the samples. the scib method, a benchmarking tool, was used to evaluate and compare the results obtained from these integration methods. in the scpoli model, we configured the following parameters for effective training and integration: embedding _ dim was set to 3; hidden _ layer _ sizes were determined as the square root of the total number of cells.