Transcriptional heterogeneity of tumor-associated high endothelial venules defines inflammatory and stress-metabolic states with distinct prognostic associations

Article type: Research Article

Keywords: High endothelial venules, Single-cell RNA sequencing, Transcriptomics, Prognosis

Authors: Hwal-Seok Choi, Yongjun Lee, Jin-Woo Choi, So-Jeong Kim, Hyunji Kim, Hyo-Kyung Pak, Jin Roh, Chan-Sik Park

Affiliations: https://ror.org/02c2f8975grid.267370.70000 0004 0533 4667Department of Pathology, Asan Medical Center, University of Ulsan College of Medicine, Seoul, Republic of Korea; https://ror.org/02c2f8975grid.267370.70000 0004 0533 4667Department of Biochemistry and Molecular Biology, Asan Medical Center, University of Ulsan College of Medicine, Brain Korea 21 Project, Seoul, Republic of Korea; https://ror.org/02c2f8975grid.267370.70000 0004 0533 4667Asan Institute for Life Science, Asan Medical Center, University of Ulsan College of Medicine, Seoul, Republic of Korea; https://ror.org/03tzb2h73grid.251916.80000 0004 0532 3933Department of Pathology, Ajou University School of Medicine, Suwon, Republic of Korea

License: © The Author(s) 2026 CC BY 4.0 Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

Article links: DOI: 10.1007/s12672-026-05162-2  |  PubMed: 42104153  |  PMC: PMC13323434

Relevance: Moderate: mentioned 3+ times in text

Full text: PDF (2.4 MB)

Introduction

The tumor microenvironment plays a pivotal role in driving cancer progression, shaping therapeutic responses, and determining patient prognosis [ref. 1]. The tumor vasculature not only provides oxygen and nutrients but also serves as a critical gateway for regulating the entry of immune cells into tumor tissue [ref. 2]. In particular, high endothelial venules (HEVs) are specialized post-capillary venules that mediate lymphocyte trafficking into secondary lymphoid organs, and they have been observed across a variety of human cancers [ref. 3]. As evidence accumulated that HEVs within tumors may modulate the microenvironment through lymphocyte recruitment, studies of tumor-associated HEVs (TA-HEVs) have surged [ref. 4–ref. 7].

Early histological studies revealed that TA-HEVs are associated with favorable clinical outcomes across several cancer types, including breast cancer, gastric cancer, and head and neck cancers [ref. 2, ref. 8, ref. 9]. However, others have reported neutral or even adverse associations depending on tumor type, anatomic localization, and microenvironmental context [ref. 8, ref. 10]. These discordant associations raise the possibility that TA-HEVs exist in distinct functional states shaped by local cues. Most prior studies have evaluated TA-HEVs within a single tumor type, with integrative efforts mostly limited to cross-study or meta-analytic approaches [ref. 10]. In addition, a systematic, state-resolved characterization of TA-HEVs across cancers has been lacking.

In the present study, we integrated publicly available single-cell RNA sequencing (scRNA-seq) datasets from 11 cancer types to construct a comprehensive atlas of tumor-associated endothelial cells (TA-ECs) and to identify TA-HEV populations, thereby enabling a state-resolved view of their transcriptional and functional heterogeneity.

Results

Transcriptional classification of tumor-associated endothelial cells

To reliably identify TA-HEVs across heterogeneous scRNA-seq datasets (Fig. 1A), we profiled 47,748 TA-ECs from 11 solid tumor types using canonical markers (PECAM1, VWF, CDH5, PLVAP). Unsupervised clustering revealed six EC subtypes (Fig. 1B)—venous, arterial, capillary, tip-like, lymphatic, and proliferative—defined by key differentially expressed genes (Supplementary Table S1): ACKR1 and SELP (venous), SEMA3G and GJA4 (arterial), CD36 and CA4 (capillary), ESM1 and TP53I11 (tip-like), CCL21 and PDPN (lymphatic), and MKI67 and ASPM (proliferative) (Fig. 1C-E). Venous ECs were the most abundant (50.7%), while proliferative ECs were rare (2.5%) (Fig. S1). At the tumor-type level, venous ECs were the dominant population in most cancers, whereas lymphatic and proliferative ECs were enriched in specific tumor contexts (Fig. 1F). These results indicated that TA-ECs can be clearly classified into six transcriptionally coherent subtypes, providing a transcriptional framework for the identification of the TA-HEV populations.

Identification of TA-HEV subclusters within venous ECs

Next, we focused on venous ECs from the initial analysis to identify the TA-HEV EC population. Secondary unsupervised clustering identified 20 transcriptionally coherent subclusters (Fig. 2A–B and Supplementary Table S2). The top DEG heatmap (Fig. 2B) demonstrates that HEV-associated marker genes, including LRG1 and MADCAM1, are coherently co-enriched in the TA-HEV subclusters, with substantially reduced enrichment in the remaining venous EC subclusters. UMAP visualization showed that these subclusters were evenly distributed across tumor types without apparent dominance by a specific tumor type (Fig. S2). Based on HEV marker expression and HEV signature gene scores, five subclusters (0, 1, 5, 6, 12) were defined as classified as TA-HEVs, four (7, 9, 13, 14) as immature venous ECs, and the rest as mature venous ECs (Fig. 2C-D and Fig. S3). Pathway analysis based on differentially expressed genes (DEGs) of each venous subgroup showed that TA-HEVs were enriched in immune-related pathways, including leukocyte cell-cell adhesion (Fig. 2E and Supplementary Table S3). Immature ECs were enriched in developmental and extracellular matrix organization, and mature ECs were enriched in cytoskeletal regulation. These were consistent with their annotated biological identities. Pseudotime analysis revealed a trajectory from immature to mature ECs, with TA-HEVs occupying a coherent region along the venous endothelial trajectory (Fig. 2F), consistent with their annotation as a specialized venous endothelial state.

Functional heterogeneity of TA-HEVs reveals inflammatory and stress-metabolic subtypes

Pathway analysis based on DEGs of the five TA-HEV subclusters defined above revealed two functional groups: inflammatory and stress-metabolic (non-inflammatory) (Fig. 3A and Supplementary Table S4). Among inflammatory TA-HEVs, subcluster 0 was enriched in pathways related to innate immune stimulation and leukocyte recruitment and expressed high IL6 (Fig. 3B and Fig. S4), suggesting early immune activation. Subcluster 12 was enriched in leukocyte migration and cytokine signaling, consistent with a role in immune activation and functional modulation (Fig. 3C). Subcluster 1 was enriched in MHC-II–mediated antigen presentation, supported by higher gene set scores for MHC-II (p < 0.001, t-test; Fig. 3D). Non-inflammatory TA-HEVs exhibited stress and metabolic signatures. Subcluster 5 showed higher gene set scores for endoplasmic reticulum (ER) and proteotoxic stress pathways, including heat stress and unfolded protein responses (p < 0.001 for both, t-test; Fig. 3E), and subcluster 6 showed higher gene set scores for oxidative phosphorylation and ATP biosynthesis (p < 0.001 for both, t-test; Fig. 3F). Thus, inflammatory TA-HEVs (0, 1, 12) are associated with immune-interactive transcriptional programs, while stress–metabolic TA-HEVs (5, 6) engage in stress adaptation and energy metabolism, demonstrating functional heterogeneity within TA-HEVs.

Prognostic associations of TA-HEV states across cancer types

To evaluate the potential clinical relevance of TA-HEV functional heterogeneity, we generated subcluster-specific gene signatures and applied them to The Cancer Genome Atlas (TCGA) bulk RNA-seq data. Univariate and multivariate Cox regression analyses revealed cancer type-specific prognostic differences among TA-HEV subclusters (Fig. 4A–D). In univariate analysis (Fig. 4A, B), the inflammatory TA-HEV score was significantly associated with favorable survival in SKCM (HR = 0.81, 95% CI 0.70–0.94, p = 0.004) and BRCA (HR = 0.82, 95% CI 0.68–0.96, p = 0.017), but with worse outcomes in STAD (HR = 1.20, 95% CI 1.02–1.41, p = 0.027) and PAAD (HR = 1.26, 95% CI 1.04–1.53, p = 0.019). Multivariate analysis incorporating age, sex, and tumor stage (Fig. 4C, D) confirmed independent prognostic associations in SKCM (adjusted HR = 0.756, 95% CI 0.650–0.879, p < 0.001) and KIRC (adjusted HR = 0.790, 95% CI 0.682–0.914, p = 0.002), with directional consistency maintained across 10 of 11 cancer types. In contrast, the stress-metabolic TA-HEV score showed statistically significant adverse effects in STAD (HR = 1.19, 95% CI 1.02–1.41, p = 0.028) and KIRC (HR = 1.30, 95% CI 1.11–1.53, p = 0.001) in univariate analysis (Fig. 4A, B). Multivariate analysis (Fig. 4C, D) confirmed independent adverse associations in KIRC (adjusted HR = 1.272, 95% CI 1.084–1.493, p = 0.003) and ESCA (adjusted HR = 1.410, 95% CI 1.057–1.881, p = 0.020), with directional consistency across 9 of 11 cancer types. Although several cancer types did not reach nominal significance, effect estimates were directionally adverse in most cohorts, indicating a consistent trend toward poorer survival associated with the stress-metabolic state (detailed HRs, 95% CIs and p-values in Supplementary Table S5).

Stratification of patients by high or low inflammatory and stress-metabolic TA-HEV scores revealed distinct survival patterns (Fig. 4C). In SKCM, patients with high inflammatory and low stress-metabolic scores exhibited significantly improved survival, whereas those with high stress-metabolic and low inflammatory scores showed poor outcomes. The high inflammatory TA-HEV group also exhibited an overall increase in estimated immune cell fractions, with particularly prominent enrichment of CD8⁺ T cells (Fig. 4D). In contrast, prognosis in lung squamous cell carcinoma (LUSC) was primarily driven by the stress-metabolic TA-HEV score, and a high inflammatory TA-HEV score was not accompanied by a clear increase in estimated immune cell fractions (Fig. 4E, F).

Collectively, these results define TA-HEVs as a transcriptionally and functionally coherent endothelial population.

Discussion

TA-HEVs are often described as beneficial gateways for lymphocyte entry into tumors. However, comprehensive analyses of their functional roles across cancers remain limited. Here, by integrating single-cell and bulk transcriptomic data from multiple cancer types, we delineate two coherent TA-HEV programs–an inflammatory state and a stress-metabolic state–that exhibit opposing, tumor-context-dependent associations with patient outcomes.

To identify TA-HEVs, we constrained subclustering to the tumor-associated venous endothelial lineage. According to the classical definition, HEVs are specialized post-capillary venules. Thus, our analytic strategy is consistent with the developmental and anatomical identity of TA-HEVs [ref. 11, ref. 12]. Beyond canonical HEV marker expression, we also applied complementary criteria, including signature scores, pathway enrichment, and a pseudotime analysis. This biologically coherent and dataset-agnostic framework suggests that the observed TA-HEV states reflect genuine biological change rather than artifacts. Accordingly, the cancer-specific differences in TA-HEV prevalence and state usage can be interpreted as a valid, tumor-context-dependent heterogeneity.

The inflammatory TA-HEV program encompasses innate immune stimulation, cytokine/chemokine signaling, and MHC-II antigen presentation. In SKCM, the inflammatory TA-HEV state coincides with higher immune cell infiltration, most prominently CD8⁺ T-cells, and is associated with improved survival. These observations are consistent with reports that TA-HEVs are associated with enhanced lymphocyte recruitment, correlating with clinical benefit [ref. 13]. In BRCA, STAD, and LUSC, a high inflammatory TA-HEV score was not accompanied by a clear increase in estimated immune cell fractions and was associated with adverse prognosis. These findings suggest that similar inflammatory transcriptional programs may not translate into effective immune infiltration within immunosuppressive microenvironments [ref. 14]. The stress-metabolic TA-HEV program characterized by the unfolded protein response (UPR), heat shock pathways, oxidative phosphorylation, and ATP biosynthesis is associated with poor outcomes across multiple cancers. This observation aligns with literature suggesting endothelial metabolic reprogramming and ER stress/UPR pathways can be associated with immunosuppressive vascular features and reduced therapy response [ref. 15]. Overall, TA-HEV function is not uniform but varies by state and tumor context. Although the inflammatory and stress-metabolic programs described here share features with previously reported endothelial activation and stress-response states, they were identified within a HEV-defined venous endothelial compartment. We therefore interpret them not as generic endothelial programs, but as HEV-resolved functional states that retain core HEV identity.

Our results also underscore the clinical significance of TA-HEV state subtype. Collapsing TA-HEV subclusters into state-specific TA-HEV scores provides a practical, cohort-scale means of risk stratification in cancer patients. In immunogenic tumors such as SKCM, inflammatory-high/stress-metabolic-low patients showed more favorable survivals, whereas a higher stress-metabolic TA-HEV score consistently marked high-risk patients across tumor contexts. These findings further suggest that vascular-targeted or HEV-modulating strategies should be context-tailored. These observations suggest that TA-HEV states may inform context-dependent vascular or endothelial-focused strategies, particularly in settings where stress-metabolic programs are associated with adverse immune landscape and microenvironments that may be less conducive to effective antitumor immunity.

Although scRNA-seq data limit sample numbers, we mitigated this by pan-cancer integration and by validating single-cell findings with bulk RNA-seq deconvolution. Deconvolution cannot reproduce single-cell resolution, but it offers cohort-scale readout. To our knowledge, this is the first pan-cancer definition of TA-HEV states with prognostic relevance.

Several limitations of the present study warrant consideration. The transcriptional characterization of the stress-metabolic TA-HEV subgroup is based on scRNA-seq data, which carries an inherent risk of dissociation-induced artifact. To address this, we compared the top 20 marker genes of each TA-HEV subgroup against the dissociation-induced gene reference set published by van den Brink et al. [ref. 16], and calculated artifact module scores across TA-HEV cells. The scores were similarly low across both subgroups (Fig. S5). These results suggest that the identified transcriptional programs appear to be reflective of biological states. Validation using snRNA-seq would provide additional confidence.

Additionally, the prognostic significance of TA-HEV subtype signatures was validated using TCGA bulk RNA-seq data, which is currently the most comprehensive resource for pan-cancer clinical outcome analysis at this scale. Cell-type specificity analysis supported the endothelial enrichment of both TA-HEV signatures across major tumor-resident cell types, though the degree of specificity varied between subprograms. These findings do not alter the overall prognostic associations of the TA-HEV states identified in this study, and the bulk-based validation approach remains appropriate given the current availability of pan-cancer clinical outcome data.

TA-HEV states differ transcriptionally and functionally, vary by tumor context, and are associated with opposing survival outcomes. These results propose TA-HEV states as potential biomarkers and therapeutic targets and highlight their context-dependent dual roles.

Methods

Single-cell RNA-seq data collection and preprocessing

scRNA-seq data from 317 primary tumor tissue samples (307 patients) across 11 solid cancer types were analyzed, totaling ~ 1.3 million cells after quality control. The cohorts included: stomach adenocarcinoma (31 patients, total 140,217 cells); breast invasive carcinoma (52 patients, total 347,945 cells); colorectal cancer (29 patients, total 73,777 cells); esophageal carcinoma (60 patients, total 192,078 cells); thyroid carcinoma (7 patients, total 65,580 cells); non-small cell lung cancer (50 patients, total 105,435 cells); melanoma (10 patients, total 132,464 cells); cutaneous squamous cell carcinoma (10 patients, total 26,299 cells); pancreatic adenocarcinoma (16 patients, total 42,272 cells); head and neck squamous cell carcinoma (21 patients, total 100,499 cells); and renal cell carcinoma (8 patients, total 100,916 cells). Only datasets without prior cell-type enrichment were included, while peripheral blood, metastatic, matched normal, and healthy donor samples were excluded. This study focused exclusively on tumor-associated endothelial cells, as HEVs are not present under homeostatic conditions in non-lymphoid peripheral tissues and emerge specifically in pathological contexts. The molecular characterization of normal tissue endothelial cells was therefore outside the scope of the present study. Detailed information for each dataset is provided in Supplementary Table S6.

Raw gene count matrices were processed using the Seurat R package [17]. Cells with < 200 detected genes, > 20% mitochondrial reads, or identified as doublets by DoubletFinder were excluded [18]. Gene counts were normalized with NormalizeData() (scale factor = 10,000), highly variable genes (HVGs) were identified using FindVariableFeatures(), and data were scaled with ScaleData(). Dimensionality reduction was performed by principal component analysis (PCA) with RunPCA(). A shared nearest-neighbor graph was constructed with FindNeighbors(), followed by clustering using FindClusters() with a resolution parameter of 0.8. Uniform manifold approximation and projection (UMAP) was performed with RunUMAP(), with dataset-specific parameter settings for the number of principal components (dims = 10–20). Unless otherwise specified, default parameters of the Seurat functions were applied.

Data integration and unsupervised clustering

ECs were characterized from the total cells of each dataset based on the expression of pan-EC markers (PECAM1, VWF, CDH5, and PLVAP). Count matrices of the identified ECs were then integrated and processed as described above, including normalization, identification of HVGs, scaling, and dimensionality reduction by PCA. Batch correction was performed using the Harmony R package with dims.use = 20 parameter [ref. 19].

For visualization, RunUMAP() and FindNeighbors() were applied with reduction = “harmony” and dims = 20. Unsupervised clustering was performed using FindClusters() with a resolution of 0.8. Major EC types were annotated according to well-established markers, and clusters expressing signatures of non-endothelial cell types were excluded, after which the above steps were repeated.

The venous endothelial subgroup was further reclustered at multiple resolutions, and cluster purity was evaluated using the ratio of global unshifted entropy (ROGUE) value [ref. 20]. The resolution yielding the highest ROGUE score (resolution = 1.35) was selected and applied in FindClusters(). TA-HEV subclusters were identified based on the enrichment of a previously reported 34-gene HEV signature [ref. 21, ref. 22]. The HEV signature score was calculated using AddModuleScore() in Seurat.

Differential expression, trajectory, and functional analysis

Cluster-specific DEGs were identified with Seurat FindAllMarkers() (adjusted p < 0.05). Gene Ontology (GO) enrichment was performed with clusterProfiler [ref. 23]. Differentiation states were inferred by CytoTRACE2 [ref. 24], and pseudotime trajectories constructed using Monocle3 [ref. 25]. Functional module scores, including MHC class II (MHC-II, KEGG NETWORK: N00590), Antigen processing (GO:0019886), Response to heat stress (R-HSA-3371556), Response to unfolded protein (GO:0034620), Oxidative phosphorylation (HALLMARK_OXIDATIVE_PHOSPHORYLATION), and ATP metabolic process (GO:0046034) were calculated with Seurat AddModuleScore(). The gene lists used for signature score calculation, including HEV signature genes, were provided in Supplementary Table S9. To assess potential dissociation-induced artifact in the identified TA-HEV subgroups, the top 20 marker genes of each subgroup (inflammatory: subclusters 0, 1, and 12; stress-metabolic: subclusters 5 and 6) were compared against the published dissociation-induced gene list [ref. 16], and artifact module scores were calculated across all TA-HEV cells; results are provided in Fig. S5 and Supplementary Table S11.

Bulk RNA-seq data validation

Bulk RNA-seq data from TCGA pan-cancer cohort were analyzed by single-sample gene set enrichment analysis (ssGSEA) using the top 20 upregulated DEGs from each TA-HEV subcluster. Inflammatory (subclusters 0, 1, 12) and stress-metabolic (subclusters 5, 6) TA-HEV scores were derived as the averages of ssGSEA scores from their respective subclusters. For survival analysis, these scores were standardized by Z-score transformation within each cancer type and analyzed by univariate Cox proportional hazards regression for overall survival (OS), censored at 10 years. Kaplan–Meier (KM) curves were generated after stratifying patients into four groups by the median values of both scores, and survival differences were tested using the log-rank method. To assess the independent prognostic value of TA-HEV scores, multivariate Cox proportional hazards regression was additionally performed incorporating age, sex, and stage as covariates. Cell-type specificity of TA-HEV subcluster-specific signatures was assessed by computing average expression scores for HEV across major tumor-resident cell types in publicly available scRNA-seq datasets (GSE131907, GSE132465, GSE144236, GSE148071, GSE155698, GSE161529, GSE167297, GSE172577, GSE178481, GSE183904, GSE184362, GSE188711, GSE277165); results are provided in Fig.S6.

Data visualization and statistical analysis

Gene expression heatmaps, dot plots, and box plots were generated using the SeuratExtend package, whereas UMAP plots were produced with Seurat. GO heatmaps were created with ComplexHeatmap, and bar plots as well as survival analysis plots were generated using ggplot2. Statistical analyses were performed using R 4.2.2. Statistical tests included Kruskal–Wallis, Wilcoxon, t-test, and log-rank tests, with p < 0.05 considered significant. HRs and 95% confidence intervals (CIs) were estimated using Cox proportional hazards models. Immune cell composition in TCGA samples was inferred using CIBERSORTx with the LM22 signature matrix, yielding estimated immune cell fractions for each sample [ref. 26]. All software and tools used in this study are summarized in Supplementary Table S8.

References

1. DF Quail, JA Joyce. Microenvironmental regulation of tumor progression and metastasis. Nat Med, 2013. [DOI | PubMed]
2. L Martinet. Human solid tumors contain high endothelial venules: association with T- and B-lymphocyte infiltration and favorable prognosis in breast cancer. Cancer Res, 2011. [DOI | PubMed]
3. JP Girard, TA Springer. High endothelial venules (HEVs): specialized endothelium for lymphocyte migration. Immunol Today, 1995. [DOI | PubMed]
4. D Schrama. Targeting of lymphotoxin-alpha to the tumor elicits an efficient immune response associated with induction of peripheral lymphoid-like tissue. Immunity, 2001. [DOI | PubMed]
5. P Yu. Priming of naive T cells inside tumors leads to eradication of established tumors. Nat Immunol, 2004. [DOI | PubMed]
6. P Yu. Targeting the primary tumor to generate CTL for the effective eradication of spontaneous metastases. J Immunol, 2007. [DOI | PubMed]
7. ED Thompson, HL Enriquez, YX Fu, VH Engelhard. Tumor masses support naive T cell infiltration, activation, and differentiation into effectors. J Exp Med, 2010. [DOI | PubMed]
8. B Wang. Intratumoral high endothelial venules in solid tumors: a pooled study. Front Immunol, 2024. [DOI | PubMed]
9. G Vella, Y Hua, G Bergers. High endothelial venules in cancer: regulation, function, and therapeutic implication. Cancer Cell, 2023. [DOI | PubMed]
10. S Milutinovic, J Abe, A Godkin, JV Stein, A Gallimore. The dual role of high endothelial venules in cancer progression versus immunity. Trends Cancer, 2021. [DOI | PubMed]
11. A Ager. High endothelial venules and other blood vessels: critical regulators of lymphoid organ development and function. Front Immunol, 2017. [DOI | PubMed]
12. J Amersfoort, G Eelen, P Carmeliet. Immunomodulation by endothelial cells – partnering up with the immune system?. Nat Rev Immunol, 2022. [DOI | PubMed]
13. A Asrir. Tumor-associated high endothelial venules mediate lymphocyte entry into tumors and predict response to PD-1 plus CTLA-4 combination immunotherapy. Cancer Cell, 2022. [DOI | PubMed]
14. D Klein. The tumor vascular endothelium as decision maker in cancer therapy. Front Oncol, 2018. [DOI | PubMed]
15. K Kane, D Edwards, J Chen. The influence of endothelial metabolic reprogramming on the tumor microenvironment. Oncogene, 2025. [DOI | PubMed]
16. SC van den Brink. Single-cell sequencing reveals dissociation-induced gene expression in tissue subpopulations. Nat Methods, 2017. [DOI | PubMed]
17. Y Hao. Dictionary learning for integrative, multimodal and scalable single-cell analysis. Nat Biotechnol, 2024. [DOI | PubMed]
18. CS McGinnis, LM Murrow, ZJ Gartner. DoubletFinder: doublet detection in single-cell RNA sequencing data using artificial nearest neighbors. Cell Syst, 2019. [DOI | PubMed]
19. I Korsunsky. Fast, sensitive and accurate integration of single-cell data with Harmony. Nat Methods, 2019. [DOI | PubMed]
20. B Liu. An entropy-based metric for assessing the purity of single cell populations. Nat Commun, 2020. [DOI | PubMed]
21. Y Hua. Cancer immunotherapies transition endothelial cells into HEVs that generate TCF1+ T lymphocyte niches through a feed-forward loop. Cancer Cell, 2022. [DOI | PubMed]
22. A Asrir. Tumor-associated high endothelial venules mediate lymphocyte entry into tumors and predict response to PD-1 plus CTLA-4 combination immunotherapy. Cancer Cell, 2022. [DOI | PubMed]
23. G Yu, LG Wang, Y Han, QY He. clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS, 2012. [DOI | PubMed]
24. M Kang. Improved reconstruction of single-cell developmental potential with CytoTRACE 2. Nat Methods, 2025. [DOI | PubMed]
25. J Cao. The single-cell transcriptional landscape of mammalian organogenesis. Nature, 2019. [DOI | PubMed]
26. AM Newman. Robust enumeration of cell subsets from tissue expression profiles. Nat Methods, 2015. [DOI | PubMed]