The research protocol was approved by the Ethics Review Committee of the Affiliated Hospital of Yunnan University, and all study participants gave written informed consent (approval number: 2022198). The study was carried out in accordance with the Declaration of Helsinki. All surgical procedures were performed under local anesthesia by the same surgeon. Thirty-four patients underwent elective pterygium surgery, including 5 males and 29 females. The average age was 56.8 years and the average duration of illness was 4.19 years (Supplementary Figure S1). During the surgical procedure, the affected pterygium tissue was removed alongside the surrounding excess loose conjunctiva. The integrity of the total RNA extracted from both pterygium and bulbar conjunctiva samples was evaluated using an Agilent 2,100 bioanalyzer. Subsequently, rRNA was eliminated from the total RNA to isolate the sample mRNA. This mRNA was then subjected to random fragmentation using divalent cations in NEB Fragmentation Buffer, followed by chain-specific fragmentation for mRNA construction. Initial library quantification was carried out with a Qubit2.0 Fluorometer, and the library was subsequently diluted to a concentration of 1.5 ng/ul. The size of the library inserts was determined using an Agilent 2,100 bioanalyzer, and QRT-PCR was utilized to accurately quantify the library’s effective concentration, which needed to exceed 2 nM. After assessing the quality of the genomic DNA, it was fragmented through mechanical interruption (ultrasound). The fragmented DNA then underwent purification, end-repair, 3′end adenylation, ligation to a sequencing adapter, and size-selection using agarose gel electrophoresis. The resulting polymerase chain reaction (PCR) product was amplified to generate the sequencing library. Sequencing was conducted on the Illumina NovaSeq 6,000 platform with a read length of 150 bp. Through quality control, trimming, deduplication, and alignment of the original fastq data by a high-throughput sequencing service provider, we obtained the gene expression matrix of transcriptome sequencing. The data presented in this publication has been archived in NCBI’s Sequence Read Archive (SRA) database (accession number: PRJNA1147595, URL: https://dataview.ncbi.nlm.nih.gov/object/PRJNA1147595?reviewer=3k3nnr66jke53qo77la1sbt91b). Supplementary datasets can be found in NCBI’s Gene Expression Omnibus (GEO) (accession numbers: GSE2513, URL: https://www.ncbi.nlm.nih.gov/geo/). GSE2513 includes four conjunctival samples and eight matched pairs of pterygium and control conjunctival samples (Wong et al., 2006). GSE51995 includes four conjunctival samples and four pairs of pterygium and control conjunctival samples (Hou et al., 2014).

To detect DEGs, we employed the “DESeq2” R package. The DESeq2 package was used for difference analysis of the original Counts matrix, and the analysis was carried out according to the standard process. The Variance Stabilizing Transformations (VST) method provided by package DESeq2 was used to Normalize the original Counts matrix. A Wilcoxon rank sum test was conducted to assess gene expression differences between the conjunctiva and pterygium samples. Genes were considered significantly differentially expressed if they met the criteria of adjusted p < 0.05 and |logFC| >1. For data visualization, we utilized the “pheatmap”, “ggpubr” and “ggplot2” R packages.

To analyze the DEGs, we employed the R tool known as “clusterProfiler.” (Yu et al., 2012) Using the “clusterProfiler” R package, we carried out both Gene Ontology (GO) (Ashburner et al., 2000) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses (Ogata et al., 1999). Gene set enrichment analysis (GSEA) was conducted on all genes (previously ranked based on their log2FC between analyzed groups) using the cluster profiler package. Enrichment was considered significant if the nominal false discovery rate (FDR) was <0.25 and the P-value was <0.05, referencing the ‘c2. cp.all.v2022.1. Hs.symbols.gmt’ gene set. Additionally, we conducted Gene Set Variation Analysis (GSVA) by utilizing the “GSEABase” and “GSVA” R packages (Hänzelmann et al., 2013) and referencing the ‘h.all.v2023.2. Hs.symbols.gmt’ gene set. These packages offer utilities for evaluating the enrichment or variability of gene sets within gene expression data, thereby enabling a thorough examination of pathway or gene set activities across different samples. P-values were calculated with the Benjamini–Hochberg method, and the terms with P-values <0.05 were considered to be significant.

For the construction of expression networks within the PRJNA1147595 dataset, we utilized the weighted gene co-expression network analysis (WGCNA) approach. This method was steered by adhering to the scale-free topology criterion, ensuring that the resultant network displayed a scale-free structure characterized by a power-law distribution of node connections. To ascertain the appropriate soft threshold power and establish adjacencies, we employed the “pickSoftThreshold” function from the WGCNA package. Once determined, the adjacency matrix was transformed into a topological overlap matrix (TOM). Hierarchical clustering analysis was then conducted based on the dissimilarity derived from the TOM.

To delineate co-expressed gene modules, we adopted the dynamic tree cutting technique, setting a minimum module size of 50. This approach allowed us to categorize genes into distinct modules according to their expression patterns. Subsequently, we evaluated the relationship between these gene modules and the trait of interest, pterygium, by considering both gene significance (GS) and module membership (MM) values. GS reflects the correlation between gene expression and pterygium, whereas MM measures the extent of co-expression within a given module. By integrating GS and MM values, we pinpointed the key modules linked to pterygium, suggesting their potential functional importance in the context of the disease.

The gene sets of the significant difference modules analyzed by WCGNA were intersected with the expression difference genes obtained by DESeq2 analysis. To identify the key hub genes, two distinct machine learning algorithms were employed: support vector machine (SVM) and random forest (RF). The SVM algorithm was executed using the “e1071”package, offering a comprehensive toolkit for SVM model training and classification. The RF algorithm was applied using the “randomForest” package, where the error rate was computed across 1 to 500 trees. The optimal tree count was determined by selecting the value that minimized the error rate while ensuring stability. Subsequently, the RF classifier was leveraged to compute feature importance scores. This selection utilized the Gini coefficient method to pinpoint the most pertinent genes. The intersection of the top five feature genes from both the SVM and RF models was used to pinpoint crucial hub genes involved in the disease process.

Utilizing the R programming language, receiver operating characteristic (ROC) curves were generated for the hub genes. To evaluate the distinctive capabilities between pterygium tissues and conjunctiva tissues, the area under the curve (AUC) was calculated for the respective ROC curves of these hub genes.

According to the CIBERSORT deconvolution algorithm (Newman et al., 2015), we analyzed the infiltration of 22 kinds of immune cells in pterygium and conjunctival tissues, and compared whether there was a significant difference between the two groups. To assess the degrees of immune cell infiltration in both conjunctiva and pterygium samples, we utilized the CIBERSORT technique. This method enabled us to measure gene expression levels quantitatively and ascertain the infiltration levels among patients. The immune cell populations encompassed both immune-stimulating and immune-inhibiting cells. To examine the relationships between these immune cells, we conducted a correlation analysis using the Spearman coefficient and depicted the results via correlation heatmaps. At the same time, we also analyzed the correlation between hub gene and immune cells. We considered findings statistically significant if they met the threshold of p < 0.05.

Patient Sample Preparation. There were 4 primary pterygium specimens obtained in Affiliated Hospital of Yunnan University during pterygium surgery. The control conjunctival tissues were derived from healthy conjunctival tissues on the temporal side of the surgical eye of the same patient, with a size of 1.5 mm × 1.5 mm.

Tissue RNA Extraction and RT-qPCR. All specimens were collected during the operation and immediately placed in RNA protective agent, stored at −20°C for usage. RNA was extracted from 80 mg specimen tissue using TRIZOL reagent. The relative expression of hub genes was measured using a reverse kit (HiScript® III RT SuperMix for qPCR (+gDNA wiper)) and ChamQ™ Universal SYBR® qPCR Master Mix (Vazyme company, Nanjing, China) on CFX96 fluorescence qPCR instrument. And the PCR reaction was carried as follows: predenaturation at 95°C for 30 s, denaturation at 95°C for 3–10 s, and annealing extension at 60°C for 10–30 s, for 40 cycles. Using GAPDH as internal references, the relative expression of genes was calculated by the 2-△Ct method.

The identification of key miRNAs targeting the two characteristic genes, as well as the screening of potential drug candidates, was conducted through a rigorous and comprehensive process. We utilized the Targetscan database (version 7.2, accessed via the Enrichr platform at http://amp.pharm.mssm.edu/Enrichr/) to forecast miRNAs that could potentially target the selected genes (Agarwal et al., 2015). This involved inputting the gene symbols into the Targetscan interface and retrieving a list of predicted miRNAs based on the presence of conserved target sites within the 3′UTRs of the genes. The predicted miRNAs were then filtered based on their conservation across species, expression levels in relevant tissues, and previous reports of their involvement in related biological processes. Subsequently, to identify potential drug candidates that could modulate the expression or activity of these genes, we accessed the DSigDB database through the Enrichr platform. This database integrates drug signatures from various sources, allowing us to search for drugs that have been reported to affect the expression or activity of the selected genes.

The statistical analyses were performed using R software (version 4.2.1). All p values were two-sided, and significance was indicated by p < 0.05.

Based on the RNA-Seq data, a total of 718 DEGs were obtained, and these genes reached the threshold adjust p-value of <0.05 (Supplementary Table S1). Among these DEGs, 254 transcripts were upregulated in the disease group against the health one, while 464 transcripts were found to decrease in pterygium tissues (Figure 2A). In addition, the relative expression levels are shown in a heatmap plot and the top 50 upregulated and downregulated genes were shown using a volcano plot (Figure 2B).

GO and KEGG pathway enrichment analyses were conducted to assess the function of these 718 DEGs. The clusterProfiler package was used for GO and KEGG pathway enrichment analysis to find potential biological pathways. The top 5 BP (biological process), CC (cellular component) and MF (molecular function) terms with the lowest p-values in each category were selected and visualized using bubble plots (Figure 2C). The DEGs were primarily enriched in BP terms such as “leukocyte migration (GO:0050900)”, “leukocyte chemotaxis (GO:0030595)”, “cell chemotaxis (GO:0060326)”, “leukocyte mediated immunity (GO:0002443)” and “regulation of inflammatory response (GO:0050727)”. Genes were also enriched in CC terms like “external side of plasma membrane (GO:0009897)”, “tertiary granule (GO:0070820)”, “secretory granule membrane (GO:0030667)”, “ficolin-1-rich granule (GO:0101002)” and “cornified envelope (GO:0001533)”. Additionally, the DEGs were enriched in MF terms including “cytokine activity (GO:0005125)”, “glycosaminoglycan binding (GO:0005539)”, “immune receptor activity (GO:0140375)”, “chemokine activity (GO:0008009)” and “chemokine receptor binding (GO:0042379)”. In the KEGG pathway enrichment results, many pathways related to inflammation and signal transduction were enriched, such as IL-17 signaling pathway (hsa04657), chemokine signaling pathway (hsa04062), TNF signaling pathway (hsa04668), MAPK signaling pathway (hsa04010) and cytokine-cytokine receptor interaction (hsa04060) (Figure 2D) and so on (Supplementary Table S2).

Our reference gene set was ‘c2. cp.all.v2022.1. Hs.symbols.gmt’. The datasets underwent GSEA enrichment analysis to identify significant enrichment based on the criteria of FDR <0.25 and p < 0.05. In gene set enrichment analysis (GSEA), the significant enrichment in upregulated pathways including REACTOME_MET_ACTIVATES_PTK2_SIGNALING, REACTOME_COLLAGEN_FORMATION, REACTOME_MET_PROMOTES_CELL_MOTILITY, KEGG_OXIDATIVE_PHOSPHORYLATION and so on (Figures 3A–D). The significant enrichment in downregulated pathways including WP_MAPK_SIGNALING_PATHWAY (Figure 3E) and so on (Supplementary Table S3). GSVA enrichment analysis was conducted on the PRJNA1147595 dataset, revealing distinct pathways (Supplementary Table S4). By utilizing the gene set variation analysis (GSVA) package and referencing the “h.all.v2023.2. Hs.symbols.gmt” gene set (Hänzelmann et al., 2013), the gene expression matrix data were subjected to GSVA. The GSVA variance analysis was performed using the limma package of R software. Differential pathways were filtered based on an p-value <0.05 and |log2FC| > 0.25 (Supplementary Table S4). The differential pathways in the PRJNA1147595 dataset encompassed HALLMARK_TNFA_SIGNALING_VIA_NFKB, HALLMARK_INFLAMMATORY_RESPONSE, HALLMARK_OXIDATIVE_PHOSPHORYLATION, HALLMARK_IL6_JAK_STAT3_SIGNALING and HALLMARK_ALLOGRAFT_REJECTION (Figure 3F).

Initially, expression data from 34 pterygium samples and 34 conjunctiva samples were utilized to establish co-expression modules through the WGCNA algorithm. We focused on the top 25% of genes exhibiting the greatest variability for subsequent investigation. Utilizing the “flashClust” package, we conducted cluster analysis on these samples, with the findings illustrated in Figure 4A. At a power value of 11, an independence level of 0.9 was attained, accompanied by an elevated mean connectivity (Figures 4B, C). By employing the dynamic cutting technique, we discerned seven unique gene co-expression modules within pterygium, and a heatmap of the TOM was also generated (Figures 4D, E). The genes within these seven modules were concurrently employed to explore the relationship between module eigengenes and clinical characteristics (Figure 4F). Notably, 3,956 genes were divided into 7 modules. Among them, salmon module (cor = −0.66, P = 1e-9) and turquoise module (cor = 0.34, p = 0.005) were significantly correlated with pterygium phenotype (Figure 4F).

A Venn diagram was constructed to visually represent the relationships among the differentially expressed genes (DEGs) and the salmon and turquoise modules identified through Weighted Gene Co-expression Network Analysis (Figure5A). Subsequently, two distinct machine learning algorithms—SVM and RF—were leveraged to meticulously screen and pinpoint key signature genes within the pterygium dataset. For the RF classifier, an optimal number of 500 trees was selected, a decision grounded in an evaluation of error rates and classifier stability across varying tree counts. By analyzing the Mean Decrease Gini results, we identified top five pivotal hub genes: KRT10, SPRR1B, SERPINB13, PRSS3and NGEF (illustrated in Figure 5B). The SVM algorithm was also employed to rigorously select hub genes associated with pterygium, revealing five significant genes: KRT10, KRT6B, BTG2, ASPG and NGEF (Figure 5C).

To identify crucial hub genes involved in pterygium, we examined the intersection of the top five feature genes from both the SVM and RF models (Figure 5D). The model’s performance was evaluated using the area under the ROC curve (AUC) (Figures 5E, F). Impressively, the dataset yielded an AUC of 0.921 (KRT10) and 0.904 (NGEF), indicating a high level of accuracy in classifying gene expression data. Notably, KRT10 and NGEF were consistently identified as hub genes associated with pterygium across all two machine learning methods and WGCNA.

We downloaded the original data of GSE2513 and GSE51995 from GEO website, and conducted probe annotation and data standardization on the original expression data to obtain the gene expression matrix. The datasets GSE2513 and GSE51995 were merged utilizing the R packages dplyr and purrr, and the shared genes were identified based on the intersection of genes derived from sequencing analysis. Subsequently, the cleaned data were standardized using the stats, preprocessCore, and limma R packages. The resultant standardized data were then subjected to further analysis (Smyth and Gentleman, 2013). KRT10 achieved an AUC of 0.906 in GSE2513 and GSE51995 combined data and NGEF achieved an AUC of 0.625 (Supplementary Table S2).

The CIBERSORT deconvolution algorithm was employed to quantify the infiltration of 22 distinct immune cell types within pterygium and bulbar conjunctiva tissues (Figure 6A). Our analysis revealed significant differences in the levels of five immune cells between these two tissue types. Specifically, M1 macrophages, resting dendritic cells, and activated dendritic cells were notably elevated in the pterygium group, whereas resting mast cells and M2 macrophages were significantly increased in the bulbar conjunctiva group (Figure 6A). Furthermore, a general correlation was observed among the 22 immune cell types within pterygium samples (Figure 6B). Correlation analysis between these immune cells and KRT10 (a hub gene) in pterygium tissues indicated that dendritic resting cells (cor = 0.58, p = 1.63e-07) and M1 macrophages (cor = 0.43, p = 3.21e-05) were positively correlated with KRT10 expression, while Neutrophils also exhibited a significant negative correlation with KRT10 (cor = −0.38, p = 0.0016) (Figure 6C). Correlation analysis between these immune cells and NGEF in pterygium tissues indicated that Macrophages M1 (cor = 0.61, P = 2.57e-08), Dendritic cells resting (cor = 0.54, P = 1.95e-06), Macrophages M0 (cor = 0.42, P = 0.0003), T cells regulatory (cor = 0.42, P = 0.0004) and Mast cells resting (cor = 0.36, P = 0.0029) were positively correlated with NGEF expression, while Dendritic cells activated (cor = −0.35, p = 0.0034) and Mast cells activated (cor = −0.33, p = 0.0063)exhibited a significant negative correlation with NGEF (Figure 6D).

To validate our findings, we employed the combined dataset of GSE2513 and GSE51995 from the GEO database as a verification cohort. Following CIBERSORT analysis on this integrated dataset, we observed that there were no statistically significant differences in the proportions of M1 macrophages, resting dendritic cells, activated dendritic cells, resting mast cells, and M2 macrophages between the pterygium and conjunctiva groups. Additionally, a comprehensive correlation analysis revealed notable associations among the 22 immune cell types within the pterygium samples (Supplementary Table S3).

In our analysis of RNA-Seq data derived from 68 samples, we observed a statistically significant elevation in the expression levels of KRT10 and NGEF genes in the pterygium group compared to the bulbar conjunctiva group, with a P-value <0.05 (as depicted in Figures 7A, B). To validate these findings, we conducted a qPCR experiment, which confirmed a significantly higher expression of KRT10 and NGEF in the pterygium group relative to the bulbar conjunctiva group, with a P-value of 0.04 and 0.02 as illustrated in Figures 7C, D and primer sequences are shown in Supplementary Table S5.

To pinpoint crucial miRNAs and potential drug candidates that target the two identified feature genes, we conducted a comprehensive analysis using data sourced from the Targetscan and DSigDB databases. This rigorous screening process, employing a statistical significance threshold of p < 0.05, led to the identification of three pivotal miRNAs (depicted in Figure 8A). Our findings revealed intricate interaction networks: KRT10 was found to interact with one distinct miRNA (hsa-miR-4750), while NGEF was associated with hsa-miR-450a and hsa-miR-4737.

Furthermore, we undertook an extensive screening to identify the drug molecules, guided by P-value < 0.05 as per the DSigDB database (illustrated in Figures 8B, C). We conducted a rigorous screening process and identified KRT10 and NGEF as crucial genes implicated in the pathogenesis of pterygium. Subsequently, we presented only the top four drugs that significantly influenced KRT10 expression, along with benzo [a] pyrene CTD 00005488, which demonstrated the ability to jointly regulate both KRT10 and NGEF genes. In the case of the NGEF gene, our analysis revealed only four drugs with statistical significance (p < 0.05), including benzo [a] pyrene CTD 00005488, and we have accordingly reported these four drugs. Among these, benzo [a]pyrene (CTD 00005488) emerged as a prominent candidate, demonstrating interactions with all two feature genes and achieving a notable combined score of 94002.9. The remaining drug candidates exhibited interactions with KRT10 and NGEF, offering valuable insights for the advancement of pterygium treatment research and development.