the plots showing the expression of certain genes or proteins were created based on the normalized values ( denoised values ) using a quantile cutoff of 0. 99 and via either the ggplot2 r package ( v3. 2. 1 ) or the scanpy. pl. umap function of the scanpy python package ( v1. 5. 1 ). to find the conserved human and mouse kc markers we started by identifying the human kc markers. we mapped the annotation of the human myeloid umap on the human pool of techniques / patients umap to identify the real kcs in this last umap. the real kcs were identified as the top part of the mac cluster. using this new annotation we then calculated the de genes and de proteins for each cluster. some genes are listed as marker for multiple clusters, only for the cluster where the gene had the highest score ( raw _ normalized _ mean1 / raw _ normalized _ mean2 lfc _ mean ), the gene was kept as marker. this way we found 110 potential human kc markers. we then created a heatmap of these 110 genes ( using denoised gene values scaled between 0 and 1 ) and filtered this heatmap by removing the genes where the scaled normalized value was higher than 0. 50 in more than 30 % of the cells of a certain cell type other than kcs. except for the macs, we only removed a gene when it had a scaled normalized value higher than 0. 50 in more than 70 % of the macs. after this filtering we ended up with 36 human kc markers. next we converted these human gene symbols into mgi ids via the biomart tool on the hgnc website ( https: / / biomart. genenames. org / martform / # ! / default / hgnc ? datasets = hgnc _ gene _ mart ). we found a mgi id for 30 genes. we then converted these mgi ids into mouse gene symbols via the mgi webtool ( http: / / www. informatics. jax. org / batch / ). to identify the mouse kc markers we similarly mapped the annotation of the mouse myeloid umap on the mouse pool of techniques umap to identify the real kcs in this last umap. the real kcs matched with the mac cluster. similarly as in human, the de genes for each cluster was calculated and genes listed as marker for multiple clusters were dealt with in a similar way. this way we found 264 potential mouse kc markers. we then removed the genes that had a score ( raw _ normalized _ mean1 / raw _ normalized _ mean2 lfc _ mean ) lower than 10 and ended up with 214 genes. we then created a heatmap of these 214 genes ( using denoised gene values scaled between 0 and 1 ) and filtered this heatmap by removing the genes where the scaled normalized value was higher than 0. 50 in more than 30 % of the cells of a certain cell type other than kcs. after this filtering we ended up with 68 mouse kc markers. next we converted these mouse gene symbols into mgi ids via the mgi webtool ( http: / / www. informatics. jax. org / batch / ). we then converted these mgi ids into human gene symbols via the biomart tool on the hgnc website ( https: / / biomart. genenames. org / martform / # ! / default / hgnc ? datasets = hgnc _ gene _ mart ) and ended up with 60 genes. at this point we found 30 human kc markers and 60 mouse kc markers. in a next step, we only kept the human kc markers that we identified as a highly variable gene ( hvg ) in the mouse pool of techniques umap ( 20 genes ) and the mouse kc markers that were identified as hvgs in the human pool of the techniques umap ( 30 genes ). we next put these 20 mouse kc markers in singlecellsignatureexplorer ( pont et al., 2019 ) to see where these genes are enriched in the mouse pool of techniques umap. in order to only get an enrichment in the kcs we decided to only use top 10 mouse kc markers ( ordered on score ), together with slc40a1 and hmox1. we then started to add the top human kc markers as long as we keep the enrichment solely in the kcs. this way we ended up with final list of 15 human - mouse conserved kc markers. we next converted these kc markers into the monkey, pig, chicken or zebrafish orthologs by looking up the human gene symbol on ncbi ( https: / / www. ncbi. nlm. nih. gov / search / ) and checking if there is an ortholog of the species of interest listed under the ‘ ortholog ’ tab. the found orthologs were then used as input for the singlecellsignatureexplorer tool. the protein normalized values ( denoised values; calculated in the totalvi workflow ) were converted into an fcs file using the write. fcs function of the flowcore r package ( v1. 50. 0 ). we first removed per sample all spots that were clear outliers compared to the location of the tissue. each sample was then normalized individually using the sctransform function of the seurat r package ( v3. 2. 3 ) with default parameters. all samples were then merged with the merge function of the seurat r package ( v3. 2. 3 ) with default parameters. next, we determined the hvgs, created a pca plot, performed clustering and created an umap plot as described in the spatial workflow available on the seurat website ( https: / / satijalab. org / seurat / articles / spatial _ vignette. html ). clusters which showed high mitochondrial gene expression were removed. spots located at the darker parts of the tissue were also removed as these parts are considered to be dead tissue or of bad quality. for modelling the cell type composition and zonation, spatial cite - seq and transcriptomics data were analyzed using probabilistic graphical models, similar to what is used in tools such as cell2location and scvi. in brief, transcriptomics data was modelled as a negativebinomial distribution, parameterized with a mean and dispersion, the latter optimized as a free parameter for each gene. visium highly multiplexed protein data was modelled as a mixture of negativebinomials, with a and and a shared dispersion. the actual foreground / background signal within a modality was modelled as a that depends on the latent space, and which is multiplied with the empirical library size to get. for visium highly multiplexed protein, was modelled as a latent variable specific for each gene. for visium highly multiplexed protein and for rna - seq were modelled as deterministic functions depending on the use case as described in the following paragraphs. the posterior of the probabilistic graphical model was inferred using black - box variational inference ( ranganath et al., 2013 ), in which the variational distribution was specified as a diag onal normal distribution, transformed into the correct domain using transforms (,, ). free parameters within this model were optimized using gra dien t descent, with the elbo as loss function and adam as optimizer as implemented in pytorch ( paszke et al., 2019 ) ( pytorch. org ). we used a learning rate of 0. 01 for variational parameters, and 0. 001 for parameters of the amortization functions. to calculate the average expression of each gene within a cell type, we used a linear model in which both and were modelled as a latent variable specific for each gene and cell type. the for nuclei were multiplied with a gene - specific correction factor ( optimized as a latent variable ) that corrected for differences between scrna - seq and snrna - seq. given that spatial transcriptomics data sequences the whole cell, the uncorrected values were used for spatial deconvolution. to infer the proportions of each cell type within a spot, we used a model in which the gene expression is modelled as a linear combination of cell type proportions and average expression in each cell type: source paper: pmc8809252 for we adapted the values from the reference, but included - a capture bias per gene, which corrects for technical and biological differences between spatial and sc / sn - rna - seq. the capture bias was modelled as a latent variable with prior - a red blood cell cell type, which was not included in the reference dataset but nonetheless had a dominant presence in the spatial data.