The study flowchart is shown in Figure 1. Three gene expression datasets (GSE27155, GSE60542 and GSE33630) were collected from the GEO database. GSE27155 and GSE60542 were integrated into the training set, and GSE33630 was selected as the validation set. Differentially expressed genes (DEGs) were defined by the R package “limma.” Gene Ontology (GO) function enrichment analyses and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses were conducted through the “clusterProfiler” package of R based on the DEGs in the training set. Also, immune cell infiltration analysis was performed by R. We determined 10 important genes through random forest from the DEGs by the “randomForest” package. Further, an artificial neural network diagnostic model was established through the “neuralnet” package using the 10 important genes and was assessed by AUC. Finally, the validity of the ANN diagnostic model was validated with the performance of a validation set.

The GSE27155, GSE60542, and GSE33630 datasets were collected via the Gene Expression Omnibus database (GEO; https://www.ncbi.nlm.nih.gov/geo/). Then, gene names were obtained by transforming probe names by R software. If several probes could be matched to identical gene, the expression data with gene was replaced by its mean expression value. Finally, 55 samples (4 normal samples and 51 PTC samples) in GSE27155, 63 samples (30 normal samples and 33 PTC samples) in GSE60542 and 94 samples (45 normal samples and 49 PTC samples) in GSE33630 were utilized in this research (Table 1). GSE27155 and GSE60542 were combined into a training set, and GSE33630 as a validation set.

Differentially expressed genes (DEGs) were identified in 34 normal samples and 84 PTC samples of the training set through the R package “limma” using the classic Bayesian data analysis. |log2FC| >1.5 and adjusted p Value <0.05 were set as threshold. Then, we obtained 94 DEGs involving 53 up-regulated genes and 41 down-regulated genes. R package “pheatmap” and “ggplot2” were utilized to conduct the heatmap and volcano plot.

To explore the biological significance of these DEGs, Gene Ontology (GO) enrichment analysis (adjusted p value < 0.05) categorizing genes into biological process (BP), cellular component (CC), and molecular function (MF) was performed. Meantime, Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis was used to describe metabolic pathways (p value < 0.05). Enrichment analysis of DEGs was conducted through the R package “clusterProfiler.”

As a crucial component of the tumor environment, immune cells impacted the development and prognosis of tumors, their composition as well as function in various types of tumors were different. On the one hand, some immune cells have the roles to be favorable targets for immunotherapy. On the other hand, some may react negatively, even resulting in drug resistance. Given these causes, figuring out the ingredient and possible influence of immune cells in PTC was beneficial to identify valuable therapeutic targets. We downloaded a reference including human gene expression in immune cells. Then, we exploited the reference and a file including the gene expression of each sample to complete immune cell infiltration analysis and obtain the result involving the expression of immune cells of each sample by R software. Based on the result, we performed correlation analysis of immune cell by the “corrplot” package and difference analysis of immune cells by the “vioplot” package.

The important genes were identified by a random forest classifier via “randomForest” package of R. Firstly, the parameter “ntree” was set as 500 to find the best number of trees. By calculating the error of cross validation, 15 was selected as the best number of trees represented the minimum error of cross validation. Then, parameter 15 was used to reconduct random forest and the importance score of genes was obtained. Genes with an importance score greater than 2 were seen as PTC significantly related genes. The R package “heatmap” was utilized to draw a heatmap based on important genes.

Based on 10 important genes identified by RF, an artificial neural network model was developed via the R package “neuralnet.” At first, the expression of 10 important genes was converted into “gene tag” based on their expression levels. In the case of a certain sample, the expression level of a specific gene was compared to the median of all sample expression level. Among the up-regulated genes, if the expression level is higher than the median, it will be valued as 1, otherwise 0. Among the down-regulated gene, if expression level is higher than median, it will be valued as 0, otherwise 1. Then, we finished a “gene tag” sheet. Next, the hidden layers of ANN were set as 5 to obtain a gene weight calculated by “gene tag”. Finally, the ANN diagnostic model was established. We assessed the model in the training set. Also, we validated the model in the validation set, and its diagnostic performance was evaluated by AUC.

In total, 53 significantly up-regulated genes and 41 significantly down-regulated genes were determined between normal samples and PTC samples based on |log2FC| >1.5 and adjusted p value <0.05 as threshold in training sets GSE27155 and GSE60542. Heatmaps (Figure 2A) and volcano plots (Figure 2B) of DEGs showed favorable discrimination of gene expression.

To reveal the biological importance of DEGs in the mechanism of PTC, we conducted GO and KEGG enrichment analysis using the package “clusterProfiler” of R software. Amid GO enrichment analysis results (adjusted p value cutoff = 0.05), concerning BP included wound healing, regulation of body fluid levels, blood coagulation and hemostasis; CC involved collagen-containing extracellular matrix; MF contained extracellular matrix structural constituents and other important functions (Figure 3A). KEGG pathway enrichment analysis (p value cutoff = 0.05) indicated that the DEGs were significantly related to tyrosine metabolism, complement and coagulation cascades, ECM–receptor interaction, and cell adhesion molecules (Figure 3B).

We utilized a reference including human gene expression in immune cells and a file involving gene expression of each sample to perform immune cell infiltration analysis. Then, we conducted difference analysis of immune cells by the “vioplot” package and correlation analysis of immune cells by the “corrplot” package. Compared with normal samples, T cells gamma delta, macrophages M2, dendritic cells resting, dendritic cells activated and mast cells resting were higher in PTC samples (p < 0.05). Conversely, macrophages M1, mast cells activated, and eosinophils were lower in PTC(p < 0.05) (Figure 4A). And macrophages M1 were correlated with dendritic cells activated and T cells gamma delta (correlation, 0.46 and 0.47) (Figure 4B).

We performed a random forest based on 94 DEGs to determine important genes. Considering the relationship plot between error of cross validation and the number of decision trees, we chose 15 trees as the parameter of the final model, which revealed a minimum error of cross validation in the model (Figure 5A). Eventually, 10 important genes were determined on the condition that their importance score was greater than 2 for the following analysis (Figure 5B). Amid the 10 important genes, GALNT7 was the most important. With heatmap, 10 important genes could be differed from normal samples and PTC samples. Among them, CITED1, AMIGO2, PSD3, GALNT7, and PROS1 were down-regulated in normal samples and up-regulated in PTC samples. TFF3, SLC4A4, AOX1, IPCEF1, and TFCP2L1 were up-regulated in normal samples and down-regulated in PTC samples (Figure 5C).

Expression of important genes determined by RF was transformed into “gene tag” marked as 0/1. The weight of all genes was calculated for optimal discrimination between normal samples and PTC samples. Then, an ANN diagnostic model based on gene weight was established (Figure 6). Performance of the model had an AUC of 0.988 (Figure 7A) in the training set and 0.968 (Figure 7B) in the validation set, indicating that the model was very satisfactory in diagnosing PTC. The results demonstrated that we had developed a precise diagnostic model between PTC and normal samples.