a second latent space, z, is defined by adding an attention - based concatenation between u and the sample embedding space to the original u - space. another layer of attention is used to incorporate an embedding of 10× genomics chemistry and experimental site ( cambridge, uk versus columbia, ny ), and this third latent space is decoded using a linear decoder to yield the rate of a negative binomial distribution. we use a cell - type - aware gaussian mixture prior in u - space. to introduce cell type awareness, we use a bias to the mixture proportions that makes it likely for cells of the same type to be sampled from the same gaussian. for the latent embedding highlighted throughout the article and used for manual cell - type curation, we used the donor identities as the sample keys and used the output of mmochi classification ( see above ) as the cell - type prior in mrvi. we used default parameters except n _ epochs _ kl _ warmup of 25, n _ latent _ u of 20, n _ latent in z - space of 200, dropout in qz as well as pz of 0. 03 ( adopted from a previous publication ). to visualize cells ( either the total immune component or individual major lineages ), we computed nearest neighbors ( scanpy. pp. neighbors ) on the mrvi u latent space and calculated umap embeddings ( scanpy. tl. umap ) using the 15 nearest neighbors, a minimum distance of 0. 4, a spread of 1. 0 and initialization with paga after running scanpy. tl. paga. to identify additional heterogeneity in cell states within samples in addition to the cell - type annotation provided by mmochi, we performed manual annotation. for each mmochi annotate d population, a new scvi model was trained with donor as the batch key, then leiden clustering ( scanpy. tl. leiden ) was performed on the lineage - specific neighbors graph at an appropriate resolution, selected to minimize over - clustering ( ranging from 1 to 15 ). markers for each cluster were computed by scanpy. tl. rank _ genes _ groups, and clusters with similar marker expression were merged. to annotate proliferating cells, scanpy. tl. score _ genes _ cell _ cycle was run, and the output was used in combination with the gene expression of mki67 and top2a. we focused our de analysis on immune lineages and cell types with sufficient representation across experimental sites, tissues and donor ages. this included six tissue groups: blood, bm, spleen, gut ( jlp and jel ), lns ( iln, lln and mln ) and lungs ( consisting of bal and parenchyma ); six immune lineages: myeloid, cd4 t cells, cd8 t cells, invariant t cells ( that is, γδ t cells and mait cells ), b cells and ilc / nk cells; and 26 individual cell types within all lineages. covariates included 10× genomics chemistry ( 3 ′ versus 5 ′ ), sex ( male versus female ), laboratory ( cambridge, uk versus columbia, ny ) and cmv status ( positive versus negative ). for aging analyses, donors were categorized as being < 40 or > 40 years of age. variance decomposition and pseudobulk de analysis were performed using lmm through the dreamlet r package ( v1. 4. 1 ). depending on the resolution of the analysis, de was performed separately either for each immune lineage ( for example, myeloid cells, b cells, and so on ) or for each immune subset ( for example, macrophages, naive b cells and so on ) using the cluster _ id parameter in dreamlet. the raw gex count matrix was pseudobulked across samples, and each tissue in each donor was treated as a separate sample. before performing de, samples and genes with poor representation were filtered using dreamlet:: processassays. samples with fewer than 50 cells and genes not represented in at least 40 % of the samples with at least five counts were excluded. to confirm findings by mrvi counterfactual analysis ( see below ), these thresholds were reduced to a minimum of ten cells for a sample to be included, and at least 10 % of samples with at least five counts. de for a subset was not performed when fewer than three or four samples ( for tissue and age analysis, respectively ) met the minimum cell thresholds. variance decomposition was performed for age analysis for each lineage using dreamlet:: fitvarpart with sex, sequencing chemistry, cmv serostatus, age group, processing site and tissue as covariates ( supplementary table 13 ). lmm was performed using dreamet:: dreamlet with ebayes estimation enabled. tissue effects ( figs. 3 and 4, extended data figs. 3 and 4, and supplementary tables 4, 7 and 10 ) were modeled by comparing each lineage / subset in one tissue group to the same lineage / subset in the remaining tissue groups, with donor identity encoded as a random effect. age effects ( figs. 5 – 7, extended data fig. 5, and supplementary tables 15 and 24 ) were modeled across each tissue - group – age - group combination while controlling for cmv serostatus and sex as fixed effects and with sequencing chemistry and processing site as random effects. age effects within each tissue group were then measured using the contrasts parameter in dreamlet:: dreamlet between old and young for each tissue group ( for example, the effect of age in the gut was computed as ‘ old - gut − young - gut ’ ). cmv effects ( supplementary fig. 10 and supplementary table 21 ) were modeled across each tissue - group – cmv serostatus combination while controlling for age and sex as fixed effects and with sequencing chemistry and processing site as random effects. cmv effects within each tissue group were then measured using the contrasts parameter in dreamlet:: dreamlet between cmv and cmv for each tissue group. for identifying cross - tissue and cross - donor gene signatures for each major immune lineage, we constructed probabilistic factor models directly from scrna - seq count matrices using schpf. the output of schpf includes two matrices: an m × k gene score matrix containing weights for each of m genes in each of k factors and a k × n cell score matrix containing weights for each of n cells in each of k factors. in the original report of schpf, the algorithm required a user - supplied value of k, the number of factors in the model. here, we use a new consensus factorization implementation of schpf, in which the user specifies a broad range of k values from which many schpf models are generated. the gene score matrices for these models are then clustered to identify k recurrent factors, which are combined to seed a final round of training to construct a final consensus model with k factors. we constructed two types of schpf models: a tissue - level model ( extended data fig. 4 ), in which the number of cells from each of three tissue groups was balanced by random sub - sampling ( gut: jel and jlp; lung parenchyma; and lns: mln and lln ), and a donor - level model ( figs. 5 and 6 ), in which the number of cells from each organ donor was balanced. we constructed both types of models for the major immune lineages: cd4 t cells, cd8 t cells ( including all invariant t cells ), nk cells, ilcs, b cells and macrophages. for donor models, donors with fewer than 300 cells for a given lineage were removed. in both models, the count matrices were randomly downsampled such that the average number of transcripts per cell was the same for each donor to avoid coverage bias. schpf models considered only protein - coding genes ( excluding tcr and immunoglobulin cassettes ) detected in at least 1 % of cells across the final subsampled and downsampled training matrix. for all consensus schpf models, we ran schpf five times for each of 16 values of k ( 15 – 30 ), from which we selected the top three models for each value of k based on convergence criteria for clustering. we applied walktrap clustering to identify recurrent clusters, which we required to form clusters with factors from at least two different models from which we trained the final consensus model. immune subset composition within each lineage across tissues was visualized by violin plots or box plots ( using seaborn ).