Single-cell transcriptomics of human T cells reveals tissue and activation signatures in health and disease

1T lymphocytes coordinate adaptive immune responses and are essential for establishing protective immunity and maintaining immune homeostasis. Activation of naive T cells through the antigen-specific T cell receptor (TCR) initiates transcriptional programs that drive differentiation of lineage-specific effector functions; CD4+ T cells secrete cytokines to recruit and activate other immune cells while CD8+ T cells acquire cytotoxic functions to directly kill infected or tumor cells. Most of these effector cells are short-lived, although some develop into long-lived memory T cells which persist as circulating central (TCM) and effector-memory (TEM) subsets, and non-circulating tissue resident memory T cells (TRM) in diverse lymphoid and non-lymphoid sites–. Recent studies in mouse models have established an important role for CD4+ and CD8+ TRM in mediating protective immunity to diverse pathogens,–. Defining how tissue site impacts T cell function is therefore important for targeting T cell immunity.

2In humans, most of our knowledge of T cell activation and function derives from the sampling of peripheral blood. Recent studies in human tissues have revealed that the majority of human T cells are localized in lymphoid, mucosal and barrier tissues and that T cell subset composition is a function of the specific tissue site,. Human TRM cells can be defined based on their phenotypic homology to mouse TRM and are distinguished from circulating T cells in blood and tissues by a core transcriptional and protein signature–. However, the role of tissue site in determining T cell functional responses, and a deeper understanding of the relationship between blood and tissue T cells beyond composition differences are key unanswered questions in human immunology.

3The functional responses of T cells following antigen or pathogen exposure have been largely defined in mouse models, and are generally classified based on whether or not they secrete specific cytokines or effector molecules. Effector CD4+ T cells comprise different functional subtypes (Th1 cells secrete IFN-γ and IL-2; Th2 secrete IL-4, 13; Th17 secrete IL-17A, etc.), while effector CD8+ T cells secrete pro-inflammatory cytokines (IFN-γ,TNF-α) and/or cytotoxic mediators (perforin and granzymes). Certain conditions can lead to inhibition of functional responses; for example, CD4+ T cells encountering self-antigen become anergic and fail to produce IL-2, while CD8+ T cells responding to chronic infection, tumors, or lacking CD4+ T cell help become functionally exhausted, and express multiple inhibitory molecules (e.g., PD-1, LAG3)–. While human T cells can produce similar cytokines, effector and inhibitory molecules as mouse counterparts–, the full complement of functional responses for human T cells in tissues has not been elucidated. Thus, establishing a baseline of healthy T cell states in humans is essential for defining dysregulated and pathological functions of T cells in disease.

4Single cell transcriptome profiling (scRNA-seq) has enabled high resolution mapping of cellular heterogeneity, development, and activation states in diverse systems,. This approach has been applied to analyze human T cells in diseased tissues, and in response to immunotherapies in cancer; however, baseline functional profiles of human T cells in healthy blood and tissues would be an important reference dataset. We have established a tissue resource where we obtain multiple lymphoid, mucosal, and other peripheral tissue sites from human organ donors–,,,, enabling study of T cells across different anatomical spaces.

5Here, we used scRNA-seq of over 50,000 resting and activated T cells from lung (LG), lymph nodes (LN), bone marrow (BM) and blood, along with integrated computational analysis to define cellular states of homeostasis and activation of human blood and tissue-derived T cells. We reveal how human T cells in tissues relate to those in blood, and identify a conserved tissue signature and activation states for human CD4+ and CD8+ T cells conserved across all sites. We further show how scRNA-seq profiles of T cells associated with human tumors can be projected onto this healthy baseline dataset, revealing their functional state. Our results establish a high dimensional reference map of human T cell homeostasis and function in multiple sites, from which to define the origin, composition and function of T cells in disease.

6We obtained BM, LN, and LG as representative primary lymphoid, secondary lymphoid and mucosal tissue sites, respectively, from two deceased adult organ donors who met the criteria of health for donation of physiologically healthy tissues for lifesaving transplantation, being free of chronic disease and cancer (Supplementary Table 1). For comparison, we obtained blood from two healthy adult volunteers. CD3+ T cells isolated from tissues and blood were cultured in media alone (“resting”) or in the presence of anti-CD3/anti-CD28 antibodies (“activated”) (Fig. 1a). Single cells were encapsulated for cDNA synthesis and barcoded using the 10x Genomics Chromium system, followed by library construction, sequencing, and computational identification of T cells (Supplementary Fig. 1, Supplementary Table 2, Supplementary Data 1).

7We initially analyzed tissue T cell populations from the two individual donors, comprising six samples per donor (resting and activated samples from three tissue sites). We merged all data for each donor, performed unsupervised community detection to cluster the data based on highly variable genes (Supplementary Data 2), and projected cells in two dimensions using Uniform Manifold Approximation and Projection (UMAP). For both donors, the dominant sources of variation between cells were activation state (vertical axis) and CD4/CD8 lineage (horizontal axis) (Fig. 1b). Tissue site was also a source of variability; T cells from BM and LN co-localized while LG T cells were more distinct (Fig. 1b), consistent with a greater proportion of CD8+ T cells and TRM phenotype cells in LG relative to the two lymphoid sites (Supplementary Fig. 2 and previous studies,,).

8Differential gene expression from the scRNA-seq data resolved T cell subsets and functional states within and between sites and lineages into 10–11 clusters (Fig. 1c, Supplementary Data 3, 4). CD4+ T cells comprised 6–7 clusters: resting cells expressing CCR7, SELL and TCF7, (corresponding to naive or TCM cells); three activation-associated clusters expressing IL2, TNF, and IL4R at different levels; TRM-like resting and activated clusters expressing canonical TRM markers CXCR6 and ITGA1,; and a distinct regulatory T cell (Treg) cluster expressing Treg-defining genes FOXP3, IL2RA, and CTLA4 (Fig. 1c). CD8+ T cells comprised four clusters distinct from CD4+ T cells and included: two TEM/TRM-like clusters expressing CCL5, cytotoxicity-associated genes (GZMB, GZMK), and TRM markers (CXCR6, ITGA1); an activated TRM/TEM cluster expressing IFNG, CCL4, CCL3; and clusters representing terminally differentiated effector cells (TEMRA) expressing cytotoxic markers PRF1 and NKG7 (Fig. 1c). In terms of tissue distribution, TRM cells were largely in the lung, Tregs were primarily identified in LN, while TEMRA cells were enriched in BM (consistent with phenotype analysis, Supplementary Fig. 2); the remaining resting and activated CD4+ and CD8+ T cell clusters derived from all sites (Fig. 1b, c). These results show subset-specific profiles in human tissues, but suggest similar activation profiles across sites.

9To assess how blood T cells relate to those in tissue, we performed scRNA-seq analysis of resting and activated blood T cells from two adult donors, and projected the merged data onto the UMAP embeddings of T cells from each tissue donor (Fig. 2a, b). The majority of blood T cells co-localized with resting or activated T cells from BM but did not exhibit substantial overlap with LG or LN T cells from either donor, particularly in the resting state (Fig. 2a, b). We also quantified the number of blood T cells that were transcriptionally similar to CD4+ and CD8+ T cells from each tissue within resting or activated samples (Fig. 2c, d). Resting blood T cells were highly represented among CD4+ and CD8+ T cells in BM (Fig. 2c, d). Interestingly, a substantial number of unstimulated blood T cells projected onto activated CD4+ T cells in BM for both donors (Fig. 2c, d, left panels). In contrast, activated blood T cells were strongly represented among activated CD4+ T cells for all tissue sites and in LN for CD8+ T cells (Fig. 2c, d; right panels). Similar results were obtained when each blood sample was compared separately to each tissue donor (Supplementary Fig. 3), and when blood T cells were projected onto tissue T cells using scmap, an alternative scRNA-seq data projection package (Supplementary Fig. 4). These results indicate that resting blood T cells are most similar to those in the BM, while activated blood and tissue-derived T cells share common signatures.

10The major transcriptional differences between tissue and blood T cells based on population-level RNAseq originate from the presence of TRM in tissues. Because scRNA-seq enables high-resolution detection of gene expression differences that can be unambiguously traced to individual T cells, we investigated whether there were intrinsic features of tissue T cells that distinguished them from blood. Resting memory T cells in tissues and blood express high levels of CCL5 (Supplementary Fig. 5), a marker of CD8+ TEM cells, enabling direct comparison of gene expression between similar subsets. We identified a similar complement of genes that were highly expressed in TEM cells from each tissue compared to blood (Fig. 3a–c). Interestingly, these tissue-intrinsic genes include those associated with microtubules and cytoskeleton (tubulin-encoding genes TUBA1A, TUBA1B, TUBB, TUBB4B; S100A4) as well as genes encoding cell matrix, membrane scaffolding, and adhesion molecules (VIM or vimentin, galectins LGALS1/LGALS3, AMICA1, ITM2C, EZR, annexins ANXA1/ANXA2) (Fig. 3a–c). TRM signature genes including ITGA1 and ITGAE were also upregulated in tissues compared to blood, particularly in the lung (Fig. 3a–c). These findings suggest that localization of T cells in tissues likely involves structural changes in the cell that facilitate interactions with tissue matrix.

11We next compared the single-cell distribution of average expression of tissue signature genes in the blood and tissues (Fig. 3d). CCL5+ TEM cells from all three tissues (both donors) express higher levels of tissue signature genes compared to blood, though LG and LN T cells have higher expression than those from BM (Fig. 3d). Notably, a minute fraction of blood TEM cells (<0.5%) express this tissue signature at levels comparable to that in LN (within one standard deviation of the mean for all tissues). Shown in a heat map are the relative expression levels for genes within the tissue signature, including genes enriched in human TRM cells,, and genes associated with cytoskeletal, cell-matrix interactions, cell division, apoptotis, and signaling (Fig. 3e). Expression of the tissue signature genes is highest in LG, followed by LN and BM expressing only a subset of tissue-associated genes; the outlier subpopulation from blood expresses a fraction (<40%) of tissue signature genes at levels comparable to those in tissues (Fig. 3e). When resting TEM cells from all sites and donors were visualized by UMAP using the tissue-associated signature genes, blood T cells clustered distinctly from all tissues, while LG T cells clustered distinctly from LN and BM (Fig. 3f). Notably, a subset of T cells from BM and LG clustered more closely to blood T cells (Fig. 3f), indicating the presence of circulating T cells within these sites. Together, these results show that tissue T cells express genes associated with infiltration and localization in tissues along with residency markers, while blood contains only trace numbers of cells expressing these genes.

12The tissue signature identified in Fig. 3 compared tissue from deceased organ donors to blood from living individuals. To establish that the observed differences were not due to tissue processing and/or T cells from organ donors versus living individuals, we analyzed our tissue signature in scRNAseq data from several available datasets, including BM from living individuals and additional blood data (see Methods). We found that the tissue signature was significantly enriched in BM from living individuals, compared to blood (Supplementary Fig. 6a). Similarly, we found the tissue signature was enriched in all organ donor sites compared to blood from additional living donors (Supplementary Fig. 6b). Together, these results indicate that the tissue signature is an intrinsic feature of T cells from non-blood sites, and that our results from blood and BM are representative of T cells in these sites and representative of diverse individuals.

13The clustering analysis in Fig. 2 suggested that activated T cells were more similar across sites than resting counterparts. To uncover gene expression patterns that were conserved across T cell populations in different tissues throughout activation, we applied a new analytical method called single-cell Hierarchical Poisson Factorization (scHPF). The scHPF algorithm identifies a small number of expression patterns, called factors that vary coherently across cells. These factors can represent discrete, subpopulation-specific programs or continuous programs like T cell activation that are expressed as a gradient across cells in different stages of a biological process. We applied scHPF to merged resting and activated T cells from each tissue and donor separately and hierarchically clustered the resulting factors (Fig. 4a, Supplementary Figs. 7, 8a). This analysis revealed seven gene expression modules (three resting and four activated/functional) that were highly conserved across tissues and donors, for which the highest scoring genes formed interpretable gene signatures (Fig. 4a, Supplementary Fig. 8a, Supplementary Data 5). Modules were annotated based on known markers among their highest scoring genes, association with resting or activated states, and CD4:CD8 experssion ratio. The three modules associated with a resting state (Fig. 4a) included a Treg module defined by canonical genes (FOXP3, CTLA4, IRF4, TNFRSF4 (OX40)); a putative resting CD4+ Naive/Central memory (NV/CM) module enriched in CD4+ T cells and defined by genes associated with lymphoid homing, egress and quiescence (SELL, KLF2, LEF1, respectively); and a CD4+/CD8+ Resting module, distinguished by expression of IL7R, a receptor required for T cell survival,, and AQP3, which encodes a water channel protein of unclear function in lymphocytes. Importantly, the CD4+/CD8+ Resting module did not contain factors from blood and had the highest enrichment for the tissue signature identified in Fig. 3 (Supplementary Fig. 9).

14Four modules were associated with T cell activation and/or function, some of which were lineage-specific. A Proliferation module expressed by activated CD4+ and CD8+ lineages included genes associated with T cell activation/proliferation (IL2, LIF) and cell division (CENPV, G0S2, ORC6) (Fig. 4a). This module was also marked by expression of NME1, a metastasis suppressor/endonuclease-encoding gene not previously associated with T cells (Fig. 4a). An Interferon (IFN) Response module enriched among activated CD4+ T cells included multiple gene families associated with canonical IFN responses– (IFIT3, IFIT2, STAT1, MX1, IRF7, and JAK2). In contrast, CD8+ T cell-enriched modules included a Cytotoxic module, containing genes associated with cytotoxicity (GNLY, GZMK) and transcription factors associated with effector/memory differentiation (ZEB2, EOMES, ZNF683)–, and a Cytokine module with genes encoding chemokines and cytokines (CCL3, CCL4, CCL20, IFNG, IL10, TNF), inhibitory molecules (LAG3, CD226 (TIGIT), HAVCR2 (TIM3)), and the widely expressed homeobox protein HOPX. These results indicate a limited spectrum of functional states for human T cells across blood and tissue sites.

15To understand how these gene modules correspond to resting and activated states in CD4+ and CD8+ T cells, we visualized the average expression of their top-ranked genes on diffusion maps for each donor and tissue (Fig. 4b–e). This visualization defined activation trajectories with resting T cells on the left (blue) and activated T cells projecting to the right (red; Fig. 4b, c). In all four sites in both individuals, module expression for CD4+ T cell was positioned along activation trajectories from CD4 NV/CM Resting (left) to IFN-Response (middle) to Proliferation (right) (Fig. 4d). Expression of genes within the Proliferation module co-localized with peak expression of NME1 and IL2RA (Supplementary Fig. 8b, c), while the IFN Response module genes exhibited peak expression at the middle of the trajectory as exemplified by IFIT3 expression (top ranked gene) (Supplementary Fig. 8d), suggesting a potential intermediate activation state. In CD8+ T cells, the Cytokine module localized in the most activated cells for all sites also shown by IFNG expression (Fig. 4e, Supplementary Fig. 8e), while the Cytotoxic module was expressed among resting and activated cells (Fig. 4e). Therefore, scHPF takes an unbiased approach to uncover major functional states, reference signatures and activation trajectories for human T cells that are conserved across sites.

16The functional states identified for human CD8+ T cells in Fig. 4 were consistent in with those seen in vivo in mouse infection models. By contrast, the modules identified for CD4+ T cell activation revealed markers and functional states not typically associated with effector CD4+ T cells. We therefore assessed expression kinetics of the top-scoring genes in the Proliferation and IFN Response modules, NME1 and IFIT3, respectively, during the course of T cell activation ex vivo by qPCR. Expression of NME1 transcripts rapidly increased after TCR-stimulation, peaking between 16 and 24 h and remaining elevated for up to 72 h, for both CD4+ and CD8+ T cells compared to unstimulated controls, a pattern of expression similar to the canonical T cell activation marker IL2RA (Fig. 5a). Notably, the extent of activation-associated upregulation of NME1 transcripts was greater in CD4+ compared to CD8+ T cells, while IL2RA was more upregulated in CD8+ T cells (Fig. 5a). At the protein level, NME1 expression increased in CD4+ and CD8+ T cells after TCR-mediated stimulation from 24 to 120 h (Fig. 5b, upper), and with each successive round of T cell proliferation, while CD25 was expressed similarly, independent of cell division (Fig. 5b, lower). These results establish NME1 expression as a marker of T cell activation, coupled to the extent of proliferation.

17In contrast to NME1/IL2RA upregulation, expression of the interferon-inducible transcript IFIT3, showed transient upregulation by CD4+ T cells following TCR-stimulation, peaking at 16 h and returning to near baseline levels by 48 h post-stimulation (Fig. 5c). By contrast, induction of IFIT3 by culturing T cells with IFN-α (type I) or IFN-γ (type II) occurred rapidly (within 2 h) and persisted throughout the culture period (Fig. 5d, e). To identify the contribution of type I or type II IFN signaling to TCR-triggered IFIT3 induction, we included blocking antibodies to type I or Type II IFN in the cultures. While neutralizing antibodies for type I IFNs and IFNαR2 completely inhibited IFIT3 induction by Type I IFN, TCR-mediated upregulation of IFIT3 was unaffected (Fig. 5d). However, blockade of type II IFN signaling via a combination of anti-IFNγ and anti-IFNγR1 antibodies inhibited upregulation of IFIT3 by both exogenous IFN-γ and TCR-mediated stimulation (Fig. 5e). Importantly, blocking type II (or type I) IFN signaling did not inhibit T cell activation as assessed by induction of NME1 transcript expression, and addition of IFN-α or −γ did not induce NME1 expression (Fig. 5d, e). These results establish that the IFN-responsive state suggested by the scRNA-seq trajectories is recapitulated in real-time as part of an intermediate activation state driven by TCR-triggered IFN-γ production.

18We further assessed whether CD4+ T cells express IFIT3 and NME1 in vivo using a published scRNA-seq dataset of T cells isolated from the blood of dengue virus-infected patients, which contains a fraction of activated CD4+ T cells. Both NME1 and IFIT3 were expressed by CD4+ T cells from dengue-infected patients (Supplementary Fig. 10). These results show that genes associated with functional modules identified for CD4+ T cell activation are expressed in vivo.

19Although there have been several large-scale scRNA-seq studies of disease-associated T cells, these data are generally not placed in the context of T cell activation in healthy individuals. To demonstrate the utility of our resource as a reference point for human disease, we used UMAP to project recently reported scRNA-seq profiles of tumor-associated T cells from four different human cancers onto our map of T cell activation states. We merged all of our T cell data from four donors and four sites in a single UMAP embedding (Fig. 6a), colored by tissue site, donor, stimulation, cluster-level CD4/CD8 status, and CCL5 expression, indicative of effector status. We projected scRNA-seq profiles of tumor-associated T cells from four different human cancers,– (non-small cell lung cancer (NSCLC), colorectal cancer (CRC), breast cancer (BC), and melanoma (MEL)) onto this embedding to compare each tumor-associated T cell to healthy T cells (Fig. 6b, c). We also investigated expression of activation state and lineage markers in the healthy T cell embedding and tumor projections (Fig. 6c). Tumor-associated CD8+ T cells project onto healthy CD8+ T cells from all sites in both resting and activated states (Fig. 6b). Moreover, genes associated with TRM (CXCR6) and the Cytotoxic and Cytokine modules are all represented among tumor-associated CD8+ T cells (Figs. 6c, 7). By contrast, tumor-associated CD4+ T cells projected mostly onto resting blood and tissue T cells (Fig. 6b), while CD4+ T cell activation states and associated markers (NME1, IFIT3) were largely absent (Fig. 6c). Projecting tumor-associated T cells onto each individual tissue and blood donor yielded results consistent with projection onto the combined dataset (Supplementary Figs. 11–14). We note that projecting tumor-associated T cells onto our reference map using the alternate projection algoritham scmap showed similar results (Supplementary Fig. 15). This analysis reveals that tumor-associated T cells contain activated CD8+ T cell states, but lack the presence of functionally activated CD4+ T cell states.

20A hallmark of tumor-associated T cells is a state of hyporesponsiveness or functional exhaustion, marked by persistent expression of surface inhibitory markers including PD-1, CTLA4, LAG3, TIM3 and others, many of which are expressed following T cell activation,,. Some of these molecules (PD-1, CTLA4) are important targets for immunotherapy to promote anti-tumor immunity–. We compared expression of exhaustion and functional markers across healthy and tumor-associated T cells (Fig. 7; Supplementary Figs. 16, 17). Tumor-associated CD8+ T cells expressing exhaustion markers across all four tumor types project onto activated CD8+ T cells in our map, and express genes within the Cytokine module (CCL3, CCL4, XCL1, XCL2, and IFNG), and to a lesser extent Cytotoxic module (Fig. 7; Supplementary Figs. 16, 17). Interestingly, a subset of these tumor-associated CD8+ T cells, but not healthy T cells, express high levels of MKI67, associated with proliferating cells and other cell cycle control markers (Fig. 7, Supplementary Fig. 17). Therefore, tumor-associated T cells expressing exhaustion markers also express genes associated with normal CD8+ effector T cell function and ongoing proliferation.

21Human T cells persist in distinct anatomic sites, maintain protective immunity and surveillance, and are key targets for immune modulation in tumor immunotherapy, transplantation, and autoimmunity. Here, we used scRNA-seq profiling of resting and TCR-stimulated T cells from blood, lymphoid and mucosal tissues to generate a reference map of human T cells and understand how T cell homeostasis and function are related to the tissue site. Our findings demonstrate fundamental differences between T cells from tissues and blood, but similar functional and activation states across sites that are intrinsic to lineage; human CD4 T cell activation is defined by response to cytokines and proliferation while CD8+ T cells are defined by effector function. We further demonstrate that this high-resolution map of T cell homeostasis and activation across sites, lineages, and individuals can serve as a new baseline for defining human T cell states in disease.

22The study of healthy human T cells has largely focused on blood, while the majority of T cells persist in diverse lymphoid, mucosal and barrier sites,. Human tissue T cells are largely memory subsets, comprising tissue-resident (TRM) and non-resident (TEM, TCM) populations; TRM predominate in mucosal sites, while TEM are found in spleen, LN and BM,,. The transcriptional differences and functional relationship of these tissue-localized TEM to blood TEM has been unclear. Importantly, profiling using scRNAseq enabled unambiguous assessment of T cell-intrinsic differences in tissue versus blood T cells. We show here that TEM from all tissue sites examined (LG, LN, BM) exhibit fundamental changes in expression of cytoskeletal, cell-matrix interaction, and proliferative genes compared to blood TEM cells, indicating alterations in cellular structure. These tissue-intrinsic expression patterns are in addition to TRM-associated genes identified in previous studies, or functional adaptations of T cells to specific tissue sites,. Whether T cells require these changes in gene expression to enter or reside within the tissue architecture, and if their loss of expression enables tissue T cell egress to circulation remains to be established.

23Our results reveal conserved functional states for human blood and tissue-derived T cells. CD8+ T cells segregate into two major effector subsets based on expression of genes involved in cellular cytotoxicity (Cytotoxic module) and myriad cytokines and chemokines (Cytokine module). These predominant effector states within activated human CD8+ T cells are consistent with results showing that mouse CD8+ T cell activation triggers an effector differentiation program,. We identified two major activation states that were not associated with effector function: one associated with proliferation and IL-2 production, and a second state enriched in CD4+ T cells which is characterized by induction of multiple IFN-responsive genes and gene families including IFIT3, MX1, IRF7, and others. Induction of this IFN-response state is due to TCR-mediated IFN-γ production (likely autocrine responses), and appears as a kinetic intermediate early after CD4+ T cell activation, and prior to induction of the proliferative program. Identification of a functional state for T cells based on cytokine responses is distinct from T cell functional states that are typically defined based on cytokine secretion profile. We propose that the IFN-responsive state for human CD4+ T cells may serve an autoregulatory function to temper high IFN levels produced by predominant memory responses, and ongoing responses to persistent viruses.

24This scRNA-seq analysis provides a high-resolution map for human T cells from which to define T cell states in disease. We demonstrate this approach by projecting T cell profiles from human tumors onto our reference map. We identify predominant CD8+ T cell effector populations, Tregs, and resting (but not activated) CD4+ T cells in datasets derived from diverse tumor types (breast, lung, skin, colon). Interestingly, the tumor-associated CD8+ T cells exhibited transcriptional features similar to healthy activated CD8+ T cells including expression of multiple effector molecules such as perforin, IFN-γ and chemokines. We also examined the expression of multiple markers associated with exhaustion, a functionally hyporesponsive state found in tumor-infiltrating T cells targeted by checkpoint blockade immunotherapies,,. Interestingly, exhaustion markers were upregulated along with CD8-associated cytokines in activated T cells from both healthy and tumor tissues, emphasizing the importance of obtaining baseline healthy profiles for high resolution analysis of T cells on the single cell level. Moreover, subsets of these CD8+ T cells in all four tumors expressed higher levels of proliferation markers compared to healthy T cells, consistent with a recent report that T cells expressing exhaustion markers in melanoma exhibit aberrant proliferation. This analysis can therefore enable precise identification of features of resting and activated T cells that are associated with tissues, activation and disease.

25Our high-resolution analysis of human T cells across sites, lineages, and activation states provides insights into human T cell adaptations to tissues and their intrinsic activation properties. Limitations of the study include that the select tissues and donors profiled here may not include the full diversity of T cell transcriptional programs throughout the body, and that quantification of cell types may be subject to dissociation biases between the individual tissues,. Importantly, our dataset establishes a starting point for the integration of other T cell scRNA-seq datasets to ultimately capture the full breadth of T cells states in humans. International collaborative efforts like the Human Cell Atlas are now underway, generating comprehensive scRNA-seq datasets profiling a diverse range of cells, including T cells and their transcriptional states. Recently developed computational tools including scVI, mutual nearest neighbors, Seurat v3, Conos, and Scanorama will be useful for this integration and as a guide for future studies. In this way, our novel reference map can serve as a valuable resource for the ongoing study of human T cell immunity in disease, immunotherapies, vaccines and infections, with the ultimate goal of diagnosing, screening and monitoring immune responses.

26We obtained human tissues from deceased, brain-dead donors at the time of organ acquisition for clinical transplantation through an approved research protocol and MTA with LiveOnNY, the organ procurement organization for the New York metropolitan area. Obtaining tissue samples from deceased organ donors does not qualify as “human subjects” research, as confirmed by the Columbia University Institutional Review Board (IRB). Donors were free of chronic disease, cancer and chronic infections such as Hepatitis B, C, and HIV. Clinical and demographic data regarding organ donors used in this study are summarized in Supplementary Table 1. We obtained peripheral blood from healthy consenting adult volunteers by venipuncture, through a protocol approved by the Columbia University IRB and have complied with all relevant ethical regulations for work with human participants.

27Tissues acquired from donors were maintained in cold saline during transport to the laboratory, typically within 2–4 h of procurement. We isolated mononuclear cells from donor lungs, lung-draining lymph nodes (LN) and bone marrow (BM) as previously described,. Briefly, lungs were flushed with cold complete medium (RPMI 1640, 10% FBS, 100 U/ml penicillin, 100 μg/ml streptomycin, 2 mM L-glutamine) and left lateral basal segment of the lung was isolated. LN were isolated from the hilum, near the intersections of major bronchi and pulmonary veins and arteries, removing all fat. To obtain mononuclear cell suspensions, LN and lung tissues were mechanically processed using a gentleMACS tissue dissociator (Miltenyi Biotec), enzymatically digested (complete medium with 1 mg/ml collagenase D, 1 mg/ml trypsin inhibitor and 0.1 mg/ml DNase for 1 h at 37 °C in a mechanical shaker) and centrifuged on a density gradient using 30% Percoll Plus (GE Healthcare). BM was aspirated from the superior iliac crest. For BM and peripheral blood, we isolated mononuclear cells by density gradient centrifugation using Lymphocyte Separation Medium (Corning). T cells were enriched from all samples using magnetic negative selection for CD3+ T cells (MojoSort Human CD3+ T cell Isolation Kit; BioLegend), followed by a dead cell removal kit (Miltenyi Biotec), resulting in 80–90% purity. We cultured 0.5–1 × 106 CD3+ enriched cells from each donor tissue for 16 h at 37 °C in complete medium, with or without TCR stimulation using Human CD3/CD28 T Cell Activator (STEMCELL Technologies). After stimulation, dead cells were removed as above before cell isolation for single-cell RNA-seq.

28Single-cell suspensions were loaded onto a Chromium Single Cell Chip (10x Genomics) according to the manufacturer’s instructions for co-encapsulation with barcoded Gel Beads at a target capture rate of ~5000 individual cells per sample. We barcoded captured mRNA was barcoded during cDNA synthesis and converted the barcoded cDNA into pooled single-cell RNA-seq libraries for Illumina sequencing using the Chromium Single Cell 3′ Solution (10x Genomics) according to the manufacturer’s instructions. All samples for a given donor were processed simultaneously with the Chromium Controller (10x Genomics) and the resulting libraries were prepared in parallel in a single batch. We pooled all of the libraries for a given donor, each of which was barcoded with a unique Illumina sample index, for sequencing in a single Illumina flow cell. All of the libraries were sequenced with an 8-base index read, a 26-base read 1 containing cell-identifying barcodes and unique molecular identifiers (UMIs), and a 98-base read 2 containing transcript sequences on an Illumina HiSeq 4000. Cell counts and transcript detection rates are summarized in Supplementary Table 2.

29Prior to gene expression analysis, raw sequencing data was corrected for index swapping, a phenomenon that occurs during solid-phase clonal amplification on the Illumina HiSeq 4000 platform and results in cross-talk between sample index sequences. We corrected index swapping using the algorithm proposed by Griffiths et al. First, we aligned the reads associated with each sample index to GRCh38 (GENCODE v.24) using STAR v.2.5.0 after trimming read 2 to remove 3′ poly(A) tails (>7 A’s) and discarding fragments with fewer than 24 remaining nucleotides as described in Yuan et al.. For each read with a unique, strand-specific alignment to exonic sequence, we constructed an address comprised of the cell-identifying barcode, unique molecular identifier (UMI) barcode, and gene identifier. Next, we counted the number of reads associated with each address in each sample. Because of index swapping, we found that some addresses occurred in multiple samples at much higher frequencies than one would expect by chance. For the vast majority of addresses, there was a single sample containing most of the associated reads. If >80% of reads for a given address were associated with a single sample (e.g., a single index sequence), we kept all of the reads corresponding to that address in that sample and removed all of the reads associated with that address from all other samples. We also identified addresses for which no sample contained >80% of the corresponding reads and removed all of these reads from all samples. After correcting for index swapping, we collapsed amplification duplicates using the UMIs and corrected errors in both the cell-identifying and UMI barcodes to generate a preliminary matrix of molecular counts for each cell as described previously.

30We filtered the cell-identifying barcodes to avoid dead cells and other artifacts as described in Yuan et al.. Briefly, we removed all cell-identifying barcodes where >10% of molecules aligned to genes expressed from the mitochondrial genome or for which the ratio of molecules aligning to whole gene bodies (including introns) to molecules aligning exclusively to exons was >1.5. Finally, we also removed cell-identifying barcodes for which the average number of reads per molecule or average number of molecules per gene deviated by >2.5 standard deviations from the mean for a given sample.

31Thoroughly removing non-T cells from the data set is complicated by technical issues such as molecular cross-talk, multiplet capture, and a broad coverage distribution. We developed a procedure to remove non-T cells that accounts for these issues by identifying both individual cells and clusters of cells that are enriched in expression of a blacklisted gene set that is highly specific to contaminating cell types. We began by clustering the single-cell profiles within each sample using a pipeline that we reported previously,. Briefly, we identified highly variable genes that are likely markers of specific subpopulations by normalizing the molecular counts for each cell to sum to one, ordering all genes by their normalized expression values, and computing a drop-out score dsg for each gene g defined as:where fg is the fraction of cells in which we detected g and fgmax is the maximum fg in a 25-gene rolling window centered on g. We selected genes with dsg > 0.15 or with dsg > 6σds + < dsg >, where σds and < dsg > are the standard deviation and mean of the dropout score distribution. Using these genes, we computed a cell-by-cell Spearman’s correlation, from which we constructed a k-nearest neighbor’s graph (k = 20) and used this as input for the Phenograph implementation of Louvain clustering to identify cellular subpopulations.

32Next, we used the pooled normalization approach described by Lun et al. as implemented in the scran package with the computeSumFactors function to compute size factors for each cell,. We supplied the computeSumFactors function with the cluster identifiers obtained from Phenograph to account for cell type-specific coverage differences. Using the resulting normalized expression profiles, we identified Phenograph clusters with positive enrichment of average CD3D and TRAC expression and labeled these clusters as T cell clusters (Supplementary Fig. 1). Within each sample, we conducted differential expression analysis between all pairs of T cell and non-T cell clusters via the Wilcoxon rank-sum test using the SciPy function ranksums and Benjamini-Hochberg corrected p-values with the StatsModels function multipletests in Python, yielding p values padj. Finally, we established an initial blacklist of genes that are highly specific to the non-T cell clusters by taking any gene with padj < 0.001 and greater than 10 fold-enrichment in a non-T cell cluster for any of the above pairwise comparisons in any sample. To refine the blacklist and avoid including genes that are specific to T cell subsets found in only a limited set of samples or clusters, we also generated a whitelist of genes with positive enrichment in any T cell cluster. We removed any member of this whitelist from the initial blacklist to produce a final, refined blacklist containing 744 genes highly specific to contaminating cell types (Supplementary Data 1). As expected, genes on the final blacklist included markers of epithelial cells, dendritic cells, mast cells, B cells, neutrophils, and red blood cells.

33The blacklist of genes was used to remove cells from the T cell clusters that are either improperly clustered (unlikely to be T cells) or potentially multiplets (a cell-identifying barcode co-encapsulated both T cells and non-T cells). Importantly, because of molecular cross-talk in scRNA-seq libraries from PCR recombination, we only considered a cell to be expressing a blacklisted gene if the average number of reads supporting the detected molecules was above a certain threshold. This threshold depends on the average depth to which we sequenced the libraries in a given sample. The distributions of the number of reads-per-molecule are generally bimodal for a given sample. We assume that the mode with lower read counts per molecule arises from PCR recombination in which a molecule originating from one cell receives the cell-identifying barcode of a different cell at an intermediate point in PCR, thereby resulting in a detected molecule supported by an unusually small number of reads (i.e., amplicons). We therefore considered the sample h with the highest coverage (and therefore the clearest separation between the two modes) and took the minimum point between the two modes in the reads-per-molecule distribution to be the threshold number of reads per molecule, Th, below which a detected molecule would be considered to arise from cross-talk. We extrapolated a reads-per-molecule threshold for each of the other samples s as:where RPMs is the average number of reads per molecule detected in sample s.

34Finally, for each cell c in a sample with threshold Ts, we computed bc, the per-cell fraction of blacklisted genes detected with an average number of reads per molecule above Ts. As expected, bc was typically bimodally distributed within each sample (Supplementary Fig. 1e). The vast majority of cells in the lower mode were in the T cell clusters described above, while the high mode was composed mainly, but not exclusively, of cells from non-T cell clusters (Supplementary Fig. 1e). In each sample, we fit a Gaussian to bc’s distribution across cells assigned to T cell clusters and established a threshold at two standard deviations above the fitted mean. We considered any cell with bc above this threshold and any cell that clustered among the non-T cell clusters to be a non-T cell and discarded these cells from all downstream analysis.

35Once we had identified the T cells from each sample using the methodology described above, we merged resting and activated samples from all of the tissues in each donor and clustered the T cells from the two donors separately to generate Fig. 1b, c. We used the methodology described above to identify a set of highly variable genes for each sample (including the blood samples), and then merged those sets to generate a large list of 315 highly variable genes (Supplementary Data 2) with which we clustered the merged samples from both donors. We computed Louvain clusters from the two merged data sets with k = 12 and a minimum cluster size of 100 cells using a k-nearest neighbors graph constructed from the Spearman’s correlation matrix calculated using the 315 highly variable genes. We used the Python implementation Uniform Manifold Approximation and Projection (UMAP) to produce the two-dimensional projections shown in Fig. 1b, c. To obtain CD4/CD8 ratios for each cluster, we first computed the expression level of CD4 and CD8A in each cell using the normalized counts from computeSumFactors as described above. For both CD4 and CD8A, we then computed the average log2(normalized counts + 1) for each cluster and normalized this value by the average log2(normalized counts + 1) for all cells. We then took the log-ratio of these values for CD4 and CD8A to generate Fig. 1b, where all the cells in each cluster are labeled with the cluster’s log-ratio. Differentially expressed genes for cells in each cluster versus all other cells were determined using a binomial test (Supplementary Data 3, 4).

36To project the data obtained from blood T cells onto the tissue-derived profiles from each organ donor, we first merged the scRNA-seq profiles from both blood donors. We note that the scRNA-seq data from blood were subjected to the same computational procedure described above for eliminating non-T cell profiles. We used the same highly variable gene set (Supplementary Data 2) that was used in the original UMAP model of each organ donor to compute a Spearman’s correlation matrix between the blood and tissue profiles. We then projected the blood T cell profiles onto the UMAP embeddings for each of the two organ donors using the transform function in UMAP. We note that the organ donor UMAP embeddings used for this analysis are slightly different from what appears in Fig. 1b, c, because a small number of genes in the highly variable gene set were eliminated due to lack of expression in the blood. We also note that a small modification to the UMAP source code was needed to accommodate the use of Spearman’s correlation as a similarity metric (available at https://github.com/simslab/umap_projection).

37We confirmed our findings using scmap, a previously published scRNA-seq data projection algorithm. When projecting blood T cells onto the tissue T cells from Tissue Donors 1 and 2, scmap yielded projections that were consistent with UMAP (Supplementary Fig. 4a, b, e, f) with projection coordinates that were highly correlated across both data sets (Supplementary Fig. 4c, d, g, h).

38To generate the cell number heatmaps in Fig. 2 and Supplementary Fig. 3, we first computed a centroid position in the UMAP embedding for each condition, subset and tissue combination in the tissue data based on the Louvain clustering described above for Fig. 1b, c. For example, for Donor 2, we computed the average position of LG-, BM-, and LN-derived resting. We then identified the nearest condition, subset and tissue combination for each cell in the blood samples based on the Euclidean distance between a given blood-derived cell’s position in the UMAP model (following projection of the blood data onto the tissue UMAP model) and each centroid position. The heatmaps summarize the results of these calculations, providing the number of blood-derived cells that are closest to each condition, subset and tissue combination in the organ donor data.

39To identify a tissue-specific T cell signature, we compared the expression profiles of effector memory cells from resting LG, BM, and LN T cells from the two tissue donors to resting blood T cells from the two blood donors. We found CCL5 to be an extremely highly expressed marker of effector-memory (TEM) cells that exhibited strong anti-correlation with SELL, a marker of non-effector memory cells, in all of our resting samples (Supplementary Fig. 5a). We also found that the average number of reads per molecule for CCL5 was bimodally distributed, consistent with spurious detection of CCL5 in a population of cells due to PCR recombination (Supplementary Fig. 5b). For each sample, we used the point between these two modes where the probability density was minimal as a threshold for the minimum average number of reads per molecule of CCL5 required for a cell to be considered positive for CCL5. For each sample, we normalized the matrix of molecular counts for the CCL5+ TEM cells using the computeSumFactors function in scran to compute size factors for each cell,. For each tissue site, we then identified differentially expressed genes for all four pairwise comparisons of resting tissue to resting blood CCL5+ T cells (tissue donor 1 vs. blood donor A, tissue donor 2 vs. blood donor A, etc.) using the Wilcoxon rank-sum test with the SciPy function ranksums and computed Benjamini-Hochberg corrected p-values with the StatsModels function multipletests in Python after removing genes from the blacklist described above (Supplementary Data 1). For each tissue, we took all genes with padj < 0.05 and fold-change > 2 in all 4 pairwise comparisons to comprise a tissue-specific effector memory T cell signature (Fig. 3).

40Next, all of the genes in the tissue-specific effector memory T cell signature and computed the average normalized expression of the resulting gene set to obtain Fig. 3d. Z-scored normalized expression for each of these genes appears in the heatmap in Fig. 3e for each site/donor, which also includes a set of blood T cells with outlier expression of the tissue-enriched gene signature (blood T cells with average expression within one standard deviation of that of the tissue T cells as indicated by the dashed line in Fig. 3d). For Fig. 3f, we took all genes with padj < 0.05 (Wilcoxon and correction described above) and two-fold higher expression in either the tissue- or blood-associated T cells in all four pairwise comparisons. We then constructed a Spearman’s correlation coefficient between pairs of resting CCL5+ cells in the dataset across these genes and used this to generate a UMAP embedding.

41We used Gene Set Enrichment Analysis (GSEA) to assess the enrichment of the tissue-associated T cell signature in T cells profiled from the bone marrow of living donors and its depletion in T cells profiled from two additional blood samples in Supplementary Fig. 6. We obtained 10x Genomics Chromium scRNA-seq profiles of bone marrow from 20 individuals from GEO accession GSE120221. We merged the data from all of the samples and clustered the merged scRNA-seq profiles using Phenograph as described above. To computationally isolate TEM cells, we took all CCL5+ cells that occurred in any Phenograph cluster that was positively enriched in TRAC expression. We then compared the expression profiles of CCL5+ T cells from the bone marrow of living donors to the resting, CCL5+ T cells from the two blood donors from this study (Donors A and B) using the Wilcoxon rank-sum test as described above. After ranking all genes for which a test statistic could be computed by fold-change (comparing bone marrow to blood), we used GSEA (pre-ranked, “classic” mode with 10,000 permutations) to calculate the enrichment of the tissue-associated T cells signature among the differentially expressed genes. We used the Java implementation of GSEA that is freely available from http://software.broadinstitute.org/gsea/index.jsp.

42To assess the depletion of the tissue-associated T cell signature in other blood data sets relative to the tissue samples collected here, we obtained 10x Genomic Chromium scRNA-seq data from PBMCs of a healthy donor and purified T cells of a healthy donor from 10x Genomics (Donor C PBMCs: pbmc8k data set from https://support.10xgenomics.com/single-cell-gene-expression/datasets/2.1.0/pbmc8k; Donor D purified T cells: t_4k data set from https://support.10xgenomics.com/single-cell-gene-expression/datasets/2.1.0/t_4k). We performed clustering separately on the two samples using Phenograph and computationally isolated TEM cells as described above for the live donor bone marrow samples. We then performed differential expression analysis using the Wilcoxon rank-sum test to compare TEMs from Donors C and D to those in each tissue from the resting samples collected for this study. For each tissue site, we merged the resting TEMs from Donors 1 and 2. Finally, we performed GSEA as described above to assess the depletion of the tissue-associated T cell signature in Donors C and D for each tissue site comparison.

43We applied Single-cell Hierarchical Poisson Factorization (scHPF), a method that we recently reported for de novo discovery of gene expression signatures in scRNA-seq data, to the merged activated and resting cells for each tissue and donor. Given a molecular count matrix, scHPF identifies a small number of latent factors that explain both continuous and discrete expression patterns across cells. Each gene has a score for each factor, quantifying the gene’s contribution to the associated expression pattern. Likewise, each cell assigns a score to each factor, which reflects the contribution of the factor to the observed expression in the cell.

44We applied scHPF to each tissue and blood sample after merging their respective resting and activated datasets. We considered only genes with GENCODE protein coding, T cell receptor constant or immunoglobulin constant biotypes, excluded genes on the previously described blacklist, and removed genes detected in fewer than 0.1% of cells in a given merged dataset. scHPF (version 0.1) was run with default parameters for seven values of K, the number of factors, equal to all values between 6–12, inclusively. This resulted in seven candidate scHPF factorizations per merged dataset. We then selected K to avoid factors with significant overlap in their gene signatures. For each dataset and value of K, we calculated NK: the maximum pairwise overlap of the 300 highest-scoring genes in each factor for the corresponding scHPF model. We considered overlap significant if p < 0.05 by a hypergeometric test with a population size equal to the number of unfiltered genes in the tissue sample and NK observed successes. Finally, for each dataset, we selected the model with maximum K such that p >= 0.05 (Supplementary Fig. 7). This procedure resulted in eight factorizations: six from tissue donors (lung, BM, and LN from each of two organ donors) and two factorizations from the blood of living donors. We defined each factors’ CD4/CD8 bias as the log2 ratio of its mean cell score in CD4+ and CD8+ T cells.

45To discover common patterns of expression across tissues and donors, we performed unsupervised clustering of all factors for tissue-derived cells. First, we calculated Pearson correlation on the union of the fifty highest and lowest scoring genes in each factor for each tissue factorization (2291 genes total) using the Python pandas package’s DataFrame.corr function. Next, we hierarchically clustered the factor-factor correlation matrix using scipy.cluster.hierarchy.linkage with method = ‘average’ and scipy.cluster.hierarchy.dendrogram (Supplementary Fig. 8a). This defined clusters of tightly correlated expression patterns, which we call expression modules. We focused on seven modules (out of nine) whose factors had mean pairwise correlations greater than 0.25. Most modules contained at least one factor from each tissue and donor. To identify the top genes in each module (Fig. 4a, Supplementary Data 5), we ranked genes by their mean gene score across all constituent factors. The CD4 IFN response module contained two factors from Donor 2 BM; however, one of the two factors was far more tightly correlated with the rest of the factors in the module than the other. As the top genes in the module were nearly identical with and without the less tightly-correlated factor, we excluded it from the module in downstream analyses for clarity.

46We used the factorizations described above to compute T cell activation trajectories by diffusion component analysis. We first converted the cell score matrix obtained from the factorization of each resting/activated merged tissue or blood sample into a cell-by-cell Euclidean distance matrix. We then extracted the distance submatrices corresponding to the CD4 and CD8 clusters in each sample as defined from the merged analysis of all samples from each donor described above. We used the two resulting distance submatrices to compute diffusion components for CD4 and CD8 activation with the C++ Accelerated Python Diffusion Maps Library (DMAPS) with a kernel bandwidth of four. The diffusion maps shown in Fig. 4b–e each show the first two diffusion components which we define as the two diffusion eigenvectors with the second- and third-highest eigenvalues scaled by the diffusion eigenvector with the largest eigenvalue.

47To evaluate the expression of T cell surface markers by flow cytometry, we incubated tissue and blood cell suspensions with Human TruStain FcX (BioLegend) and stained with following fluorochrome-conjugated antibodies: CD3 (UCHT1, BD Biosciences; OKT3, BioLegend), CD4 (SK3, BD Biosciences; SK3, Tonbo Biosciences), CD8 (SK1, BioLegend; RPA-T8, BD Biosciences), CCR7 (G043H7; BioLegend), CD45RA (HI100; BioLegend), CD25 (BC96; BioLegend), CD127 (A019D5; BioLegend), CD69 (FN50; BioLegend), CD103 (Ber-ACT8; BioLegend), CD45 (HI30; BioLegend), and Fixable Viability Dye eFluor 780 (eBioscience). For stimulation/proliferation assays, we magnetically enriched for CD3+ T cells from single cell suspensions, stained cells with Cell Proliferation Dye eFluor 450 (eBioscience), and cultured cells for up to 120 h with or without TCR stimulation as above. At indicated time points, we performed intercellular staining of NME1 (11615-H07E; Sino Biological) using a Foxp3/Transcription Factor Staining Buffer Kit (Tonbo Biosciences) for fixation and permeabilization of cells according to manufacturer’s instructions. We acquired cell fluorescence data using a BD LSR II flow cytometer and used FCS Express (De Novo Software) for analysis. The results are summarized in Supplementary Fig. 2 and the gating strategy is shown in Supplementary Fig. 18a.

48PBMC were magnetically enriched for CD3+ T cells, and sorted for live CD4+ and CD8+ T cells (singlets, FSClowSSClow, and Viability Dye-) using a BD Influx cell sorter (Supplementary Fig. 18b). Sorted cells were cultured in complete medium with or without anti-CD3/anti-CD28 stimulation as above for 2–72 h. For dissecting the contribution of type I and type II IFN signaling to gene expression, cells were pre-incubated with Human Type 1 IFN Neutralizing Antibody Mixture (PBL Assay Science, Cat# 39000-1) according to manufacturer’s instructions, or 1 µg/mL of both anti-IFNγ (R&D Systems, MAB285, clone # 25718) and anti-IFNγR1 (R&D Systems, MAB6731, clone # 92101). As a control, CD4+ T cells were activated with 1000 units/mL of recombinant human IFNα2 (PBL Assay Science, Cat#11101-1) or 10 ng/mL recombinant human IFNγ (Peprotech, Cat# 300-02). Control and stimulated CD4+ and CD8+ T cells were harvested at indicated time points and RNA isolated using a RNeasy Micro Kit (Qiagen) with on-column DNase digestion. We converted RNA to cDNA via SuperScript IV VILO Master Mix (Invitrogen) and performed quantitative real-time PCR (qPCR) on a Viia 7 Real-Time PCR system (Applied Biosystems) using TaqMan Gene Expression Assays (NME1 Hs00264824_m1; IL2RA Hs00907777_m1; IFIT3 Hs00155468_m1; TBP Hs00427620_m1) and TaqMan Fast Advanced Master Mix, all from ThermoFisher Scientific. Quntitative PCR (qPCR) reactions were set up according to manufacturer’s instructions and fold changes between stimulated and unstimulated cells at each time point were calculated using the ∆∆ cycle threshold method in ExpressionSuite Software (ThermoFisher Scientific) with TBP as a reference gene.

49We projected scRNA-seq profiles of tumor-associated T cells from four different tumor types onto a UMAP embedding of resting and activated T cells from our combined tissue and blood data set using the methods described above for projecting blood T cells onto embeddings of the tissues. Briefly, we used the highly variable gene set from Supplementary Data 2 to generate a UMAP embedding of our tissue/blood data from a Spearman’s correlation matrix. We did not find any qualitative differences between UMAP embeddings when donor-specific genes were removed (Supplementary Fig. 19). We then projected the tumor-associated T cell profiles onto this embedding using the transform function in UMAP. Tumor-associated T cells from non-small cell lung cancer (NSCLC) and breast cancer (BC), which were profiled using the 10x Genomics Chromium platform, were obtained from https://gbiomed.kuleuven.be/scRNAseq-NSCLC and GEO accession GSE114724 (samples BC09, BC10, and BC11), respectively. For these two data sets, we used the UMI-corrected molecular counts provided by the authors. T cells from colorectal cancer (CRC) and melanoma (MEL), which were profiled using SMART-seq, were obtained from GEO accessions GSE108989 and GSE120575 (pre-treated samples only). For these two data sets, we used the TPM values provided by the authors. We note that the tissue/blood embedding was re-computed for each projection and is therefore slightly different in each case because not all of the processed data sets from the tumor studies contained all of the genes in Supplementary Data 2.

50The resulting projections are displayed in Fig. 6 in three different ways. In the top row, the projections are displayed as contour plots of estimated probability density (kernel density estimates) with a maximum of 14 contours. In the second row, we used a hexbin two-dimensional histogram of the number of cells in each bin with the colorbars normalized such that the intensity can be compared across samples (e.g., scaled so that the melanoma projection can be compared to the CRC projection). Finally, we also show where individual tumor-associated T cells project in subsequent rows along with gene expression values for several key markers. In Fig. 7, we show the average expression of several canonical exhaustion markers in individual cells. The markers used for this analysis were PDCD1, CTLA4, LAG3, LAYN, TIM-3, CD244, and CD160. We applied the same methodology to project the tumor-associated T cell profiles onto independent UMAP embeddings for each donor as shown in Supplementary Figs. 11–14.

51As above in Fig. 2, we validated the tumor-associated T cell projection analysis using scmap. Mapping the tumor-associated T cells onto our reference tissue and blood T cell dataset using scmap generated projections that were consistent with UMAP, with coordinates that were highly correlated (Supplementary Fig. 15).

52To analyze the expression of IFIT3 and NME1 in the context of virus infection, we analyzed scRNA-seq profiles of peripheral blood from dengue virus-infected patients (GSE116672). We clustered the data using the methodology described above and isolated T cell clusters based on enrichment of TRAC expression. We then generated the UMAP embedding shown in Supplementary Fig. 10 using the methodology described above and the same 315-gene set used throughout this study (Supplementary Data 2).

53Further information on research design is available in the Nature Research Reporting Summary linked to this article.