Germline alterations were filtered out from the tumor samples using the Genome Analysis Toolkit (GATK), focusing on retaining somatic alterations. WES revealed that our local patients’ dominant variant classification, variant type, and single nucleotide variant (SNV) class were consistent with the profiles observed in the TCGA-THCA (Thyroid Carcinoma, THCA)cohorts ( Figure S2 ). However, there were notable differences in the top mutated genes. While canonical BRAF mutations were observed in over 60% of our patients, other top mutated genes in our cases included OR51M1, MAGEB16, EBLN2, ZNF714, SGIP1, PCSK9, NPAP1, KRTAP5-7, and A1BG ( Figures 1A , S2A ), instead of NRAS and HRAS observed in TCGA-THCA ( Figure S2B ). The frequency of different types of mutations was analyzed, revealing a higher frequency of transitions compared to transversions ( Figure 1B ).

Furthermore, genomic data indicated a higher tumor mutation burden (TMB) in our local patients compared to the TCGA-THCA dataset ( Figure 1C ). Subsequently, we examined the driver mutations in our local patients. BRAF (Z-score=5.546, p<0.001) and PNPLA5 (Z-score=5.546, p<0.001), a gene closely associated with lipid metabolism (17), were the only two genes implicated in tumor pathogenesis ( Table 1 ). Regarding pathways related to tumor activities, only the RTK-RAS pathway was significantly affected due to BRAF mutations ( Figure 1D ). However, it is essential to note that our sequencing technique focused on the exonic regions and may have provided limited insights into gene fusions and chromosomal aberrations occurring in intronic or other non-exonic regions. Further investigation using complementary techniques, such as whole genome sequencing or Long-read-sequencing analysis, may be necessary to better understand these genomic alterations in PTC.

A total of 25,978 non-low expression genes’ transcriptomic profile was analyzed using an R script after filtering out genes with transcript per million mapped (RPKM) values less than 3. Principal component analysis (PCA), using the FactoMineR R package, revealed a clear distinction in gene expression between normal and tumor samples ( Figure 2A ). Except for one paired sample (EA-002), the normal samples exhibited a convergence compared to the tumor samples, while tumor samples showed relative heterogeneity. The log2 fold changes and corresponding false discovery rates (FDR) to identify differentially expressed genes (DEGs), employing criteria of |log2 (fold change) | > 2 between normal and tumor samples and FDR < 0.001. 572 unique genes were identified, using the edgeR R package, as DEGs, consisting of 353 upregulated genes and 219 downregulated genes, as illustrated in the volcano plot ( Figure 2B ).

To validate the ability of the identified DEGs to accurately reflect the gene expression differences between normal tissue and tumors, we selected samples from the TCGA-THCA dataset. Specifically, based on their clinical information, we included 59 normal samples, 339 papillary adenocarcinoma-NOS (not otherwise specified), 101 papillary carcinoma-follicular variant, and 35 papillary carcinoma-columnar cell samples ( Figure 2C ). The heatmap demonstrated that the DEGs we identified successfully distinguished normal samples from most tumor samples ( Figure 2C ). However, it is noteworthy that the transcriptomic profiles of normal samples and follicular variants of papillary thyroid carcinoma (FVPTC) appeared similar in the heatmap. This observation may be attributed to the low risk of adverse outcomes associated with FVPTC (20).

We performed Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses to explore the potential mechanisms underlying PTC ( Figures 2D, E ). The GO analysis identified the top eight significant pathways in biological processes, cellular components, and molecular functions. Notably, these pathways were predominantly related to the extracellular matrix (ECM), including extracellular matrix organization, collagen-containing extracellular matrix, and extracellular matrix structural constituent ( Figure 2D ). Previous studies have highlighted the dysregulation of ECM in PTC, which is closely associated with tumor activities such as migration and invasion (21). In addition, the KEGG analysis revealed seven significant pathways, including protein digestion and absorption, cell adhesion molecules, ECM-receptor interaction, malaria, AGE-RAGE signaling pathway in diabetic complications, focal adhesion, and axon guidance. Of particular interest is the interaction between ECM and cell adhesion ( Figure 2E ), the latter of which has been reported in several studies to play a crucial role in tumor growth and invasion (22, 23). Furthermore, it is intriguing that the GO and KEGG analyses unveiled some neuronal system-related pathways, potentially indicating the presence of neuronal cell-adhesion molecules in PTC, such as NrCAM (24). The gene-concept network helps illustrate how specific genes are linked to various pathways and concepts ( Figures 2F, G ), highlighting their potential roles in PTC development and progression.

To identify potential hub genes for further validation, a loose screening criterion was applied for DEGs, selecting genes with an FDR < 0.05 in the differential analysis. A total of 4,828 genes from the 6 tumor samples were then used to construct a co-expression network. The analysis was performed using the R package WGCNA, where a soft-threshold power of 16 was chosen to ensure appropriate scale independence and mean connectivity of the network. Based on the resulting cluster dendrogram, genes with similar co-expression patterns were grouped into 20 modules ( Figure 3A ). Eigengenes were correlated with the clinical data to assess the correlation between modules and clinical characteristics ( Table 1 ). Notably, the MEred module exhibited a high correlation with tumor size (r=0.89, p=0.02), while MEgrey60 (r=-0.81, p=0.05) and MElightcyan (r=-0.84, p=0.03) showed negative correlations ( Figure 3B ). Scatter plots were generated to depict the correlations between module membership (MM) and gene significance (GS) for the MEred (r=0.66, p<0.001), MEgrey60 (r=0.43, p<0.001), and MElightcyan (r=0.23, p=0.05) modules ( Figures 3C–E ), confirming the significant correlation of these modules with tumor size. Finally, we exported the edge data of genes within these three modules and utilized Cytoscape v3.10.0 for further analysis and visualization. This enabled a more detailed examination of the interactions and relationships among the genes within these modules.

DDA was utilized for quantitative proteomic analysis. Before analysis, proteins with undetectable abundance or those without annotations in the NCBI Reference Sequence Database were excluded from the dataset. A total of 3,624 proteins were identified for further analysis.

The PCA plot of the proteomic data demonstrated a distinct separation between normal and tumor samples, except sample EA-001_C, which appeared to be closer to the normal sample cluster ( Figure 4A ). This heterogeneity among the tumor samples was consistent with the findings in the transcriptomic PCA plot, indicating potential molecular variations within the PTC samples at the protein level.

To identify differentially expressed proteins (DEPs), we employed a similar strategy as used for DEGs. DEPs were defined using a threshold of |log2 (fold change) | > 2 and an acceptable FDR cutoff of < 0.05. Using the limma R package, we identified 79 DEPs, with 35 proteins upregulated and 44 proteins downregulated, as shown in the volcano plot ( Figure 4B ). Similar to the transcriptomic analysis, we extracted gene expression data corresponding to the DEPs from the TCGA-THCA dataset. As anticipated, the DEPs also demonstrated the capability to distinguish tumor samples from normal samples, although not as efficiently as the DEGs ( Figure 4C ). This observation could be attributed to the fact that direct protein abundance data were not utilized, which may have affected the discriminatory power of the DEPs. GO and KEGG enrichment analyses were conducted to investigate the functional implications of the DEPs further. Besides enriching ECM-related pathways, such as collagen-containing extracellular matrix, we observed frequent enrichment of detoxification-associated pathways, like cellular oxidant detoxification and lipid metabolism, including fatty acid binding and long-chain fatty acid binding ( Figure 4D ). KEGG analysis revealed five pathways in which the DEPs may participate: ECM-receptor interaction, tyrosine metabolism, hematopoietic cell lineage, malaria, and cholesterol metabolism ( Figure 4E ). These findings were partly consistent with the results from the GO analysis, suggesting potential roles of ECM, lipid metabolism, and detoxification processes in PTC. The gene-concept networks were also exhibited to visualize the link between proteins and specific pathways ( Figures 4F, G ).

PPI analysis was performed by using the STRING database to explore protein-protein interactions (PPI) and identify potential hub genes. The resulting PPI networks were visualized ( Figure 5A ), and the edge file was downloaded and imported into Cytoscape. Similarly, we imported the edge file obtained from the previous WGCNA into Cytoscape to search for hub genes within the network.

To identify the hub genes, we utilized the cytoHubba app in Cytoscape and ranked the nodes based on maximal clique centrality (MCC), a measure of network connectivity. From the PPI analysis, we selected AEBP1, CD34, COL12A1, GYPA, ITGA2B, KRT19, POSTN, and THBS2 as hub genes ( Figure 5B ). From the WGCNA analysis, we selected ABR, C2CD2L, CDC42EP5, CTSA, GPX1, LASP1, and SFN as hub genes ( Figure 5C ). However, ENSG00000234943, a lncRNA, was excluded from further analysis.

We examined the overlap among somatic alterations, DEGs, and DEPs to integrate the results from multiple omics analyses. Through this analysis, we identified two additional hub genes, TPO and AHNAK2, which had somatic alterations and were also differentially expressed at the transcriptional and protein levels ( Figure 5D ). Therefore, we obtained a set of 17 hub genes through the multi-omics analysis. These hub genes represent key network elements and merit further investigation to elucidate their functions in predicting the prognosis and diagnosis of PTC. Among the identified hub genes, AEBP1, CD34, COL12A1, ITGA2B, THBS2, and TPO exhibit lower expression levels in thyroid tumors. On the other hand, AHNAK2, GPX1, KRT19, and SFN show higher expression levels in thyroid tumors ( Figure S3 ). These differential expression patterns suggest their potential as diagnostic markers for distinguishing thyroid tumors from normal tissue.

Based on the Kaplan-Meier plotter employing TCGA-THCA (n=502), the associations between the expression of hub genes and patient survival were found ( Figures 6 , S4 , S5 ):

Overall Survival (OS) ( Figures 6A–H , S4 ):

Recurrence-Free Survival (RFS) ( Figures 6I–P , S5 ):

These findings suggest that the expression levels of these hub genes may have predictive value in determining patient survival outcomes in papillary thyroid carcinoma.

RNA-seq data from 534 samples (475 PTC samples and 59 normal samples) obtained from the TCGA-THCA cohorts were analyzed to evaluate the diagnostic potential of the hub genes in PTC. First, we examined whether the expression levels of individual genes could distinguish between tumor and normal tissue. Using receiver operating characteristic curve (ROC) analysis, we assessed the performance of each gene in differentiating PTC from normal tissue to select genes qualified to be employed to construct a clinical predictive model ( Figures 7 , S6 ). According to the events per variable (EPV) rule, we considered the area under the curve (AUC) values above 0.85, suggesting these genes show promising diagnostic efficacy in controlling the number of variables. Five genes, ABR (AUC=0.945), AHNAK2 (AUC=0.877), CTSA (AUC=0.901), GPX1 (AUC=0.914), TPO (AUC=0.894), performed significantly in the ROC analysis ( Figure 7 ). The high AUC values obtained from the ROC analysis support the potential utility of these genes in constructing a robust diagnostic model for PTC.

To develop a diagnostic model for PTC, logistic regression analysis was performed using the glmnet R package. The gene expression data from the TCGA-THCA cohort was preprocessed by assigning scores of 1-10 based on percentile ranking within each gene’s expression data, with higher expression levels receiving lower scores. Tumor samples were given a score of 1, while normal samples were given 0. The logistic regression model was built using the expression levels of the following genes: TPO, GPX1, ABR, AHNAK2, and CTSA. The glmnet package with the “family=‘binomial’” parameter was used for logistic regression analysis. The odds ratio (OR) with its 95% confidence interval (95% CI) and p-values were calculated and presented in a forest plot ( Figure 8A ). Stepwise regression was performed using the step function. Among the genes, ABR, AHNAK2, GPX1, and TPO significantly associated with PTC diagnosis (p<0.05) and were identified as pivotal elements in the clinical model. Then, we achieved the formula of “y=0.4211*TPO-0.5256*GPX1-0.7043*ABR-0.4178*AHNAK2 + 13.8133”. A nomogram was constructed using the rms R package, incorporating these four genes as a reference for diagnosing PTC ( Figure 8B ). The calibration curve demonstrated good prediction efficacy of the model ( Figure 8C ). The diagnostic performance of the model was evaluated using ROC analysis, yielding an AUC value of 0.969 ( Figure 8D ), which outperformed the individual genes’ diagnostic potential ( Figures 7 , 8E ). External validation was conducted using the GSE33630 dataset (49 PTC and 45 normal samples), yielding an AUC value of 0.959 ( Figure 8F ), indicating the model’s capability to diagnose PTC.

Pan-cancer expression data of the four genes were exhibited for a further selection of other tumors to assess the model’s specificity for PTC ( Figure S7 ). The model was applied to different types of tumors from the GSE63514 (cervical squamous cell carcinoma), GSE132305 (cholangiocarcinoma tumor), GSE53757 (kidney renal clear cell carcinoma), GSE121248 (liver hepatocellular carcinoma), and GSE43458 (lung adenocarcinoma) datasets. The model exhibited inferior diagnostic potential in these tumor types compared to PTC, confirming its specificity for PTC diagnosis ( Figure S8 ). These results suggest that the diagnostic model based on the expression levels of ABR, AHNAK2, GPX1, and TPO has strong diagnostic efficacy for PTC and demonstrates specificity for PTC compared to other tumors.