Brain tumor transcriptomes of 3810 patients were downloaded from the Genomic Data Commons, the NCBI Gene Expression Omnibus, the International Cancer Genome Consortium, the Chinese Glioma Genome Atlas, the European Bioinformatics Institute, and GlioVis ( Table S1 ). Each dataset has more than 100 patients accompanied by their vital status and period of follow-up. Histology, tumor grade, radio-chemotherapy treatment, recurrent/secondary status at sampling time, IDH mutation, MGMT methylation, and co-deletion of 1p and 19q data were collected if available ( Table S1 ). Meanwhile, we collected 93,293 single-cell RNA profiles subjected to Smart-Seq2 sequencing protocol from 16 previous studies ( Table S2 ) from brain cancer, lung cancer, colorectal cancer, ovarian cancer, melanoma, and head and neck squamous cell carcinoma. These single cell datasets encompass T cells, B cells, monocytes, macrophages, natural killer cells, dendritic cells, cancer cells, and other nonmalignant cells including fibroblast, epithelia cells, gliocytes, and neurons. In addition, we manually curated a list of genes related to tumor microenvironment, immune cells, immune checkpoint blockade therapy response, and prognosis ( Table S3 ).

For each single-cell dataset, we performed logarithmic transformation as log₂(TPM/10 + 1). We clipped gene expression values at 99% quantile values of all genes. Subsequently, we employed R package preprocessCore (version 1.40.0) (12) to perform quantile normalization. We applied ComBat routine implemented in R package sva (version 3.26.0) (13) to perform batch effect correction for the normalized expression data of bulk brain tumors.

We developed a feature encoder through self-supervised feature representation learning. The feature encoder could learn nonlinear feature representations of transcriptomes in a reduced dimensional space.

The feature encoder was trained with a self-supervised deep learning algorithm based on contrastive learning (14). Specifically, contrastive learning allows the feature encoder to learn representations in a label-free manner ( Figure S1 ). Positive pairs are defined as two different noise-adding views of the same transcriptome (V_q , V_k+ ). Two different transcriptomes form a negative pair (V_q , V_k- ). V_k- came from a dictionary of transcriptomes {V_k1-, V_k2-,…, V_n1-}, which was defined on-the-fly by a set of trained data. Contrastive loss aims to minimize the distance between the positive pair and maximize the distance between the negative pair. The function of contrastive loss is defined as:

where τ is a temperature hyper-parameter (15). In this analysis, the similarity of each pair was calculated based on the expression features extracted from the feature encoder. In this manner, the feature encoder was driven to learn features of transcriptomes by contrastive loss.

In our task, the feature encoder was trained to learn the same representation of different noise-adding views of the same single-cell transcriptome and dissimilar representation of different cells.

The 93,293 single cells were randomly divided into a training set (N = 83,964) and a validation set (N = 9,329). We logarithmically transformed the transcriptomes of the preselected genes ( Table S3 ) and scaled to a range of 0 to 1 before feeding them into the feature encoder. At each epoch, we made noise through random zero out and shuffling and added Gaussian noise (mean: 0, standard deviation: 0.1) to 20% of genes for all transcriptomes.

The feature encoder was an 18 layered ResNet (16). We replaced the convolutional layer in the original ResNet with a linear layer to allow it to process gene expression data. For each residual block, the input skips training from a few layers and is connected directly to the output. Moreover, we set the project head with 128 output neurons. The use of multi-layer perceptron (MLP) as project head was demonstrated to be beneficial for contrastive learning. The architecture of the feature encoder was provided as Figure S2 .

We employed stochastic gradient descent algorithm (17) as the optimizer. The weight decay of the optimizer is 1e-4 and the momentum is 0.9. We set batch size for each training iteration of 256. The initial learning rate was 0.03 and decay with a cosine annealing schedule. We set the contrastive learning dictionary size to 3072. The momentum and τ of contrastive loss were set to 0.999 and 0.2, respectively. The model was trained in parallel on two graphic processing units for 300 epochs. The model was developed with PyTorch (v1.3.0) package.

The developed feature encoder was applied to extract expression signatures from bulk sample transcriptomes. Specifically, we extracted feature representations from the developed feature encoder applied to the expression data of TCGA pan-cancer. The feature encoder transformed the expression profile of each bulk sample into 128 features, which was determined by the output neurons of feature encoder. The extracted features were hierarchically clustered through R package ConsensusClusterPlus (version 1.42.0) (18). The obtained clusters were further grouped into expression signatures because of the high negative correlations among these clusters ( Figure S3 ). Subsequently, we dichotomized brain tumor patients based on each expression signature and selected one that can better represent unique immune infiltration signature of brain tumor. Specifically, we used R package Ckmeans.1d.dp (version 4.2.1) (19) to perform k-means clustering. The k-means clustering cutoff value closest to the median value of signature was selected as the optimal cutoff to dichotomize samples. We used R package fgsea (version 1.6.1) (20) to perform gene set enrichment analysis (GSEA) for a gene set related to unique immune infiltration properties of brain tumor such as microglia and reactive gliosis ( Table S4 ). We kept the signature that ranked on the top to dichotomize patients as mentioned above for downstream analysis. A flowchart illustrating this procedure was provided in Figure S3 .

To examine the advantage of self-supervised learning paradigm, we employed principle component analysis (PCA) as a linear feature encoder and compared the PCA features with the deep learning features. Specifically, we performed PCA on single cell transcriptomes of the 2616 filtered genes with python package sklearn (v0.24.1). The principle components of single cells were then projected to brain tumor transcriptome. Then, brain tumor patients were dichotomized through hierarchical clustering based on R package ConsensusClusterPlus (version 1.42.0) (18).

We analyzed the association between molecular subtypes with immune and genomic alteration signatures, which include immune cellular fractions, immunomodulatory expressions, oncogenic and immune pathways, genomic alterations, driver mutations, and molecular subtypes of glioblastoma proposed by the Cancer Genome Atlas (TCGA) (21). We used CIBERSORT to estimate the proportions of 22 immune cell types based on LM22 matrix (22). We performed paired t-test for 78 genes related to immunomodulation (23) in the 11 collected datasets. In addition, we performed GSEA based on R pakage fgsea (version 1.6.1) (16) for cancer hallmark (24) and immune-related gene sets (19) ( Table S4 ). Continuous variables were evaluated by Wilcoxon rank sum test, while discrete variables were evaluated via Chi-square test if not specified. For GSEA, P-values were calculated based on 10,000 permutations. Kaplan-Meier survival analysis and multivariate Cox hazards model were utilized to analyze the association of subtypes and prognosis, which were carried out with R survival package (2.40-3).

As described above, we generally employed Wilcoxon rank sum test or Chi-squared test for the statistical analysis as appropriated, if unspecified. The paired t-test was used for the analysis of immunomodulators. The P value of enrichment analysis were calculated based on 10,000 permutations. We employed Kaplan-Meier analysis to estimate survival distribution. Cox proportional-hazards model was utilized for multi-variable survival analysis. We applied log-rank test to compare the statistic difference of survival curves between two groups. All figures and statistical analysis were conducted using R software (version 3.6.1). A P < 0.05 was considered as statistically significant. All statistical tests were two-sided. P-values were adjusted with FDR method.

RNA profiles and clinical data of 3810 brain tumor patients were collected from 11 public studies. The baseline characteristics are shown in Table S5 . Glioma and medulloblastoma account for 83% (3146) and 16% (624), respectively. Among 2951 patients with cancer type information, the primary, recurrent, secondary, and post-treatment tumor account for 71% (2105), 11% (332), 1% (40), and 16% (474), respectively. In the glioma cohort (3146), patients with tumor grade II, III, and IV respectively accounted for 24% (743), 28% (872), and 39% (1226); 10% (305) of patients did not have tumor grade information. Meanwhile, there were 2960 glioma patients with pathological information. This glioma cohort consisted of diverse pathological subtypes such as astrocytoma (27%), oligodendroglioma (21%), oligoastrocytoma (8%), and glioblastoma (44%). Among 2289 glioma patients with IDH mutation examined, 55% of them (1254) had IDH mutation. Among 1098 of these 1254 patients with co-deletion of 1p/19q tested, 42% (459) carried co-deletion of 1p and 19q. Among 1731 patients tested for MGMT methylation, 59% (1028) of them were positive for hypermethylation of MGMT promoter. Among 1226 patients with radio-chemotherapy treatment information, the proportion of patients treated with chemotherapy, radiotherapy, and a combination of both were 15% (184), 23% (281), and 62% (761), respectively.

A flowchart depicting the whole procedures of this study is shown in Figure 1 . We collected 93,293 single-cell RNA profiles from 16 published datasets. A total of 2,616 genes were selected for the analysis. These genes were associated with tumor microenvironment, immune cells, immune checkpoint blockade therapy response, and prognosis (See Methods). We developed a self-supervised deep learning model based on single-cell RNA profiles of these 2616 genes to decipher gene expression signatures from transcriptomes. Subsequently, we applied this developed feature encoder to extract expression signatures from transcriptome of bulk brain tumor samples (See Methods and Figure S3 ). We then examined the association of expression signature with immune signatures, genomic alteration, and prognosis.

The results obtained from CIBERSORT (22) showed that 18 of 22 types of immune cells were significantly different between C1/2 subtypes ( Figure 2A ; Wilcoxon rank sum test, P < 0.05). All types of TAMs (i.e. M0, M1, M2), CD4+ follicular helper T cells, and neutrophils had higher infiltration rate in C2 as compared with C1 subtype ( Figure 2A ). Contrastively, C1 had higher infiltration of CD8+ T cells, plasma cells, and dendritic cells than C2 subtype ( Figure 2A ). The infiltration of the other cell types was provided in Table S6 . Furthermore, the immune infiltration of C1 versus C2 exhibited consistent trends among different brain tumor subtypes ( Figure S4 ). For example, there are higher infiltration of plasma cells and lower infiltration of M1 and M2 macrophages in C1 subtype ( Figure S4 ).

We observed that 24 immunomodulatory genes were differentially expressed in C1 versus C2 subtype ( Figure 2B ; Paired t-test, P <0.05). Specifically, C10orf54, CX3CL1, and EDNRB were highly expressed in C1 versus C2 subtype ( Figure 2B ). CD276, CCL5, CXCL10, and HMGB1 and the other 17 immunomodulatory genes were significantly upregulated in C2 versus C1 subtype ( Figure 2B ). The detailed expression of all immunomodulatory genes were provided in Table S7 .

Enrichment analysis of 50 cancer hallmarks and 132 immune signaling modules showed that CSF-1, MYC, TGF-β, JAK/STAT3, IFN-α, and the other 29 signaling pathways were enriched in C2 versus C1 subtype ( Figure S5 , P < 0.05). For 11 signaling modules, including core serum response, proliferation, DNA repair, and E2F target pathways, the same trends were validated across brain tumor subtypes ( Figure S5 , P < 0.05).

In the TCGA low-grade glioma, non-silent mutation burden, intratumor heterogeneity, aneuploidy, and the other six types of genomic variation were significantly higher in C2 versus C1 subtype ( Figure 2C and Table S8 ; Wilcoxon rank sum test, P < 0.05). The corresponding trends were observed in astrocytoma and oliodendroglioma patients ( Figure S6 ). In the TCGA glioblastoma cohort, there was no difference among the aforementioned variations except for segments of copy number variation ( Figure S6 ).

We also examined the association of C1/2 subtypes and driver gene mutations of brain tumors that linked to prognosis and therapeutic resistance ( Table S6 ). Our finding showed that four driver events were significantly higher in C1 versus C2 subtype, including IDH mutation, hypermethylation of MGMT promoter, high CpG island methylation phenotype (G-CIMP), and co-deletion of 1p and 19q ( Figure 2D ; Chi-squared test, P < 0.05). Four driver events were significantly higher in C2 versus C1 subtype such as EGFR amplification, deletion of CDKN2A/CDKN2B and PTEN, gain of chromosome 7, and/or loss of chromosome 10 ( Figure 2D ; Chi-squared test, P < 0.05). Across brain tumor subtypes, such as glioblastoma and low-grade glioma, differences in mutation rates showed the same trend in eight driver events among C1/2 subtypes, including CDKN2A/CDKN2B and PTEN deletion, IDH mutation, and co-deletion of 1p and 19q ( Figure S7 ).

In addition, we found that C1/2 subtypes were linked to TCGA molecular subtypes, namely classical, neural, proneural, and mesenchymal subtypes (21) ( Figure 2E ). Neural [168(37%) versus 104(11%); Chi-squared test, P = 6.2e-29] and proneural subtypes [186(41%) versus 243(26%); Chi-squared test, P = 6.2e-29] were significantly enriched in C1 versus C2 subtype. C2 had higher proportions of classical [307(33%) versus 56(12%); Chi-squared test, P = 1.4e-16] and mesenchymal subtypes [265(29%) versus 44(10%); Chi-squared test, P = 2.3e-15] as compared with C1 subtype.

Clinical characteristics of brain tumor patients were provided in Table S9 . C2 subtype had lower Karnofsky scores (Median: 80 vs. 90, Wilcoxon rank sum test, P = 3.4e-6) and higher tumor microvascular infiltration rate versus C1 subtype (61/76, 80% vs. 31/65, 48%; OR: 4.2, 95% CI: 2.0 – 8.7; Chi-squared test, P = 1.8e-4). Among patients with recurrence, C1 subtype had marginally significant lower distant recurrence rate (4/23, 17% vs. 19/48, 40%; OR: 0.3, 95% CI: 0.1 – 1.1) and higher local recurrence rate (19/23, 83% vs. 29/48, 60%; OR: 3.1, 95% CI: 0.9 – 10.6) as compared with C2 subtype (Chi-squared test, P = 0.1). There were no significant differences in family history of cancer, pre-diagnostic symptoms, and tumor location between C1/2 subtypes (Chi-squared test, P > 0.5).

Kaplan-Meier survival analysis showed that C1 subtype has better survivability than C2 subtype ( Figure S8A ; Log-rank test, P = 8.2e-78) in the combined cohort of 3810 patients. This result was also observed in each individual in the 11 datasets ( Figure S8A ; Log-rank test, P < 0.05). Moreover, the difference remained significant in the combined cohort after controlling for confounding factors such as age, gender, tumor, histology, radio-chemotherapy, recurrent/secondary status, IDH mutation status, MGMT methylation status, and co-deletion of 1p and 19q ( Figures S8B , S12 ; Multivariate Cox model, HR: 2.2, 95% CI: 1.7 – 2.9; Log-rank test, P = 3.7e-10). The independent association of C1/2 subtypes with prognosis from the multivariate model remained significant in six individual datasets and exhibited the same trend in the other four datasets ( Figures S8B , S12 ; Log-rank test, P < 0.05). In the TCGA glioma cohort, surgery type was taken into consideration additively. In the medulloblastoma cohort (i.e. GSE85217), clinically relevant confounding factors, such as age, gender, and molecular subtypes, were included.

We observed that C1/2 subtypes of PCA have significantly different overall survival in seven independent datasets ( Figure S9 ; Log-rank test, P < 0.05). Cox analysis ( Table S10 ) shows that C1/2 subtypes have prognosis significance in five individual datasets (Log-rank test, P < 0.05) and did not show any trend in three datasets (i.e. E-MTAB-3892, TCGA-GBM, GSE13041). In summary, the association between prognosis and expression signatures derived from deep learning is more general as compared with PCA.

We examined the association between C1/2 subtypes and prognosis of glioma patients with respect to histology, genomic alteration, and grade. The glioma patients were divided into nine subgroups: astrocytoma, oligodendroglioma, glioma with or without IDH mutation, glioma with IDH mutation with or without co-deletion of 1p and 19q, tumor grade II, III, and IV ( Figure 3A ). The C2 subtype has significantly poorer survival outcome than C1 in all subgroups ( Figure 3A ; Log-rank test, P < 0.05). In addition, the difference remained significant in eight out of these nine subgroups and marginally significant in grade IV glioma after taking into account age, gender, histology, IDH mutation status, MGMT methylation status, and co-deletion of 1p and 19q ( Figures 3B , S12 ; Log-rank test, P < 0.05). The dataset was taken as strata variable in multivariate Cox model.

The C2 subtype of glioblastoma with IDH mutation has poor survival outcome analogous to glioblastoma without IDH mutation ( Figure 4A ; Log-rank test, P = 0.8). The C1 subtype of glioblastoma with IDH mutation, meanwhile, has a favorable survival outcome versus C2 subtype (Log-rank test, P = 1.2e-3) or glioblastoma without IDH mutation (Log-rank test, P = 1.3e-6). The result remained significant after ruling out confounding impacts of age, gender, and co-deletion of 1p and 19q ( Figure 4B ). CIBERSORT analysis demonstrated that there are high infiltration rates of regulatory T cells and dendritic cells and low infiltration rates of follicular helper T cells, M1 macrophages, and neutrophils in C1 subtype of glioblastoma with IDH mutation ( Figure 4C ; Wilcoxon rank sum test, P < 0.05). GSEA analysis showed that glycolysis, MTORC1, core serum response, proliferation, and E2F signaling pathways were enriched in the C2 subtype of glioblastoma with IDH mutation and IDH wildtype glioblastoma as compared with C1 subtype of glioblastoma with IDH mutation ( Figure 4D ; P < 0.05).

We further stratified glioblastoma patients through different treatment modalities. For patients with radiotherapy (N = 86), Kaplan-Meier survival analysis showed that C2 subtype has better overall survival than C1 subtype ( Figure S10 , Log-rank test, P = 3.3e-2). The same trend was also observed in C1/2 subgroup of patients with radio-chemotherapy (N = 49, Figure S10 , Log-rank test, P = 0.2).

Kaplan-Meier survival analysis showed that C2 subtype had worse progression-free survival as compared with C1 subtype in TCGA glioma cohort ( Figure S11 , Log-rank test, P = 6.1e-4). The difference remained significant in radio-chemotherapy patients ( Figure S11 , Log-rank test, P = 5.3e-3) and showed the same trend in radiotherapy alone patients ( Figure S11 , Log-rank test, P = 0.4). Progression-free survival was not analyzed for the chemotherapy group due to the limited sample size ( Table S11 ).