tissue - specific enrichment of immune subset frequencies in specific tissues was also assayed within each major lineage using sccoda for bayesian inference. significant enrichment of an immune subset in one tissue over the rest was determined using sccoda. util. comp _ ana. compositionalanalysis to detect credible effects, and was run sequentially, selecting each cell type as the reference. majority voting was then used to identify cell types that are credibly changing more than half the time with automatic reference - subset selection and the default false discovery rate of 0. 05. to determine tissue - specific gene expression signatures across immune lineages ( fig. 3 ), significant degs were defined as adjusted p < 0. 05 and log 2 ( fc ) > 1 by pseudobulk de across tissues at the lineage level ( see above ). mean z - score gene expression was calculated for each pairing of tissue group and lineage. genes and samples were both hierarchically clustered using scipy. cluster. hierarchy. linkage with ward ’ s method and euclidean distance. discrete clusters of genes with similar expression patterns were calculated using scipy. cluster. hierarchy. fcluster with the ‘ maxclust ’ method ( supplementary table 6 ). for each gene cluster, association with specific tissue groups or lineages could arise from de within one or more specific subsets of that lineage or from compositional shifts in subsets across tissues. to disentangle these possibilities, we first used pre - ranked gsea to compare the gene clusters identified via lineage - level de to the effect size ( that is, log ( fc ) ) of de across tissues in the subset - level de. to visualize potential effects caused by compositional shifts across tissues, we computed the average frequency of the subset ( as a proportion of the total cells within that lineage group ) within a tissue, the fc of that frequency over the frequency in the remaining tissue groups and the average expression of the gene cluster. to assess whether differential transcript expression was reflected in the surface protein profiling ( supplementary tables 11 and 12 ), we selected adts corresponding to degs in at least one tissue. to identify enrichment in one tissue group over the other tissue groups, we used scanpy. tl. rank _ genes _ groups on the normalized expression with wilcoxon and tie - correction enabled. to minimize the influence of technical staining artifacts or donor covariates, analysis was conducted separately within each donor. 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. before de analysis, the adt count matrix was subsampled to equalize cell numbers and randomly downsampled such that the average number of transcripts per cell was the same for each group to avoid coverage bias. we next sought to identify factors from the tissue - level schpf models of each major immune lineage that were shared across cell types. as described above, we first constructed consensus schpf models for cd4 t cells, cd8 t cells, macrophages, nk cells, ilcs and b cells with equal representation of cells from each of three major tissue groups ( gut, lung and lns ). from each model, we removed probable nuisance factors containing heat shock protein - encoding genes ( common dissociation artifact, > 1 gene ), ribosomal protein - encoding genes ( common coverage artifact, > 10 genes ), genes from the highly inducible metallothionein cluster ( > 1 gene ), hemoglobin transcripts ( red blood cell contamination, > 0 genes ) and genes in a previously published signature of dissociation - induced cell stress in scrna - seq ( > 7 genes ) among the 30 top - weighted genes. next, we computed the average cell score for each factor in each of the three major tissue groups and identified all factors with an average tissue - group score that was at least 80 % higher in one tissue group than the average of the remaining two. thus, the resulting set of 53 schpf factors from across all 6 lineage - specific models exhibits some degree of tissue specificity. to compare these factors to each other, we computed the pearson correlation between the gene score vectors for each pair of factors. we then identified factors with a pairwise correlation that was greater than the 95 % confidence threshold with at least two other factors, which yielded 31 schpf factors from across the six major immune lineages. finally, we performed hierarchical clustering of the pearson correlation matrix for these 31 factors ( seaborn. clustermap using euclidean distances ) to identify modules containing factors with similar gene signatures that originated from different, lineage - specific schpf models ( extended data fig. 4 ). modules of genes were further interrogate d by average gene expression and validated in specific immune subsets using pseudobulk gex de and adt de as described above. to detect shifts in the subset composition of specific lineages across the age groups and cmv serostatus ( supplementary table 14 ), we performed generalized linear modeling by fitting a statsmodels. glm model for each tissue subset, considering sex, sequencing chemistry, cmv serostatus and processing site as additional covariates. 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 estimated coefficients were used to calculate a covariate - aware log 2 ( fc ) for visualization. we used the statsmodels. multipletests function to adjust p values for multiple comparisons ( benjamini – hochberg method ), and subset – tissue combinations with adjusted p < 0. 05 were considered significantly changing across age. to depict age - associated effects on the immune system, we visualized the similarity of trending degs by age on immune subsets across tissues using t - distributed stochastic neighbor embedding ( t - sne ) ( fig. 5d ). we first calculated trending deg ( unadjusted p < 0. 05; < 40 or > 40 years old ) pairwise similarities by summing the intersection of positively regulated genes ( log 2 ( fc ) > 0. 1 ) and negatively regulated genes ( log 2 ( fc ) < −0. 1 ), divided by the overall union of both. the similarity or distance ( 1 − similarity ) was applied to cell types containing more than 70 degs ( unadjusted p < 0. 05; mean log - normalized expression, > 0. 05 ) and present in at least 3 donors per tissue and age group. the similarity levels of cell types and tissues with more than 200 degs were further clustered using the ward. d2 method and projected into a distance - based t - sne illustration. to investigate the effect of age on specific genes within each immune subset, we plotted genes that were significant in at least one tissue ( adjusted p < 0. 05 ) and within the top 50 significant genes. although our power to detect age effects by de was limited, genes that were significantly de in one subset were often trending in the same direction across multiple tissue groups. to assess the effect of age on surface protein expression ( supplementary table 17 ), we used landmark - registered protein expression data ( by mmochi, see above ) to account for donor - to - donor batch effects in adt staining quality. although landmark registration preserves the separation between positive - expressing and negative - expressing cells for thresholding, this non - parametric normalization can obscure changes in overall expression intensity between samples. therefore, we focused on shifts in percent positivity for a marker in each tissue subset. we performed automatic thresholding by mmochi, followed by manual adjustment ( as described above ), on the landmark - registered expression of all adts corresponding to a deg by age.