the percentage of cells with expression of a given adt above the positive threshold was calculated for each donor – tissue - group – subset combination. donors with fewer than 50 cells for a particular tissue - group – lineage combination or tissue - group – lineage combinations with fewer than four suitable donors were excluded from analysis. the percent positivity was used as the response variable in the same linear regression model used to detect shifts in composition across age groups. we adjusted p values for multiple comparisons as above, and adts with adjusted p < 0. 05 were considered significant. we constructed donor - level schpf models for each major immune lineage with uniform representation of cells from each donor to identify age - associated gene signatures, as described above. for each schpf model, we performed lmm to account for covariates and identify age associations. each lmm contained six categorical covariates as fixed effects: age group, sequencing chemistry, sex, processing site and cmv serostatus. we also considered three tissue types: mucosal ( bal, lung parenchyma, jlp, jel ), lns ( iln, mln and lln ) and blood - rich, including blood, bm and spleen ), which required us to select one category ( blood - rich ) as a held - out variable. thus, we have two categorical variables for tissue, which effectively represent mucosal versus blood - rich and ln versus blood - rich. we encoded donor identity as a random effect. lmm coefficients and p values were computed for each factor in a given schpf model using the cell scores as response variables by fitting a statsmodels. mixedlm model and using the statsmodels. multipletests function to adjust p values for multiple comparisons ( benjamini – hochberg method ). to cross - validate age - associated schpf factors in other datasets, we further analyzed a bone marrow atlas containing 36 age - annotate d donors with good b cell representation for a b cell aging factor, a lung atlas containing 29 age - annotate d donors with good macrophage and cd8 t cell coverage and pbmc data from the sound life cohort ( age 25 – 65 years, n = 96 ) from the human immune health atlas ( figs. 5 and 6 and extended data fig. 5 ). using the published cell - type annotations from each atlas, we extracted the appropriate scrna - seq profiles and projected them into the corresponding donor - level schpf models generated from the data reported here using the schpf project function. this resulted in cell scores for cells from the external data sets for the same factors that were generated from this data set, allowing us to compare the average cell scores for young versus older donors from the external data. as an orthogonal approach, we also performed pseudobulk de analysis between older and younger donors ( using an age cutoff of 40 years ) from the external data sets, ranked the genes by fc and used gsea to analyze the statistical enrichment of age - associated factors among young versus old donors ( supplementary table 19 ). we used the top 200 genes ( ranked by schpf gene score ) in each age - associated factor as gene sets for gsea. for the de analysis described in fig. 7, we subsetted the mrvi sample embeddings to each tissue group and modeled the predicted ε in mrvi by a linear model adjusting for covariates in sex, cmv serostatus and age group as fixed effects, considering sequencing chemistry and processing site as site covariates in mrvi. a ridge regression parameter of 0. 1, owing to collinearity of cofactors, was added. this decomposition of ε was performed for every single cell. this yields an estimated effect in z - space for each covariate. the effect vector was added to the mean cell embedding in z - space, and degs were computed based on the modified and mean embedding for each cell. for downstream analysis, this matrix of estimated log ( fc ) for each cell and gene was further processed for each immune subset. first, all cells that were represented only in fewer than three samples were filtered out. second, for each cell type, we excluded genes with less than 0. 01 raw average expression or an estimated log ( fc ) across age groups with a 95th percentile below 0. 1, retaining only genes that might be affected by age in a group of cells. to dissect predicted gene effects into modules, neighborhood smoothing was performed using 15 nearest neighbors in u - space and multiplying two times the normalized affinity matrix by the predicted gene effects. spectral co - clustering was performed with four gene clusters and four cell clusters, with mini - batch enabled using sklearn. cluster. marker genes for each module were identified by averaging the predicted log ( fc ) across all cells from the corresponding cell module, and the top 50 genes for each module were identified. we used decoupler - py to compute a module score of log ( fc ) scores using weighted means of the signs of those marker genes ( supplementary table 23 ). for the lung, we isolated a gene module in cd4 t cells that contained t rm cells. to detect similar cells in other tissues, we computed the best cutoff for the module score to identify cells in a specific cell module based on youden ’ s j statistic, computed the module score for all cells from other tissues as described above and applied the same cutoff to all other tissues as the tissue of interest. given that the gut contained t rm cells with a t h 17 phenotype and all other tissues had no module - positive cells, we selected all cells with a mrvi predicted negative log 2 ( fc ) of il17a below −0. 05. to confirm our findings on a per - gene level, we selected module - positive cells and used pseudobulk estimates of de using dreamlet ( supplementary table 24 ). samples with fewer than 5 module - positive cells or 1, 000 total counts and genes with fewer than 3 total counts were removed. aging de was performed using the contrasts method, as described above. genes within a shared functional group were manually selected from the mrvi signature for visualization. pseudobulk de analysis was performed on the classifier - predicted cells in other tissues using the same settings as above in this cell subset. enrichment of module or selected marker genes in the pseudobulk de analysis was performed using gsea implemented in decoupler. run _ gsea ( supplementary table 23 ). cell ranger - mapped tcr and bcr contigs contained in ‘ all _ contigs. fasta ’ and ‘ all _ contig _ annotations. csv ’ output files were re - annotate d using the dandelion preprocessing pipeline. this pipeline includes the following steps: ( 1 ) sample suffix or prefix assignment to each sample barcode; ( 2 ) re - annotation of contigs with igblast ( v1. 19. 0 ) against imgt ( international immunogenetics ) reference sequences ( last downloaded on 24 / 04 / 2023 ); ( 3 ) re - annotation of d and j genes separately using blastn to enable the annotation of contigs without the v gene present; and ( 4 ) identification and recovery of nonoverlapping individual j gene segments. for bcrs, three additional steps were also performed: ( 1 ) additional re - annotation of heavy - chain constant ( c ) region calls using blastn ( v2. 13. 0 ) against curate d sequences from ch1 regions of respective isotype class; ( 2 ) heavy - chain v gene allele correction using tigger ( v1. 0. 1 ); and ( 3 ) bcr mutation calling. cell - level quality control was performed using dandelion ’ s ‘ filter _ contigs ’ function, which only considers productive vdj contigs, asserts that a single cell should only have one vdj and one vj pair or only an orphan vdj chain and explicitly removes contigs that fail these checks ( except for igm / igd and trb / trd extra pairs ). contigs that did not match any cell barcodes in the gene expression data were also removed at this step. tcrs and bcrs were then grouped into clones or clonotypes.