The R codes for the analyses presented in this study are available at RAGG3D/LGG_SPANK (github.com). An overview of the methods used in this study are shown in Supplementary Figure 1 .

Gene transcript-abundance and patient clinical information were collected from TCGA through the GDC Data Portal (52) and the CGGA (53–56). Progression-free survival information was used as a measure of clinical outcome (51). The cell-type specific transcriptional signatures were derived from a large collection of RNA sequencing samples spanning a wide range of cell types. For NK cells, experimentally derived dataset for IL-2 expanded (27 biological replicates), PDGF-DD activated via NKp44 signaling (4 biological replicates), and resting (25 biological replicates from 6 datasets) (41) were included. For other cell types, the data collected was from the following datasets: BLUEPRINT (57), Monaco et al. (58), ENCODE (59), Squires et al. (60), GSE77808 (61), Tong et al. (62), PRJNA339309 (63), GSE122325 (64), FANTOM5 (65), GSE125887 (66), GSE130379 (67), GSE130286 (68).

In order to derive transcriptional signatures of 21 cell types (memory B cell, naive B cell, immature dendritic myeloid cell, immature dendritic myeloid cell, endothelial, eosinophil, epithelial, fibroblast, macrophage M1 and M2, mast cell, monocyte, neutrophil, ReNK, IL2NK, SPANK, memory CD4 T cell, memory CD8 T cell, naive CD8 T cell, gamma-delta T cell and helper T cell), a total of 592 highly curated (i.e. for which identity was confirmed in the literature), non-redundant biological replicates (including 25 ReNK samples, 27 IL2NK samples and 4 SPANK samples), have been used. Due to the sparse nature of heterogeneous data sets, the expected value and variability of gene transcription abundance was inferred for each cell type using a publicly available Bayesian statistical model (github: stemangiola/cellsig), based on a negative binomial data distribution (69). This model allows to fit incomplete data (e.g. transcript abundance of one gene for which data is available in a subset of reference biological replicates) and calculate theoretical data distributions of cell-type/gene pairs. The cell-type transcriptional marker selection was based on the pairwise comparison of each cell type within cell-type categories along a cell differentiation hierarchy ( Supplementary Figure 2 ) (70). For example, all cell-type permutations from the root node of level one (including epithelial, endothelial, fibroblasts and immune cells) were interrogated in order to select the genes for which the theoretical transcript abundance distribution (data generated from the posterior distribution) was higher for one cell type compared to another. This was executed calculating the distance of the upper and lower 95% credible intervals, respectively (obtained from cellsig). From each comparison, the top 5, 10 and 20 genes per cell-type pair were selected from levels 1, 2, and 3 ( Supplementary Figure 2 ), and the union of all genes was taken as overall marker gene list. This hierarchical approach favors the identification of marker genes that distinguish broad cell-type categories as well as specific activation phenotypes.

In order to estimate the cell type relative abundance for each biological replicate, we used the algorithm CIBERSORT (71) with our RNA sequencing-derived gene marker signature. In order to estimate the clinical relevance of NK activation phenotypes (72) (73),for each cancer-type/cell-type pair, Kaplan-Meier (KM) survival curves (74) were calculated from the median split CIBERSORT-inferred proportions through the R framework tidybulk. Percent survival vs time-to-event statistics were calculated by the Log-rank (Mantel-Cox) Test (75). Statistics of KM curves were performed by log-rank test then adjusted by the Benjamini-Hochberg (BH) procedure. A table of all p-values prior to adjustment is provided in Supplementary Table 1 .

Data analysis and visualization were performed using the R environment in RStudio (76). Packages include tidyverse (77), tidybulk (73), tidyHeatmap (78), survminer (79), survival (80, 81), foreach (82), org.Hs.eg.db (83), cowplot (77), ggsci (84), GGally (85), gridExtra (86), grid (76), reshape (87), Hmisc (88), and viridis (89).

In order to visually evaluate the ability of the marker gene selection in segregating cell types, we first performed principal component analyses (PCA) (90) for three levels of cell differentiation: (i) the NK activated states, (ii) all NK cells, and (iii) all cell types. Briefly, the raw read counts were normalized by trimmed mean of M values (TMM) using tidybulk function scale_abundance, whilst CGGA transcripts were already normalized by transcripts per kilobase of exon model per million mapped reads (TPM) via RNA-seq by expectation Maximization (RSEM). PCA analyses were performed by reduce_dimensions(method=“PCA”) (73). To directly test whether the selected signature for PDGF-DD activated NK cells was suitable to accurately detect changes in cell abundance across samples with a censored time-to-event, we implemented a test on simulated data. First, for a selected number of patients N, we sampled the progression-free survival time from the clinical annotation of the LGG patient cohort. Cell type proportions were simulated using a Dirichlet distribution, according to a linear model with a slope value S and a progression-free survival as factor of interest. For each simulated dataset, the slope S was assigned to only one cell type, and the slope of 0 to all the others (i.e. only one cell type changing for each simulated tissue mixture). The intercept (baseline proportion) was defined to be the same for all cell types. The simulated proportions were used to compose the in-silico mixtures. For each simulated dataset, the transcriptional profile of each cell type was sampled at random from the reference dataset. In order to test the accuracy of our method against the presence of foreign cell types (of which transcriptomic signature was not included in our reference set), a proportion P of neural cells was added to the mixture. The framework tidybulk was used to infer the cell type proportions through CIBERSORT and perform a multiple cox regression on the predicted proportions (logit-transformed) (81), with progression-free survival censored time as a covariate. That is, for half of the N samples, the survival time was censored to half of its value. The significance calls were compared with the ground truth to generate a receiver operating characteristic (ROC) curve. For each simulation condition (values of N, S, and P), 63 test runs were performed with one variable cell type each. A range of simulation conditions were tested, ranging N from 250 to 1000, S from 0.2 to 1, and P from 0 to 0.8.

Given the innate ability of NK cells to lyse tumor cells and secrete potent anti-tumor cytokines, such as IFN-γ and TNF, prior to immunization, we hypothesized that NK cells of unique phenotype may infiltrate different cancers and confer anti-tumor immunity. Specifically, we were motivated by the recent discovery of PDGF-DD as a ligand for the activating NK cell receptor NKp44 (41) and whether this mechanism of NK cell stimulation might constitute a clinically relevant pathway of anti-tumor immunity. In order to answer this question, we gathered publicly available RNA sequencing data from 21 different immune and stromal cell types (see Methods) in order to define marker genes that may distinguish resting (91), IL-2 expanded, and PDGF-DD activated, NK cell phenotypes (41).

We next performed a principal component analysis (PCA) to determine the ability of these marker genes to segregate all three NK cell phenotypes from each other and from other cell types, such as T cells ( Figure 1A ). NK cell activation states are associated with the first principal components ( Figure 1A , left and middle panels), and NK cells overall are defined by a unique cluster when compared with other major cell types ( Figure 1A , right panel). NK cell phenotypes segregated from other cell types and from each other by PCA.

We next performed a benchmark for the inference of changes in the relative abundance of PDGF-DD activated NK cells in association with survival information for artificial tissue mixtures built from our reference data set (see Methods). This benchmark measured the ability of the PDGF-DD activated NK cell signature to extract clinically-relevant information from TCGA whole tissue RNA sequencing data ( Figure 1B ). The benchmark showed a high accuracy (area under curve) across simulation settings including magnitude of variability, sample-size, and proportion of unknown cells (please see Methods). An accuracy of 0.75 (representing the area under the ROC curve) was reached for simulation settings that match our findings on TCGA data (slope and sample size; Figures 1C, D ). We refer to the different NK cell phenotypes as: resting NK cells (ReNK), IL-2 expanded NK cells (IL2NK), and the signature of PDGF-DD activated NK cells (SPANK), respectively, and the transcript abundance of each marker gene in these NK cell phenotypes are shown in Supplementary Table 2 .

Previous studies have implicated NK cells in immune responses to glioma (15, 92–95). We next sought the association of the ReNK, IL2NK, and SPANK phenotypes with LGG. The SPANK was more abundant than either the ReNK or IL2NK phenotypes in LGG tumors from TCGA ( Figure 2A ). LGG tumors enriched for the SPANK were associated with greater overall patient survival compared to the ReNK or IL2NK phenotypes ( Figure 2B ). These results show that tumor abundance of a distinct NK cell phenotype is associated with cancer patient survival, such as LGG. Moreover, these data also suggest that LGG tumors express PDGF-DD which may activate pro-tumorigenic pathways.

In contrast to evoking NK cell anti-tumor functions through PDGF-DD binding to NKp44 (41), PDGF-DD (encoded by PDGFD) binding to PDGFR-β (encoded by PDGFRB) induces pro-tumorigenic signaling pathways that are detrimental for cancer patient survival (42, 47, 96, 97). Expression of a three-gene signature, comprised of TGFBI, IGFBP3, and CHI3L1, has previously been associated with glioma tumor cell invasion and migration and poor patient survival (98). Tumor expression of TGFBI, IGFBP3, and CHI3L1 were positively correlated with PDGFD and PDGFRB expression in TCGA LGG dataset, respectively ( Figure 3A ) (98).

Since PDGFD and PDGFRB were associated with genes involved in glioma tumor cell invasion and migration, we next examined the relationship between tumor expression of PDGFD or PDGFRB and LGG patient survival. LGG patients with low tumor expression of PDGFD had more favorable prognosis compared to LGG patients with high tumor expression of PDGFD ( Figure 3B ). Higher LGG tumor expression levels of PDGFRB alone displayed a trend towards poor survival ( Figure 3B ) when PDGFD expression was low, but this was not statistically significant ( Figure 3C ). These data show that high tumor expression of PDGFD is primarily associated with poor prognosis compared to PDGFRB expression in TCGA LGG dataset.

Anti-tumor immunity would be expected to curtail pro-tumorigenic factors and benefit patient survival. In addition to pro-tumor functions, we hypothesized that tumors enriched for the SPANK would contribute to anti-tumor immunity resulting in a more favorable TCGA LGG prognosis (41). In order to assess whether the abundance of these NK cell phenotypes counteracted the pro-tumorigenic expression of PDGFD and thus improve TCGA LGG prognosis, we next determined the progression-free survival of LGG patients with tumors stratified for PDGFD expression and abundance of either the ReNK, IL2NK, or SPANK phenotypes. When LGG tumors were stratified for PDGFD expression, patients with tumors enriched for the SPANK had a more favorable prognosis compared to LGG patients with a lower tumor abundance of SPANK (e.g. compare H^PDGFD/H^SPANK, red KM curve, to H^PDGFD/L^SPANK, grey KM curve and compare L^PDGFD/H^SPANK, yellow curve, to L^PDGFD/L^SPANK, blue curve) ( Figure 4 ). In contrast, this was not observed for either ReNK or IL2NK ( Figure 4 ). These results show that LGG tumors enriched for the SPANK may mitigate the detrimental effect of PDGFD expression on the prognosis of TCGA LGG patient cohort compared to the ReNK or IL2NK phenotypes.

For a given cancer, it is likely that immune subsets other than NK cells infiltrate the tumor microenvironment to elicit anti-tumor immunity, particularly T cells. We were interested in knowing whether the abundance of a given T cell subset in the LGG tumor microenvironment is associated with anti-tumor immunity. TCGA LGG tumors enriched for the memory CD8⁺ T cell phenotype were associated with improved prognosis, but not the naïve, γδ, CD4⁺ memory, or Helper, T cell phenotypes ( Figure 5A ). Our analyses show that TCGA LGG patients enriched for the memory CD8⁺ T cell phenotype have improved prognosis.

In order to determine whether each T cell phenotype might counteract the detrimental expression of PDGFD on LGG prognosis, similarly to the SPANK, we next determined progression-free survival for LGG patients with tumors stratified for expression of PDGFD and the abundance of each T cell phenotype, respectively ( Figure 5B ). Using this approach, LGG tumors enriched for the memory CD8⁺ T cell phenotype were associated with improved prognosis ( Figure 5B ). These results show that LGG tumors enriched for the memory CD8⁺ T cell phenotype mitigate the pro-tumorigenic effects of PDGFD because they are associated with improved prognosis.

Given that the abundance of NK and T cell phenotypes differ markedly in LGG tumors ( Supplementary Figure 3 ), we were interested in understanding the relative contribution of SPANK and T cell phenotypes for LGG prognosis, respectively. We therefore determined patient survival for TCGA LGG tumors stratified for the abundance of SPANK and each respective T cell subset ( Figure 5C ). Interestingly, LGG tumors enriched for the SPANK and CD4⁺ T helper phenotypes (TH) had improved prognosis compared to other strata e.g. compare H^SPANK/H^TH to either H^SPANK/L^TH or L^SPANK/H^TH or L^SPANK/L^TH ( Figure 5C , column 1). In contrast, LGG tumors enriched for the SPANK and memory CD8⁺ T cell phenotypes (CD8mem) did not further improve LGG patient survival compared to other strata e.g. compare H^SPANK/H^CD8mem to either L^SPANK/H^CD8mem or H^SPANK/L^CD8mem ( Figure 5C , column 5). We conclude that LGG tumors enriched for the SPANK and CD4⁺ T helper cell phenotypes are associated with improved LGG prognosis in TCGA. These results provide new insights into the possible cooperation between different NK and T cell subsets for LGG anti-tumor immunity which may inform adoptive cell therapies.

NK cells express a family of germline-encoded activating and inhibitory surface receptors that engage in cancer immune surveillance, which can also be expressed by memory CD8⁺ T cells. However, the NK cell family receptors most critical for anti-tumor immunity in LGG remain unclear. Given that LGG patients in TCGA with tumors enriched for the SPANK or memory CD8⁺ T cells were associated with improved prognosis, we were interested in analyzing whether tumor expression of transcripts encoding NK cell family receptors was also associated with improved prognosis for TCGA LGG patients. LGG tumors with high expression of the KLRK1, KLRC1, KLRC2, KLRC3, or KLRC4 transcripts encoding the NKG2D, NKG2A, NKG2C, NKG2E, and NKG2F NK cell receptors, respectively, were associated with improved prognosis ( Figure 6A ). In contrast, high LGG tumor expression of CD226, CD244, CRTAM, KIR2DL4, NCR1, or NCR3 encoding the DNAM-1, 2B4, CRTAM, NKp46 and NKp30 NK cell receptors, respectively, were not associated with prognosis ( Supplementary Figure 4 ). Moreover, expression of the KLRK1, KLRC1, KLRC2, KLRC3, and KLRC4 receptor genes were overwhelmingly positively correlated with the SPANK and memory CD8⁺ T cell phenotypes in TCGA LGG tumors ( Figure 6B ). These results show that high expression of the Killer cell lectin-like receptor (KLR) family in LGG tumors is associated with improved prognosis, suggesting that expression of KLR receptors may be critical for regulating anti-tumor immunity in LGG.

Since LGG tumors enriched for expression of the SPANK and KLR family receptors were associated with improved prognosis in the TCGA LGG patient cohort, we next sought to validate these findings using another glioma patient dataset, such as the CGGA ( Figure 7 ). Similar to TCGA LGG patient cohort, high tumor expression of PDGFD in CGGA LGG patients was associated with a poor prognosis ( Figure 7A ). Moreover, LGG tumors enriched for the SPANK and memory CD8⁺ T cell phenotypes were associated with improved prognosis when CGGA tumors were also stratified for PDGFD expression compared to the IL2NK phenotype ( Figure 7B ). In contrast to TCGA LGG dataset, high tumor expression of KLRC1 and KLRC2 in CGGA LGG patients were associated with improved prognosis, but not KLRC3, KLRC4 or KLRK1 ( Figure 7C ). Like TCGA, high expression of KLRC2, which encodes the activating NKG2C receptor, was also associated with the SPANK and memory CD8⁺ T cell phenotypes in LGG tumors ( Figure 7D ). However, in contrast to TCGA, expression of KLRC1, which encodes the inhibitory NKG2A receptor, was associated with the ReNK phenotype and not the SPANK or memory CD8⁺ T cell phenotypes in CGGA LGG tumors ( Figure 7D ). Similar to TCGA, these results show that high tumor expression of PDGFD is associated with poor CGGA LGG prognosis, and tumors enriched for the SPANK and memory CD8⁺ T cell phenotypes have improved prognosis. Moreover, like TCGA, high LGG tumor expression of KLRC2 is also associated with the SPANK and memory CD8⁺ T cell phenotypes and improved prognosis of CGGA LGG patients.