The transcriptome data and clinical information about PAAD and 32 other cancer types were downloaded from the UCSC website (https://xena.ucsc.edu/). The PAAD data set from UCSC was selected as the training set. We select ICGC-PAAD (N= 182) and GSE15471 (N = 36) as the validation set. ICGC-PAAD (N = 182) included 182 pancreatic cancer samples from different patients. GSE15471 (N = 36) included 36 pairs of pancreatic cancer and normal tissue samples. The full clinical information about these two cohorts is shown in Supplementary Table 4 .

Based on previous studies, we obtained a list of eRNAs and its target genes predicted by the PresSTIGE method (16, 17). In the training set, the patients’ eRNA transcriptional data and clinical data were matched together. We used R software package “survival” and “survminer” to investigate the relationship between eRNAs and OS in PAAD patients. The correlation level of eRNAs and its target gene was evaluated through co-expression analysis methods. The most significant survival-associated eRNA and its target gene were selected for further analysis.

Based on the website Human Protein Atlas (HPA, http://www.proteinatlas.org/) (18, 19), we explored the mRNA and protein levels of PHF10 in various cancer tissues and normal tissues. Meanwhile, we downloaded the immunohistochemical images of PHF10 protein in pancreatic cancer from this website. The differential mRNA expression level of LINC00242 and PHF10 between PAAD tissues and normal tissues were analyzed via the Gene Expression Profiling Interactive Analysis (GEPIA) web tool (http://gepia.cancer-pku.cn/) (20). GSE15471 was selected to verify the difference in the expression of LINC00242 and PHF10 between PAAD tissues and normal tissues. Through univariate and multivariate Cox regression analyses, we determined the independent prognostic value of LINC00242. The correlation between LINC00242 expression level and clinicopathological was explored. ICGC-PAAD was selected to verify the relationship between LINC00242 expression and clinical parameter. Besides, we selected the differentially expressed genes in tumor compared with normal tissues and co-expressed genes to do the GO and KEGG analyses. Next, the prognostic value of LINC00242 and its co-expression relationship with PHF10 were verified in TCGA (32 other types of cancer) pan-cancer data.

The human pancreatic cancer cell line (PANC-1) was purchased from the Cell Bank of the Chinese Academy of Sciences (Shanghai, China). We use DMEM medium supplemented with 10% fetal bovine serum (Thermo Fisher Scientific, USA) and 1% penicillin–streptomycin to culture the cells at 37°C, 5% CO2 atmosphere. LINC00242 siRNAs were synthesized by GeneChem (Shanghai, China), and Lipofectamine 2000 (Invitrogen, California, USA) was employed to achieve the cell transfection.

Differently treated PAAD cell lines were collected in the logarithmic growth phase and mixed with 1 ml TRIzol. Chloroform was added and kept at room temperature for 15 min, then isopropanol was added to the solution and centrifuged to obtain RNA precipitate. qRT-PCR analyses of LINC00242 and PHF10 were performed as described previously. The transcription level was calculated by the 2^−ΔΔCt method. GAPDH was used as an internal reference gene. The primer sequences of LINC00242 were 5’-GCGGGAAGATTTCAGGCGCTTT-3’ (forward) and 5’-CAGGTGGTGAAGTGAGGAACAG-3’ (reverse). The primer sequences of PHF10 were 5’-GCACTCTAGGCTTAACAGCATT (forward) and 5’-AGCATGTTTGGCTGGATATTCTT-3’ (reverse). The primer sequences of GAPDH were 5’-GTCTCCTCTGACTTCAACAGCG-3’ (forward) and 5’-ACCACCCTGTTGCTGTAGCCAA-3’ (reverse).

The abundance of 22 immune cells in UCSC-PAAD and ICGC-PAAD was estimated via the CIBERSOFT website (https://cibersort.stanford.edu/index.php) (21). Monte Carlo sampling was used to calculate the p value. Only samples with CIBERSORT p < 0.05 were included in the further study. Based on the median value of PHF10 expression, we divided PAAD patients into a high-expression group (PHF10-H) and a low-expression group (PHF10-L). We used violin diagrams to show the difference in the content of 22 immune cells between PHF10-H and PHF10-L groups. From the previous studies, we obtained the immune-related gene set with 29 immune cell types and immune-related functions (22). Then, we used the ssGSEA algorithm to obtain the enrichment score with the “GSVA” package in R (23). PAAD samples from UCSC-PAAD and ICGC-PAAD were divided into three immune subgroups by the unsupervised clustering method based on ssGSEA scores. In order to further explore the difference of the immune microenvironment between three immune subgroups, the “estimate” package in R (24) was used to calculate the immune score, stromal score, and ESTIMATE score. Next, we use the “pheatmap” package in R to visualize the results. Additionally, the expression levels of PHF10 and the human leukocyte antigen (HLA) gene family were compared between three immune subgroups.

The Tumor Immune Estimation Resource (TIMER, https://cistrome.shinyapps.io/timer/) website (25) was used to evaluate the relationship between the expression level of PHF10 and immune cell infiltration. TIMER applies a statistical method called deconvolution to infer the abundance of TIICs according to RNA sequencing data. The relationship between immune subtype and mRNA level of PHF10 was explored via tumor–immune system interactions and DrugBank (TISIDB, http://cis.hku.hk/TISIDB/index.php) (26) database by the Subtype module.

All statistical analyses were operated via R software (version 4.0.3), GraphPad Prism (version 8.0), and IBM SPSS Statistics (version 25.0). The correlations between PHF10 and clinicopathological features were analyzed by the χ2 test. Univariate and multivariate analyses were used to determine prognostic factors. Spearman rank correlation analysis was used to estimate the correlation strength. Kaplan–Meier analysis was used to compare the survival rate among different groups. Unless noted otherwise, p < 0.05 was considered statistically significant.

After removing samples with missing values, a total of 177 patients were included in the study. The clinicopathological characteristics of the patients are shown in Supplementary Table 1 . We obtained 2,695 eRNA transcripts and 2,303 predicted target genes from previous studies. Next, the 2,695 eRNA transcript IDs were mapped to their corresponding 1,559 genes based on the Ensembl Gene ID from the Ensembl database (http://asia.ensembl.org/info/data/ftp/index.html). According to gene expression profile and clinical information from 177 PAAD patients, we identified 15 of the 1,559 candidate eRNAs, which were significantly correlated with patients’ OS (Kaplan–Meier, p < 0.01; Supplementary Table 2 ). Then, we found that all 15 eRNAs were significantly correlated with their target genes (R > 0.4, p < 0.01; Supplementary Table 2 ). The relationship between top 10 eRNAs associated with survival and its target gene is shown in Figure 1A .

LINC00242 is the most survival related eRNA, and it is significantly positively correlated with the predicted target PHF10 ( Figure 1A and Supplementary Table 2 ). Then, LINC00242 and its target PHF10 were selected for further analysis. According to the gene expression data of 179 PAAD tissues and 171 normal pancreatic tissues from the TCGA and GTEx databases, the mRNA level of LINC00242 was significantly upregulated in PAAD samples ( Figure 1B , p < 0.05). However, there was no significant difference in LINC00242 expression between pancreatic carcinoma and normal tissues in our validation set ( Supplementary Figure 3C ). The expression level of LINC00242 was significantly correlated with some clinical features of PAAD, including cancer status (p = 0.005), race (p = 0.02), and history of chronic pancreatitis (p = 0.004) ( Supplementary Table 1 ). Compared with G1, the expression level of LINC00242 of G2 was lower ( Figure 1C , p < 0.05). In addition, patients with advanced tumor stage were also correlated with low expression level of LINC00242 ( Figure 1D , p < 0.05). Furthermore, cancer status also associated with LINC00242 expression level. Tumor free status had higher gene expression level when compared with tumor status ( Figure 1E , p < 0.05). Interestingly, we found that patients of Asian descent had lower gene expression than patients of the other two races, including black or African American and white ( Figure 1F , p < 0.05).

Based on the median expression level of LINC00242, we divided 177 patients into high- and low-expression groups. High-expression LINC00242 was significantly correlated with better OS (Kaplan–Meier, p < 0.01, Figure 2A ). Unfortunately, we did not find a significant relationship between LINC00242 expression and prognosis as well as clinical parameters (TNM stage, grade, data not displayed) in ICGC-PAAD (N = 182), which may be related to the low abundance of non-coding RNA and the different sequencing methods ( Supplementary Figure 3A ). High PHF10 expression was also associated with better OS (p < 0.01, Figure 2B ). The analysis based on ICGC-PAAD showed that high expression of PHF10 was closely related to the better prognosis of the patients ( Supplementary Figure 3B ). What is more, the prognostic value of LINC00242 was explored in other cancer types from the TCGA database. LINC00242 also plays a prognostic role in KIRC and ACC. A high expression level of LINC00242 was significantly associated with poor OS in KIRC and ACC (p < 0.01, Figures 2C, D and Supplementary Table 3 ).

In addition, we found a significant co-expression correlation between LINC00242 and its target gene PHF10 in multiple cancer types, including PAAD, KIRC, ACC, TGCT, UVM, LIHC, and UCS (p < 0.01, Figures 3A–E and Supplementary Table 3 ). To further verify the regulatory relationship between LINC00242 and PHF10, we knocked down the expression of LINC00242 with transfection of LINC00242 siRNA. We used qRT-PCR to examine the efficiency of knockdown. As depicted in Figure 3F , the expression of LINC00242 was significantly downregulated after si-LINC00242 transfection (p < 0.05, Figure 3F ). Of note, compared with the si-control group, the expression of PHF10 was significantly downregulated in si-LINC00242 group (p < 0.05, Figure 3G ).

The univariate and multivariate Cox analyses were used to estimate the prognostic value of LINC00242 in PAAD. The univariate analysis showed that age (p < 0.05), grade (p < 0.05), PHF10 (p < 0.01), and LINC00242 (p < 0.01) were significantly associated with OS ( Figure 4A ). In the multivariate analysis, LINC00242 (p < 0.01) was still an independent prognostic factor of OS ( Figure 4B ).

Through the intersection of differentially expressed genes and co-expressed genes, we identified a total of 103 genes. GO and KEGG enrichment analyses showed that these genes were involved in many kinds of immune-related pathways and signaling pathways, such as T cell activation and lymphocyte differentiation ( Figures 4C, D ).

The HPA website was used to analyze the expression level of PHF10 in different tissues. The RNA expression of PHF10 in normal pancreatic tissue was relatively low ( Figure 5A ), but the protein level of PHF10 was relatively high among 43 normal human tissue types ( Figure 5B ). Additionally, the PHF10 mRNA level in PAAD samples was relatively low among 17 cancer types ( Figure 5C ); However, the expression of PHF10 was significantly upregulated in PAAD compared with normal tissues ( Figure 5D ). Based on GSE15471 analysis, the expression of PHF10 in pancreatic cancer tissues was significantly higher than that in normal tissues ( Supplementary Figure 3D ). In terms of the expression level of PHF10 protein, pancreatic cancer showed a low positive rate, ranking seventh from the bottom among 43 common cancer types ( Figure 5E ). The representative IHC images with different PHF10 expression levels, including low, medium, and high, are shown in Figure 5F . Based on web tool TIMER, we further examined PHF10 expression in multiple human cancers with RNA-seq data from the TCGA database. Interestingly, the PHF10 expression level varied significantly among different cancer types. PHF10 expression was significantly upregulated in CHOL, COAD, LIHC, and STAD and significantly downregulated in BLCA, BRCA, KICH, KIRC, KIRP, LUAD, THCA, and UCEC, than in normal tissues ( Supplementary Figure 1 ).

After performing the CIBERSORT algorithm, 134 tumor samples with p< 0.05 in the TCGA cohort were enrolled in this study. The landscape of immune infiltrations consisted by 22 immune cells is shown in Figure 6A . There was a significant positive correlation between activated CD4⁺ memory T cells and CD8⁺ T cells, and a significant negative correlation between CD8⁺ T cells and M0 macrophage ( Figure 6B ). To further investigate the effect of PHF10 on TIICs, all PAAD samples were assigned into high-expression (PHF10-H) and low-expression (PHF10-L) groups based on the median expression level of PHF10. Obviously, there was a significant difference in the proportion of TIICs between the two groups ( Figure 6C ). Compared with the PHF10-H group, the PHF10-L group contained a higher proportion of activated NK cells and M0 macrophages, but the proportion of naive B cells and resting CD4⁺ memory T cells was relatively lower (all p < 0.05, Figure 6C ).

We used the ssGSEA algorithm to obtain the enrichment score. Using an unsupervised clustering algorithm, PAAD samples were assigned into three immune infiltration clusters, including high immune cell infiltration cluster (Immunity-H, n = 97), medium immune cell infiltration cluster (Immunity-M, n = 28), and low immune cell infiltration cluster (Immunity-L, n = 53). Then, we calculated the Immune score, Stromal score, and ESTIMATE score based on the normalized gene expression data. The Immune score, Stromal score, and ESTIMATE score of the Immunity-H cluster are higher than the other two clusters ( Figure 7A ). The box plot also showed that with the increase in PHF10 expression, the expression levels of several HLA related genes were upregulated ( Figure 7B ). Additionally, the expression level of PHF10 in the Immunity-H cluster was significantly upregulated compared with the Immunity-M and Immunity-L clusters (all p< 0.05, Figure 7C ). Based on ssGESA analysis, pancreatic cancer patients in the ICGC-PAAD cohort can be divided into three categories, which was consistent with the results of the training set ( Supplementary Figure 4A ). In addition, there was infiltration of more B cell-naive and fewer M0 macrophages in the PHF10 high-expression group, which was consistent with the training set ( Supplementary Figure 4B ). At the same time, the expression of PHF10 in the Immunity-H group was significantly higher than that in the Immunity-M group ( Supplementary Figure 4C ).

Concomitantly, we further investigate the relationships between PHF10 expression and immune checkpoint genes (PD-L1 and CTLA-4) in PAAD via GEPIA. The results demonstrated that PHF10 was positively correlated with the expression of PD-L1 (R = 0.32, p < 0.01, Figure 7D ) and CTLA4 (R = 0.35, p < 0.01, Figure 7E ).

Then, we use the TISIDB database to explore whether PHF10 was correlated with the immune subtype of PAAD. All tumor samples in the TISIDB database were divided into six immune subtypes, including C1: wound healing, C2: IFN-gamma dominant, C3: inflammatory, C4: lymphocyte depleted, C5: immunologically quiet, and C6: TGF-b dominant). We found that the expression level of PHF10 varied significantly among the five immune subtypes (p <0.01, Figure 8A ). Next, we investigated the correlations between PHF10 expression and immune subtypes in human pan-cancer. Specifically, PHF10 expression was significantly correlated with immune subtypes in BLCA, STAD, LGG, BRCA, COAD, and KIRC (all p < 0.01, Figures 8B–G ). The relationship between PHF10 and immune infiltration in PAAD and human pan-cancer was explored via the TIMER database. PHF10 expression was significantly associated with the content of immune cells in PAAD, including CD4+ T cells (R = 0.47), CD8+ T cells (R = 0.52), B cells (R = 0.35), neutrophils (R = 0.46), macrophages (R = 0.54), and dendritic cells (R = 0.48) (all p < 0.001, Figures 9A–F ). Moreover, we further investigated the relationship between PHF10 and immune infiltration in human pan-cancer. The results showed that PHF10 mRNA level was associated with the immune infiltration in various types of cancer, including KICH, KIRC, SKCM, and THCA (all p < 0.05, R > 0.3, Supplementary Figures 2A-L ). Overall, these findings strongly suggest that PHF10 plays an important role in tumor immunity.