the umap embedding was transformed using the primary tissue atlas umap model. for each cell, the system selected its 100 nearest neighbors from the primary tissue dataset. the predicted cell type for the cell was assigned based on the tissue that was most frequently observed among its 100 nearest neighbors. to identify degs in primary tissue stem cells and enterocytes, we subsetted these cell types and used a linear model to calculate the covariance between sample age and gene expression for each gene. the top 100 genes with the highest coefficients were selected as degs. the gseapy method was then applied to identify the top go - enriched terms associated with these genes. to identify the heterogeneity of intestinal organoid stem cells and enterocytes, cells from the hioca were subsetted. integration was performed using the css method based on 1, 000 highly variable genes across all cells. leiden clustering with a resolution of 0. 1 was applied to identify subclusters. the wilcoxon rank - sum test was used to identify degs among subclusters, and go enrichment analysis was conducted on the top 500 degs of each group using gseapy. the scrna - seq data from both duodenum fetal and adult primary tissues were obtained from two research papers. the top 3, 000 highly variable genes were subsetted for data integration. to integrate all the cells, we applied the scpoli method with the same parameters used in the heoca atlas integration. the scpoli model was saved for the downstream comparison. for each organoid sample, the same set of variable genes used in the primary tissue atlas was chosen, and the scpoli query was executed using identical parameters to those used in the primary tissue atlas training model. the umap embedding was transformed using the primary tissue atlas umap model. for each cell, the system selected its 100 nearest neighbors from the primary tissue dataset. the predicted cell type for the cell was assigned based on the tissue that was most frequently observed among its 100 nearest neighbors. samples of scrna - seq raw reads were mapped to the human genome, and counts of the matrix were obtained. the same set of variable genes used in heoca was chosen, and the scpoli query was executed using identical parameters to those used in the heoca training model. the umap embedding was transformed using the heoca umap model. for each cell, the system selected its 100 nearest neighbors from the heoca dataset. the predicted cell type for the cell was determined by assigning it the cell type that was most frequently observed among its 100 nearest neighbors at the level 2 cell - type classification. similarly, the predicted tissue for the cell was assigned based on the tissue that was most frequently observed among its 100 nearest neighbors. for each cell in the organoid protocols, organoid perturbation, and disease samples, a matched heoca cell was reconstructed using the top ten k nn in heoca. the mean expression of these ten neighbors was calculated to represent the expression profile of the matched sample reference in heoca. in addition, the mean k nn distance of these ten neighbors was used to represent the cell ’ s distance to the heoca. to compare expression levels of the samples, the above - mentioned matched sample reference in heoca was identified. the expression difference per gene for each cell pair was calculated based on the log - normalized expression values. for each gene, the variance over the calculated expression difference per cell pair was compared with the sum of squared expression differences normalized by the number of cell pairs. an f test was applied to test for differential expression for each gene. scrna - seq reads for each sample were mapped to the human genome, and gene counts were generated using cellranger. these counts served as input for sc2heoca, with default settings used to map all samples to heoca. during mapping, each cell was annotate d with a level 2 cell type, and the distance to heoca was calculated. cell proportions were determined based on the mapping annotations. for the raw integration of perturbation samples, three samples were merged, highly variable genes were identified using scanpy with default settings, and the samples were integrated using the combat method. the sc2heoca package with default settings was used to identify degs, and go enrichment analysis was performed using gseapy on the degs of each group. degs between each treatment sample and control sample were calculated using the wilcoxon rank - sum test in scanpy ( v. 1. 9. 3 ). the scibd database was downloaded, and samples from healthy individuals, and patients with colitis, crohn ’ s disease and ulcerative colitis were extracted. only epithelial cells were selected for downstream analysis. pseudo - bulk gene expression was calculated for each individual, and the degs identified in the previous step were subsetted. the mean expression of these genes across patients was used to compare gene expression between inflammatory and viral response conditions. the disease sample analysis is similar to the sample incorporation step. the raw count matrices were downloaded from the original papers. the same set of variable genes used in heoca was chosen, and the scpoli query was executed using identical parameters to those used in the heoca training model. the umap embedding was transformed using the heoca umap model. for each cell, the system selected its ten nearest neighbors from the heoca dataset. the predicted cell type for the cell was determined by assigning it the cell type that was most frequently observed among its ten nearest neighbors at the level 2 cell - type classification. the predicted tissue for the cell was assigned based on the tissue that was most frequently observed among its ten nearest neighbors. for each cell, the mean distance of its ten nearest neighbors was assigned as its mean distance to heoca. in the analysis of degs between colon cancer organoid colonocytes and bronchial copd organoid basal cells, we performed separate subsetting for all colonocytes and basal cells. for each dataset, we isolated the top 3, 000 highly variable genes. we then integrated these subsets of cells using the bbknn method. to cluster the two datasets, we applied the leiden method with resolutions of 1 and 2 in two datasets. the clusters predominantly associated with the disease were selected as disease state cells, while the remaining clusters were categorized as normal state cells. we used d2c to score the expression levels of 2, 395 drug target signatures in single cells from the human organoid cell atlas and the human lung cell atlas. d2c scores were scaled to mitigate scale differences between different datasets in the atlas. we used the r package scdecaf ( v. 0. 99. 0 ) to select drug target signatures that exhibited global covariation in one or more cell types in heoca and hlca primary tissue atlases. the inputs to scdecaf were the scaled d2c z - scores and the cell embeddings from the atlases. the shrinkage operator in scdecaf was set to lambda = exp ( −1. 3 ) based on reconstruction error plots made available in the scdecaf package. we assigned a drug signature to a cell type if more than 50 % of the cells from the cell type had a signature score above median across all cell types. multicellular drug target signatures were identified whether a drug signature was selected in at least two cell types. to assess druggability potential of organoid cell types, we computed the cosine similarity for cell - type pairs in organoid and primary tissue based on multicellular drug signatures identified in primary tissue and organoid models. further information on research design is available in the nature portfol io reporting summary linked to this article.