A variety of microarray studies based on thyroid tissue were collected by searching the key word “thyroid nodules” in the Gene Expression Omnibus (GEO) database (Barrett et al., 2013). These collected datasets should meet the following criteria (Yang et al., 2020b): 1) the gene expression profiling was conducted using cDNA microarray for “Homo Sapiens”; 2) the tissues analyzed were thyroid nodules; 3) raw data could be available for further analysis; and 4) the collected datasets should consist of one group of malignant samples and another group of benign ones. As a result, five independent transcriptomic datasets were collected, and each comprised both benign and malignant thyroid nodules. The detailed information of these five collected datasets is provided in Table 1, including dataset ID, number of samples, microarray platform, and tissue indicated in the original publication and references.

To enhance the consistency and classification capacity, all datasets in this study (Table 1) were combined to discover the key genes of thyroid nodules. The combination of multiple datasets was carried out in R environment (v3.4.3, http://www.r-project.org) (Sepulveda, 2020). The raw data (CEL file) of all datasets were read, log-transformed, and normalized using the corresponding R package, and all parameters were set as default. All probe sets were then mapped to their corresponding gene names using Bioconductor (Tippmann, 2015). The average expression value was retained if one gene was mapped to multiple probes (Yang et al., 2020c). To remove batch effects among five independent datasets, Z-score transformation was used to adjust the gene expression levels in each dataset (Yang Q et al., 2019b; Yang et al., 2020a). Z-score transformation for each gene could be computed by subtracting the mean of all genes and dividing the difference by the standard deviation of all genes in one experiment. After data transformation, the mean value for each experiment became zero with standard deviation equaling one.

In this study, there were five collected datasets integrated as a comprehensive dataset for discovering signatures. This comprehensive dataset consisted of 150 malignant and 93 benign samples of thyroid nodules. To the best of one’s knowledge, this integrated dataset was the largest transcriptomic dataset in the analysis of thyroid nodules. Based on this comprehensive dataset, the DEGs were discovered using feature selection methods including Student’s t test and fold change (FC). For Student’s t test, multtest package of R language was applied, and the adjusted p-value < 0.05 was selected as the cutoff (Yan et al., 2019). The fold change was used to compare the mean expression of each gene between malignant and benign thyroid nodules (Yu et al., 2020). The cutoff level of FC was set to logFC >0.58 (FC > 1.5) or logFC < -0.58 (FC < 0.67). The equation of FC was shown below (as shown in Eq. (1)). log ⁡ F C = m e a n ( log ⁡ 2 ( M a l i g n a n t   G r o u p ) ) − m e a n ( log ⁡ 2 ( B e n i g n   G r o u p ) ) . (Eq.1)

The volcano plot was applied to visualize and demonstrate the DEGs using ggplot2 package. Then the analysis of gene ontology (GO) enrichment was performed to identify the key biological processes for thyroid nodules (Yang et al., 2019a). Moreover, GOplot and clusterProfiler packages were used for visualizing the biological processes (BP) of GO enrichment (Yu et al., 2012; Yang et al., 2021). The raw p-value < 0.05 of GO terms was considered statistically significant.

The WGCNA package was applied to establish the scale-free weight gene co-expression networks for thyroid nodules (Langfelder and Horvath, 2008). The unqualified genes were screened out, and the matrix of genes’ similarity by Pearson’s correlation analysis was created. Appropriate soft threshold power (β) was applied to strengthen this matrix to a scale-free co-expression network (Yang et al., 2020b). The lowest power was chosen, so the scale-free topology fit index curve flattened out upon reaching a high value. The highly correlated genes were assigned into the same module. As a result, the intersection was obtained between DEGs identified by the feature selection method and hub genes in a key module using the WGCNA model. These genes in the intersection were regarded as the key genes for further validation.

A systematic validation was conducted by evaluating the upregulated and downregulated genes based on the independent dataset (GSE34289) (Alexander et al., 2012). This validation dataset consisted of two independent datasets from two different platforms. The first independent dataset was detected based on GPL5175 platform (Affymetrix Human Exon 1.0 ST Array). In this dataset, there were 23 malignant and 26 benign thyroid nodules. The second independent dataset was detected based on GPL14961 platform (Afirma-T Human Custom Array). There were 120 malignant and 198 benign samples in this second independent dataset. In this study, the boxplot was used to demonstrate the differential expression of these key genes between malignant and benign thyroid nodules.

To construct a classification model for thyroid nodules, four powerful classifiers, namely, support vector machine, linear discriminate analysis, partial least squares, and random forest algorithm, were applied in this study (Orru et al., 2012). The key genes between malignant and benign thyroid nodules were used to discriminate different samples. In the first step, the five-fold cross validation of the comprehensive dataset (Table 1) was performed to validate the performance of this classification model. The accuracy of five-fold cross validation could reflect the quality of the model. In the second step, the comprehensive dataset was set as the training set, and the two independent datasets from GSE34289 were set as the test sets. The performance of the independent test set could accurately reflect the classification ability of the model. This high-performance classification model based on machine learning was constructed for discriminating benign and malignant thyroid nodules.

A variety of microarray studies based on thyroid tissue were collected by searching the key word “thyroid nodules” in the GEO database. As a result, five independent transcriptomic studies were obtained, and each comprised a cohort of malignant samples and another cohort of benign samples. The detailed information of these independent datasets is provided in Table 1. Among these studies, the five datasets including 150 malignant and 93 benign thyroid nodules were combined as a comprehensive dataset. The boxplots of five datasets before and after batch effect removal are shown in Supplementary Figure S1. The intensity of all samples before batch effect removal was distributed in the range of 4–15 and fluctuated greatly. After batch effect removal, the intensity of all samples was roughly distributed in the range of -1–1. The stable distribution indicated that the batch effects were well removed in the combined dataset by Z-score transformation. After data preprocessing and batch effect removal, the comprehensive dataset with 7,265 genes from five independent studies was applied to discover the key genes of thyroid nodules.

Based on this comprehensive dataset, the DEGs were discovered using feature selection methods (both Student’s t test and fold change). The volcano plot (as shown in Figure 1) illuminated the variation of DEGs in malignant versus benign thyroid nodules. The horizon line was the cutoff (adjusted p-value < 0.05) of Student’s t test. The cutoff levels for the vertical line were set to logFC >0.58 (FC > 1.5) or logFC < -0.58 (FC < 0.67) of fold change. The blue and red dots were used to indicate the upregulated (logFC >0.58) and downregulated (logFC < -0.58) genes, respectively. In this study, 279 DEGs were finally identified by both Student’s t test and fold change. The total number of upregulated genes (172 genes) was larger than that of the downregulated ones (107 genes). The top 25 upregulated and downregulated DEGs are shown in Table 2, including the information of entrez ID, gene symbol, adjusted p-value, and fold change for each gene. The information of all DEGs is shown in Supplementary Table S1.

GO enrichment analysis is ubiquitously used for interpreting high throughput molecular data and underlying biological phenomena of experiments (Tomczak et al., 2018). For a set of genes, an enrichment analysis will find which GO terms are overrepresented using annotations for the gene set. GO enrichment analysis for the DEGs was performed in this study. Using the DEGs between malignant and benign thyroid nodules, the enrichment analysis included the BP (biological process), MF (molecular function), and CC (cell component) terms. The detailed information of GO ID, description, p-value, name, and the number of genes is shown in Supplementary Table S2.

Particularly, multiple biological processes were enriched to interpret the biological mechanism of malignant thyroid nodules. The chord diagram of BP enrichment (as interpreted in Figure 2) was applied to explain the relationship between DEGs and BP terms. It was reported that these BP terms were associated with the biological mechanism of thyroid nodules. For example, there were 15 DEGs enriched in the hormone metabolic process, and the association with thyroid cancer has been reported (Han et al., 2018). The platelet degranulation enriched by 10 DEGs was discovered in papillary thyroid carcinoma using the biomarkers (Wu et al., 2018). The concentration of the vascular endothelial growth factor was increased and stimulated endothelial cell proliferation in the cyst fluid of enlarging and recurrent thyroid nodules (Sato et al., 1997). It was reported that patients with spotty skin pigmentation had a predisposition toward the development of thyroid abnormalities (Courcoutsakis et al., 2009). It was found that low thyroid hormones might have implications for reproductive health, so the reproductive structure development and reproductive system development might be affected in thyroid nodules (Medda et al., 2017). The thyroid hormone generation reported that the significant biologic process was involved in thyroid cancers (Durante et al., 2018).

The WGCNA network was constructed to identify the gene co-expression module (as shown in Figure 3). The value of power (10) was selected as the soft-threshold power to ensure scale-free (R ² = 0.8) networks using the WGCNA package (Figure 3A) because it reached the plateau at power 10 from the scale-free topology plot and mean connectivity plot. Genes with similar expression patterns were clustered into co-expression modules. Different modules were shown in different colors, and 13 modules were identified totally (Figure 3B). The heatmap of module–trait relationships was applied for depicting correlations between module eigengenes and phenotypic traits (the label of malignant and benign thyroid nodules). As shown in Figure 3C, the numbers correspond to the correlation, and the p-values were set in parentheses. Moreover, the degree of correlation was illustrated with the color legend. Here, the gray module was the most correlated one with malignant thyroid nodules (R = 0.32, p-value = 2 × 10^–5). Hence, the gray module was used for the identification of the hub genes. Hub genes in the co-expression network were characterized by high intra-modular connectivity measured by the value of gene significance and module membership. The scatterplot of module eigengenes related to malignant thyroid nodules in the gray co-expression module (R = 0.29, p-value = 3 × 10^–10) is shown in Figure 3D. As a result, 454 genes in the gray module highly correlated with gene significance were identified as hub genes using WGCNA.

In this study, there were 19 overlapping genes in the intersection between 279 DEGs identified by the feature selection method and 454 hub genes in the gray module totally. To validate these overlapping genes, two independent datasets from GSE34289 were applied to perform the systematic validation (Alexander et al., 2012). In this validation dataset, there were 23 malignant with 26 benign samples and 120 malignant with 198 benign samples form GPL5175 and GPL14961 platforms, respectively. The boxplots (as shown in Figure 4) were used to demonstrate the key genes between malignant and benign thyroid nodules. Among the 19 overlapping genes, there were four key genes expressed in the independent dataset, and the dysregulation of these key genes was validated. As shown in Figure 4, the significant differences of three upregulated genes (ST3GAL5, NRCAM, and MT1F) and one downregulated gene (PROS1) were indicated in these boxplots obviously for the independent data detected from GPL5175 (Figure 4A) and GPL14961 platforms (Figure 4B), respectively.

As a result, these four key genes were effectively validated as the important ones participated in the pathogenesis of thyroid nodules. It was reported that the specific genetic variants of ST3GAL5 in patients with thyroid-associated ophthalmopathy were discovered (Park et al., 2017). Górka et al. provided the first evidence that NRCAM is overexpressed in papillary thyroid carcinomas, and the upregulation of NRCAM was implicated in the pathogenesis and behavior of papillary thyroid cancers (Gorka et al., 2007). It was reported that MT1F might contribute to thyroid carcinogenesis and potentially serve as a diagnostic marker in distinguishing benign from malignant lesions (Kim et al., 2010; Wojtczak et al., 2017). In the previous studies, PROS1 was reported as the biomarker significantly related to thyroid nodules’ malignancy (Griffith et al., 2006; Wu et al., 2020). In this study, these four key genes (ST3GAL5, NRCAM, MT1F, and PROS1) were discovered for distinguishing malignant from benign thyroid nodules.

To distinguish malignant from benign thyroid nodules, four popular machine learning methods were applied to construct the classification model in this study. These methods included support vector machine, linear discriminate analysis, partial least squares, and random forest algorithm. The key genes between benign and malignant thyroid nodules were used to discriminate different samples. For the comprehensive dataset in Table 1, the five-fold cross validation was first performed to validate the performance of this classification model. As shown in Figure 5A1, 5B1, 5C1, and 5D1, the values of area under the ROC curve (AUC) were 0.83, 0.82, 0.82, and 0.78 for the five-fold cross validation using four different machine learning methods, respectively. Moreover, the high performance of the independent test sets could accurately reflect the ability of the classification model. The comprehensive dataset was set as the training set, and the test sets consisted of two parts detected by GPL5175 and GPL14961 platforms from the independent dataset (GSE34289). As displayed in Figure 5A2, 5B2, 5C2, and 5D2, the AUC values of the ROC curve for the first independent test set were 0.83, 0.67, 0.74, and 0.74 by four machine learning methods, respectively. As shown in Figure 5A3, 5B3, 5C3, and 5D3, the AUC values for the second independent test set were 0.81, 0.60, 0.69, and 0.77 by four machine learning methods, respectively.

As shown in Figure 5, for the five-fold cross validation, the performances (AUC >0.8) of the classification model were outstanding using support vector machine, linear discriminate analysis, and partial least squares. However, the classification models of support vector machine and random forest (AUC >0.7) have shown more excellent performances than the other methods for the two independent test sets. Therefore, the high-performance classification model using support vector machine was recommended for discriminating malignant from benign thyroid nodules based on both five-fold cross validation and independent test.

Until now, it fails to discriminate as benign or malignant in one-third of thyroid nodules using FNA with cytologic evaluation. To save medical costs and improve the diagnostic accuracy, the high-performance classification model constructed in this study could be applied before FNA. For the thyroid nodule patients, the expression of four key genes could be detected. Then, this sample could be classified as benign or malignant thyroid nodules based on the classification model. If the patient was classified as a malignant thyroid sample, it was highly necessary to make a definite diagnosis using FNA with cytologic evaluation. If the patient was classified as a benign sample based on the classification model, the necessity of the FNA could be determined depending on the specific conditions. In the future, selection method, the high-performance classification model is expected to be applied for clinical diagnosis and management for malignant and benign thyroid nodules.