In this study, we used two cohorts to explore novel biomarkers in pituitary tumors. The first dataset included 112 patients, and another independent test dataset contained 73 patients. Tumor specimens were obtained from patients with NF-PitNEts (n=112) who underwent transsphenoidal surgical resection at Beijing Tiantan Hospital between June 2018 and July 2019. The diagnosis of NF-PitNEts is defined as the absence of clinical and biochemical evidence of overproduction of adenohypophysis hormone. The mean age of these 112 patients was 52 years (range, 21-75), and there were 63 males and 49 females. The demographics and clinicopathological features of the patients are summarized in Table 1. Tumor cavernous sinus (CS) invasion was defined as Knosp grade 3 and 4 or intraoperative evidence (23). The expression profiles and matching clinical information of independent test datasets (n=73) were described previously (24). In addition, 16 NF-PitNEts specimens (8 invasive and 8 noninvasive) were collected from the same hospital as an independent validation cohort (Table S3), and their expression levels were verified by quantitative real-time polymerase chain reaction (qRT-PCR). This study recruitment process and protocol were approved by the Medical Ethics Committee of Beijing Tiantan Hospital, and informed consent was obtained from all individual participants.

A total of 1-3 μg RNA per sample was extracted and purified from the collected specimens of NF-PitNEts. According to the instructions provided, sequencing libraries were constructed using the NEBNext^® Ultra™ RNA Library Prep Kit for Illumina^® (#E7530L, NEB, USA). After the library was successfully generated (effective concentration >10 nM), the index-coded samples were clustered on the cBot cluster generation system using HiSeq PE Cluster Kit v4-cBot-HS (Illumina). The library was then sequenced on an Illumina platform, and 150-bp paired-end reads were generated. Raw data were filtered with FAST-QC, and the clean reads were then mapped to the human genome hg19 sequence (GRCh37) using HISAT2 (25). HTseq was used to generate gene counts, and the RPKM method was used to determine gene expression (26).

The Estimation of Stromal and Immune cells in Malignant Tumor tissues using Expression data (ESTIMATE) algorithm (27) was used to obtain the immune levels of 112 pituitary tumor patients. Based on the median immune score, patients were divided into high and low groups. These 112 patients were also divided into invasion and noninvasion groups based on their clinical invasion information.

The “limma” R package was used to obtain the DEGs from the high vs. low immune score groups and invasion vs. noninvasion groups. Genes with an adjusted P value < 0.05 and |log2-fold change| ≥ 0.585 were filtered as DEGs.

First, the ImmuCellAI database was applied to estimate the abundance of 24 immune cell types in the 112 patients with NF-PitNEts (28). The abundance of these immune cells in the invasion and noninvasion groups was then used for WGCNA.

Second, the “WGCNA” R package was applied to build a coexpression network of immune cell-related genes. Next, several gene modules were detected based on their similar expression patterns. Finally, the abundance of immune cells was associated with these gene expression modules, and genes in the modules were selected for further analysis.

The “randomForest” R package was used to select the candidate genes associated with invasion. We first divided the 112 patients into a training dataset (n=78) and a self-test (n=34) dataset based on the “sample” function of the R program. Then, random forest analysis (RF) was performed through the “randomForest” function in the training dataset, the error rate curve was drawn and changes in the error rate of different numbers of genes selected were observed. Finally, the “ggplot” R package was used to show the MeanDecreaseAccuracy and the best RF model of these genes.

Genes with the top 20 MeanDecreaseAccuracy were used to build the SVM model through the “e1071” R package, and the “pROC” package was used to perform the classification efficiency of the model in the training dataset, self-test dataset and independent test dataset.

The R package “ConsensusClusterPlus” was applied to explore the classification efficiency of 8 crucial genes. Subsequently, the “pROC” package was used to show the classification efficiency of these genes.

Functional enrichment analysis was carried out by the “clusterProfiler” R package, and only functions with a p value<0.05 were selected.

The QIAGEN RNeasy Kit and High Capacity cDNA Reverse Transcription Kit were used to extract total RNA from the verified samples and conduct reverse transcription reactions. qRT-PCR was performed in a volume of 20 µl with Power SYBR™ Green PCR Master Mix on a QuantStudio 3 and 5 System (Applied Biosystems). GAPDH was used as a housekeeping control. The sequences of the primers are shown in Table 2.

All statistical analyses were performed by R software (version 3.6.3). The bar plot between two different groups was drawn by the “ggplot” package, and a T test was used to compare the differences. A P value <0.05 was considered significant.

ESTIMATE analysis was performed in the 112 NF-PitNEts expression data, and the samples were divided into high- and low-immune score groups. Next, 3152 DEGs were obtained from the high- vs. low-immune score groups, including 662 upregulated and 2490 downregulated genes (Figure 1A and Table S1). The enrichment analysis indicated that these DEGs were involved in Th1 and Th2 cell differentiation, the relaxin signaling pathway, and the FoxO signaling pathway (Figure 1B). Then, we obtained 525 DEGs (352 upregulated and 173 downregulated) associated with invasion (Figure 1C and Table S2). These genes are related to the FoxO signaling pathway and Th17 cell differentiation (Figure 1D). This finding reveals that invasion-associated DEGs are involved in the immune-related functions of pituitary tumors.

Subsequently, we predicted the immune cell abundance of these patients. Cytotoxic T cells (Tc), Th2 cells, natural killer T cells (NKT), dendritic cells (DC), B cells, monocytes and neutrophils differed between the invasion and noninvasion groups (Figure 1E).

We used the “WGCNA” package to calculate the immune cell-associated modules in pituitary tumors. The soft power of the coexpression network was calculated through the “pickSoftThreshold” function and was set to 20 (Figure 2A). At this power, the R² of the scale-free topology model under the soft threshold was 0.92, which indicates that the network conformed to the scale-free feature (Figure 2A). Then, the network was constructed, and seven coexpression modules (green, turquoise, blue, yellow, black, brown, and red) were built (Figure 2B). Subsequently, the correlation between these seven modules and immune cell abundance was calculated. The green module was significantly correlated with neutrophils; the turquoise module was significantly correlated with monocytes; the blue module was significantly correlated with cytotoxic cells; and the yellow, black, brown, and red modules were significantly associated with NKT cells (Figure 2C). These results suggest that these genes are immune-related in the tumor microenvironment of pituitary tumors.

To select the candidate crucial genes in pituitary tumors, we first set the intersection of immune-related DEGs, invasion-related DEGs and module genes and found 81 common genes (Figure 3A). Next, we used the 81 common genes to perform RF analysis. We started the RF analysis by generating 1000 decision trees, which showed a lower error rate (Figure 3B). Then, the random forest results indicated that, when the number of candidate genes was 20, the classification efficiency error rate of the model was the lowest. We showed the top 50 accuracies of these genes and selected the top 20 for further analysis (Figures 3C, D).

Subsequently, SVM analysis was used to construct the IICM (Figure 4A). We found that this model could distinguish the invasive and noninvasive patients in the training dataset (AUC=0.732, Figure 4A). It also exhibited good efficacy in the self-test dataset (AUC=0.653, Figure 4B) and independent test dataset (AUC=0.619, Figure 4C). The results showed that this model could be used to predict the invasion state of pituitary tumors.

We then analyzed the function of these 20 key genes. The results suggested that 8 of the 20 genes were enriched in multiple signaling pathways, such as the KRAS signaling pathway and PPAR signaling pathway (Figure 5A). This result suggests that the function of this model is mainly driven by these 8 genes (BMP6, CIB2, FABP5, HOMER2, MAML3, NIN, PRKG2 and SIDT2).

All 8 genes were differentially expressed in the invasion and noninvasion groups (Figure 5B), and most were differentially expressed in the high- and low-immune score groups (Figure 6). Subsequently, the consensus clustering analysis suggested that the 8 genes could well divide the patients into two groups (Figures 7A, B). We then constructed a classifier for these 8 genes, and the results showed that the 8-gene classification model showed good classification efficiency in pituitary tumors (AUC=0.671, Figure 7C). Then, qRT-PCR was performed to validate the expression level of these genes between invasive and noninvasive NF-PitNEts (Figure 7D). Consistent with the sequencing data, the mRNA expression of BMP6, CIB2, HOMER2, MAML3, NIN, PRKG2, and SIDT2 was significantly upregulated in invasive NF-PitNEts, while FABP5 was significantly downregulated. Overall, we filtered 8 genes that could predict the invasion status of pituitary tumors and could be used as predictors for further treatment of the tumors.