Identifying key genes of classic papillary thyroid cancer in women aged more than 55 years old using bioinformatics analysis

Background: The incidence rate of thyroid carcinoma (THCA) markedly increased in the recent few decades and has been likely over-diagnosed, especially papillary thyroid cancer (PTC) in women. However, the incidence of advanced-stage papillary thyroid cancer is also rising. According to earlier studies, tumors with identical pathology might have different clinical outcomes, which implies some variances in papillary thyroid cancer. Although the mortality of thyroid cancer has remained stable or declined, there is still an important problem in estimating whether it is benign or needs surgery for patients with papillary thyroid cancer.

Methods: After obtaining data from The Cancer Genome Atlas (TCGA) Project-THCA database by R package TCGA bio links, 18 samples (11 at stage IV as high-risk group and 7 at stage I as low-risk group) were obtained using survival package and edgeR to ensure differential expression; ClusterProfiler package was used to carry on gene set enrichment analysis and searched the possible pathways in the Kyoto Encyclopedia of Genes and Genomes (KEGG) database. STRING and Cytoscape were used to construct and modify the protein–protein interaction (PPI) network to get hub genes of differentially expressed genes. Next, the pROC package was used to get the receiver operating characteristic (ROC) curves of hub genes’ disease-free survival (DFS). Then, transcription factors (TFs) and miRNAs of key genes were predicted by ENCORI and AnimalTFDB. In the end, TF–target genes–miRNA regulatory network was also constructed by Cytoscape.

Results: Our research obtained the top 9 candidate genes from the whole network (IFNA1, MRC1, LGALS3, LOX, POSTN, TIMP1, CD276, SDC4, and TLR2). According to the ROC results, TIMP1, LOX, CD276, IFNA1, TLR2, and POSTN were considered to play a more critical role in malignant papillary thyroid cancer or immature cancer of papillary thyroid cancer. Our analysis concludes that TIMP1, LOX, CD276, IFNA1, TLR2, and POSTN are identified as thyroid cancer biomarkers, which lead to the different clinical courses of a woman older than 55 years old with papillary thyroid cancer. Especially CD276, POSTN, and IFNA1 may be considered as new biomarkers associated with the prognosis of thyroid cancer.

Conclusions: TIMP1, LOX, CD276, IFNA1, TLR2, and POSTN have different expressions in PTCs, which lead to the various clinical courses of a woman older than 55 years old with papillary thyroid cancer. Especially CD276, POSTN, and IFNA1 may be considered as new potential biomarkers associated with the prognosis of thyroid cancer. In addition, TF–miRNA–target gene regulatory network may help further reach for PTC.

Introduction

The incidence rate of thyroid carcinoma (THCA) markedly increased in the recent few decades (1) and has been likely over-diagnosed because ultrasonography or other modern diagnostic techniques have been widely used to discover previously undetectable thyroid nodules (2). Differentiated thyroid cancer is the most common thyroid cancer, particularly the category of papillary thyroid cancer (PTC) in women diagnosed by methods such as fine-needle aspiration (FNA) after imaging tests. However, the incidence of advanced-stage PTC (large tumors and locally advanced and/or metastatic) rises (3). The fact that tumors with identical pathologies might have distinct clinical outcomes suggests some variances in PTC, according to previous studies (4, 5). Although the mortality of THCA has remained stable or declined, there is still an important problem in estimating if it is benign or needs surgery for patients with PTC.

Controversial issues in thyroid cancer management mainly focus on treatment options for different thyroid cancers (DTCs) (6). Recently, detecting genetic changes in thyroid cancer has been generally stressed and applied as an essential means to guide treatment (7). DNA sequencing technology has been frequently used to detect a different group of diagnoses. Bioinformatics technology has also been employed to identify those different gene expressions (DEGs). The BRAF and RAS family’s mutations are usually found in THCA and mitogen-activated protein kinase (MAPK) cellular signaling pathway and PI3K/AKT/mTOR pathway (8, 9). Those findings have greatly improved our understanding of studying THCA with different clinical outcomes.

This research performed many analyses through a wide-ranging bioinformatics study to classify the critical hub gene in thyroid cancer with a poor prognosis. A total of 994 genes upregulated in the high-risk group were obtained from TCGA-THCA data by edgeR (10, 11). Then, GO and KEGG pathway enrichment and the construction of protein–protein interaction (PPI) network and survival analysis were performed in these differentially expressed genes (DEGs) (12–14). Finally, key genes, transcription factors, and miRNA relations that might be important in THCA with poor prognosis were identified.

Methods

Datasets selection and DEGs identification

The Cancer Genome Atlas (TCGA) research network has collected many public clinical and molecular profiling of more than 10,000 tumor patients across 33 different tumor types. We obtained the clinical and transcriptome profiling of THCA from the TCGA-THCA database by TCGA-biolinks packages, which provide several useful functions to search, download, and prepare TCGA samples for data analysis (15), using the following criteria: (A) PTC-classic/usual, (B) female, (C) age ≥55 years old, (D) datasets including stage I and IV samples, and (E) replicated samples value being saved with unique bcr patient barcode. Next, Cox proportional hazards model was performed by survival packages of R for those data, and the result was presented by forest plot packages (16, 17). Furthermore, the edgeR packages were used to distinguish the DEGs in each sample in R (Version 4.1.2) (10, 11). The standard for DEGs has been considered a false discovery rate of <0.05. The groups’ criteria were |log2FC (fold change) | ≥ 1.

GO and KEGG enrichment analyses of DEGs

Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) were performed to reveal the functional enrichment analysis of DEGs. BiomaRt packages that include over 800 different biological datasets spanning were used to convert gene ID from other types (18, 19). ClusterProfiler packages that automate enrichment analysis of gene clusters were used to progress GO enrichment analyses, which include cellular component (CC), biological process (BP), and molecular function (MF), and KEGG pathways analysis (20, 21). Finally, GOplot and ggplot2 packages were used to make figures of those results (22, 23).

PPI network construction

STRING (https://string-db.org/), an online database, was used to establish the protein–protein interaction (PPI) network of the candidate genes, which evaluates the interacting units with protein-coding gene loci (represented by their main, canonical protein isoform) (24), and the minimum required score was set as the medium confidence (0.400). Moreover, Cytoscape (Version 3.9.1) was used to visualize, model, and analyze molecular and genetic interaction networks (25). Next, the analysis network of Cytoscape tool, Molecular Complex Detection (MCODE), was applied to screen modules of the PPI network with a degree cutoff = 2, node score cutoff = 0.2, k-core = 2, and max. depth = 100.

Survival analysis

The receiver operating characteristic (ROC) curve is a commonly used graphical summary to evaluate the predictive value of biomarkers (26). To verify the prognosis results of hub genes, disease-free survival (DFS) data obtained from TCGA and pROC package were used to make the ROC curve figure (27).

TF and miRNA regulatory network construction

After the reference sequence (RefSeq) of the key genes were obtained from National Center for Biotechnology Information (NCBI) RefSeq database (https://www.ncbi.nlm.nih.gov/), transcription factors were predicted by AnimalTFDB (https://bioinfo.life.hust.edu.cn/AnimalTFDB/#!/) and selected by score >20 (28, 29). ENCORI (https://rna.sysu.edu.cn/encori/index.php) was used to predict miRNAs that bind to hub genes and apply standards including CLIP-Data ≥3, pan-cancer ≥1, and miRNA, with the least intersections in three databases selected (30).

Results

Patient characteristics

The clinical characteristics of 359 cases of classic PTC were obtained from the TCGA-THCA project, including days to the last follow-up, days to death, disease-free status, age at diagnosis, TNM classification, and tumor stage. After Cox analysis, age at diagnosis and tumor stage were significant for disease-free status (Figure 1A). Compared to age <55 years, the order indicated a poor prognosis that corresponds with previous research. At the same time, ages between 55 and 65 still have no noticeable difference with age <55 years old, suggesting that there might some of the patients with a preferable prognosis. Consequently, with age more than 55 years, 7 cases of stage I as a low-risk group and 11 cases of stage IV as a high-risk group were finally obtained after filtering the following criteria. The different expressions between the two groups were analyzed by edgeR and shown in the volcano plot (Figure 1B).

FIGURE 1

Figure 1 Clinical characteristics of classic PTC and volcano plot of DEGs. (A) Age, stage, and tumor size are significant. The reference of age is <55 years old. The reference of stage is stage I, and references of TNM are T1, N0, and M0. (B) Volcano plot of DEGs. Those blue dots represent downregulation in high-risk group thyroid carcinoma compared with low-risk group samples; red represents upregulation, and gray dots show no significant difference in the two groups.

Function annotation of DEGs

GO enrichment analysis was completed to explore the biological functions of all DEGs of the upregulated gene obtained after edgeR. The GO term biological process analysis corroborated that the DEGs were enriched in external encapsulating structure organization, extracellular matrix organization, and extracellular structure organization. The GO term cellular component was enriched in the collagen-containing extracellular matrix, cornified envelope, and endoplasmic reticulum lumen. Apart from that, the molecular function classification of GO was enriched in extracellular matrix structural constituent, integrin binding, and cytokine activity (Figure 2).

FIGURE 2

Figure 2 Gene Ontology analyses of differentially expressed genes. Biological process (BP) was enriched in the external encapsulating structure organization, extracellular matrix organization, and extracellular structure organization. The cellular component (CC) was enriched in the collagen-containing extracellular matrix, cornified envelope, and endoplasmic reticulum lumen. Molecular function (MF) was enriched in extracellular matrix structural constituent, integrin binding, and cytokine activity.

The enrichment analysis of the KEGG pathway demonstrated that DEGs were correlated with proteoglycans in cancer, PI3K-Akt signaling pathway, and p53 signaling pathway (Figure 3).

FIGURE 3

Figure 3 KEGG pathways of differentially expressed genes. The KEGG pathway could find proteoglycans in cancer, PI3K-Akt signaling pathway, and p53 signaling pathway related to thyroid carcinoma.

Construction of PPI network

The STRING database identified the interaction of these upregulated DEGs. The combined score >0.4 (medium confidence score) was considered statistically significant. Then, a PPI network was made with Cytoscape. A total of 905 nodes and 4,941 edges were included in that PPI network. Additionally, the modules in the PPI network were analyzed by MCODE, and >8 was set as cutoff criteria with the default parameters (degree cutoff, 2; node score cutoff, 2; K-score, 2; and max depth = 100) (Figure 4). Next, the top 3 genes of every cluster (IFNA1, MRC1, LGALS3, LOX, POSTN, TIMP1, CD276, SDC4, and TLR2) were obtained as hub genes after being ranked by MCODE score.

FIGURE 4

Figure 4 MCODE cluster. Cluster calculated by MCODE after protein–protein interaction network constructed, and different sizes and colors represent different MCODE scores. (A) Cluster 1 (score, 10.1). (B) Cluster 2 (score, 9.8). (C) Cluster 3 (score, 8.5).

Survival analysis

pROC packages made the receiver operating characteristic (ROC) curve to explore the patient’s prognosis with high-expression hub genes. According to the results of the ROC curve, TIMP1, LOX, CD276, IFNA1, TLR2, and POSTN as key genes are more relevant to THCA (Figure 5). The heatmap of key gene expression in sample patients was made after normalization. Compared to IFNA1, high expressions of other genes are more common in the high-risk group (Figure 6).

FIGURE 5

Figure 5 Receiver operating characteristic (ROC) curve. Receiver operating characteristic (ROC) curve made by pROC. TIMP1, LOX, CD276, IFNA1, TLR2, and POSTN as key genes are more relevant to THCA.

FIGURE 6

Figure 6 Heatmap of key genes. Heatmap made by pheatmap package and normalized by scale function. Different expression of genes in the high- and low-risk groups.

TF and miRNA regulatory network

A total of 47 miRNAs and 44 transcription factors that could be bound to key genes are predicted. TF, miRNA, and target gene relations were revealed in the regulatory network (Figure 7).

FIGURE 7

Figure 7 TF–miRNA–target gene network. A total of 47 miRNAs and 44 transcription factors that could be bound to key genes are predicted.

Discussion

In recent years, thyroid carcinoma has had a worldwide and dramatic incidence, especially in women. However, thyroid cancer treatment options are limited, and they mostly need to have surgery and supply of thyroxine for the whole of their later life. The high-intensity treatment for most tumor diameters over 1 cm, especially patients with differentiated thyroid cancer or papillary thyroid cancer, has also been controversial after PTC prognosis is more remarkable than that of other tumors. The fact that advanced thyroid cancer is found with the increase in cardinality estimate also promotes more research to distinguish the difference from PTC so that more suitable and precise therapeutics could be suggested to each of the patients.

For the above target, it is critical to know the prognosis of carcinoma, especially for PTC, which has a high risk of incidence and is likely over-diagnosed. Many studies have utilized bioinformatics technology to uncover biomarkers in THCA or PTC. However, the majority of them compared carcinoma and normal samples. According to the 8th American Joint Committee on Cancer (AJCC) and other experts, there may be distinct influencing factors for varying clinical outcomes in patients with PTC who are over 55 years old (31–34). In the current study, based on the hypothesis of difference in PTC, we explore the potential difference between PTC classification in women over 55 years old who have the worst prognosis and those with favorable forecasts based on data from TCGA.

In the present study, 994 upregulated genes were obtained from TCGA-THCA data using the edgeR package. GO was enriched in cell connection and extracellular matrix organization. KEGG pathway enrichment analysis of those DEGs could be linked to thyroid cancer growth and metastasis pathways such as proteoglycans in cancer, PI3K-Akt signaling pathway, and p53 signaling pathway. After using MCODE and Cytoscape network analysis tools of Cytoscape, nine hub genes with greater degrees were finally obtained (IFNA1, MRC1, LGALS3, LOX, POSTN, TIMP1, CD276, SDC4, and TLR2), which were considerably overexpressed in PTC with a poor prognosis compared to those with a better prognosis. In comparison to the survival analysis of nine hub genes, TIMP1, LOX, CD276, IFNA1, TLR2, and POSTN as key genes are more relevant to THCA.

Tissue inhibitor of metalloproteinase 1 (TIMP1) belongs to the TIMP gene family, and its encoded protein can promote cell proliferation. TIMP1 was found to be overexpressed in classic and follicular variants of PTC in different age and gender groups and in the BRAF-MUT group compared to that in BRAF-WT patients in previous research (35, 36). It binds CD63 on the cell surface membrane and activates AKT signaling pathway to make antiapoptotic activity and predict aggressive behaviors (37, 38). Lysyl oxidase (LOX) was found to be overexpressed in aggressive cancers and related to MMPs and TIMPs by regulating SNA12 expression (39). In addition, LOX upregulation is also associated with anaplastic thyroid cancer progression and aggressive tall cell variant of PTC (TC-PTC) compared with the differentiated thyroid cancers [classic PTC (cPTC) and follicular variant of PTC (FV-PTC)] that may respond to BRAF activation (40–42). Toll-like receptor 2 (TLR2) is a member of the Toll-like receptor (TLR) family; it was found to be mediated by Akt phosphorylation by GT1b (trisialoganglioside 1b) and activated PI3K/Akt signaling pathway (43). POSTN-encoded protein binds to integrins to support adhesion, migration of epithelial cells, and participation in cancer stem cell maintenance and metastasis. Previous research once found a wide variability of POSTN expression in a different histological subtype of PTC, and high stromal POSTN expression is associated with aggressive tumor behavior (44). Compared to TIMP1 and LOX, the relation between TLR2 and thyroid carcinoma was unclear, especially in PTC. However, previous research indicated that overexpression of TLR2 was associated with a high risk of colorectal cancer, cervical cancer, and oral cancer (45–48). CD276, also known as B7-H3, has upregulated expression in many cancer cells and related with tumor cells by pathological angiogenesis (49). CD276 has been found to have an upregulated expression at the protein and mRNA levels in different types of thyroid cancer, and its effect on thyroid carcinoma has also been researched recently. CD276 was overexpressed in advanced thyroid carcinoma such as ATC, which is consistent with our results; its high expression indicates a poor prognosis of different PTCs (50, 51). We found that in the classic PTC, the upregulated CD276 is also related to a poor prognosis. Periostin (POSTN) encoded a secreted extracellular matrix protein and is found to exhibit stromal deposition in invasive portions and cytoplasmic expression of undifferentiated thyroid carcinoma cells (52). Periostin has different expressions in different kinds of PTCs, and upregulated periostin could also be found in PTC compared to non-neoplastic tissues. Upregulation of periostin-introduced invasion and lymph node metastasis may be due to loss of cellular polarity/cohesiveness and epithelial–mesenchymal transition (EMT) (53). This research also was consistent with the above results and found that in classic PTC of female patients aged more than 55 years, POSTN still has a different expression and upregulated POSTN still indicated a poor prognosis. Interferon alpha 1 (IFNA1) could encode a member of the type I interferon, and expression of IFNA1 is associated with HPV-positive head and neck cancers (54). Interferon alpha 1 could be found in normal thyroid cells and differentiated thyroid cancer and also used in some diseases, but long-term interferon (IFN) therapy is frequently associated with side effects on the thyroid gland (55). Although its effect is considered antiproliferative, interferon alpha is also found to influence malignant thyroid cancer by the microenvironment or affect other factors, such as tumor necrosis factor (TNF) or epidermal growth factor (EGF) (56, 57). According to our research, IFNA affects high-risk group PTC with an upregulated expression. More facts about its effect on PTC also need to be explored by further research. Above all, in our results, TIMP1, LOX, and POSTEN had different expressions no matter what type of thyroid cancer and papillary thyroid cancer; even in pure classic PTC, high expression of those genes also predicted a poor prognosis. Although CD276, TLR2, and IFNA1 were not found to have complex functions in thyroid cancer, they all participate in different carcinomas. These results indicate that further study is needed to reveal their roles in thyroid carcinoma.

Data availability statement

The original contributions presented in the study are included in the article/supplementary material. Further inquiries can be directed to the corresponding author.

Author contributions

L-QY designed the analyses. C-CL and MU analyzed the data and wrote the manuscript. XL and S-KS collected the data, YW and FL prepared the figures, and BG and M-HZ revised the manuscript. All authors contributed to the article and approved the submitted version.

Funding

This work was supported by grants from the National Natural Science Foundation of China (Nos. 82100944, 82100494, 82070910, and 81870623), Natural Science Foundation of Hunan Province (No. 2021JJ40842), and Key R&D Plan Hunan Province (2020SK2078).

Acknowledgments

We want to thank our colleagues in the lab for their assistance.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

1. Kitahara CM, Sosa JA. The changing incidence of thyroid cancer. Nat Rev Endocrinol (2016) 12(11):646–53. doi: 10.1038/nrendo.2016.110

CrossRef Full Text | Google Scholar

2. Vaccarella S, Lortet-Tieulent J, Colombet M, Davies L, Stiller CA, Schüz J, et al. Global patterns and trends in incidence and mortality of thyroid cancer in children and adolescents: A population-based study. Lancet Diabetes Endocrinol (2021) 9(3):144–52. doi: 10.1016/S2213-8587(20)30401-0

CrossRef Full Text | Google Scholar

3. Sung H, Ferlay J, Siegel RL, Laversanne M, Soerjomataram I, Jemal A, et al. Global cancer statistics 2020: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA: A Cancer J Clin (2021) 71(3):209–49. doi: 10.3322/caac.21660

CrossRef Full Text | Google Scholar

4. Takano T. Natural history of thyroid cancer. Endocr J (2017) 64(3):237–44. doi: 10.1507/endocrj.EJ17-0026

CrossRef Full Text | Google Scholar

5. Williams D. Thyroid growth and cancer. Eur Thyroid J (2015) 4(3):164–73. doi: 10.1159/000437263

CrossRef Full Text | Google Scholar

6. Tuttle RM. Controversial issues in thyroid cancer management. J Nucl Med (2018) 59(8):1187–94. doi: 10.2967/jnumed.117.192559

CrossRef Full Text | Google Scholar

7. Wang TS, Sosa JA. Thyroid surgery for differentiated thyroid cancer - recent advances and future directions. Nat Rev Endocrinol (2018) 14(11):670–83. doi: 10.1038/s41574-018-0080-7

CrossRef Full Text | Google Scholar

8. Cabanillas ME, McFadden DG, Durante C. Thyroid cancer. Lancet (2016) 388(10061):2783–95. doi: 10.1016/S0140-6736(16)30172-6

CrossRef Full Text | Google Scholar

9. Cohen-Solal KA, Boregowda RK, Lasfar A. RUNX2 and the PI3K/AKT axis reciprocal activation as a driving force for tumor progression. Mol Cancer (2015) 14:137. doi: 10.1186/s12943-015-0404-3

CrossRef Full Text | Google Scholar

10. Robinson MD, McCarthy DJ, Smyth GK. edgeR: a bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics (2010) 26(1):139–40. doi: 10.1093/bioinformatics/btp616

CrossRef Full Text | Google Scholar

11. McCarthy DJ, Chen Y, Smyth GK. Differential expression analysis of multifactor RNA-seq experiments with respect to biological variation. Nucleic Acids Res (2012) 40(10):4288–97. doi: 10.1093/nar/gks042

CrossRef Full Text | Google Scholar

12. Thomas PD. The gene ontology and the meaning of biological function. Methods Mol Biol (2017) 1446:15–24. doi: 10.1007/978-1-4939-3743-1_2

CrossRef Full Text | Google Scholar

13. Xing Z, Chu C, Chen L, Kong X. The use of gene ontology terms and KEGG pathways for analysis and prediction of oncogenes. Biochim Biophys Acta (2016) 1860(11 Pt B):2725–34. doi: 10.1016/j.bbagen.2016.01.012

CrossRef Full Text | Google Scholar

14. Khunlertgit N, Yoon BJ. Incorporating topological information for predicting robust cancer subnetwork markers in human protein-protein interaction network. BMC Bioinf (2016) 17(Suppl 13):351. doi: 10.1186/s12859-016-1224-1

CrossRef Full Text | Google Scholar

15. Colaprico A, Silva TC, Olsen C, Garofano L, Cava C, Garolini D, et al. TCGAbiolinks: an R/Bioconductor package for integrative analysis of TCGA data. Nucleic Acids Res (2016) 44(8):e71–e. doi: 10.1093/nar/gkv1507

CrossRef Full Text | Google Scholar

16. T T. a package for survival analysis in r (2022). Available at: https://CRAN.R-project.org/package=survival.

Google Scholar

17. Lumley MGAT. Forestplot: Advanced forest plot using 'Grid' graphics (2021). Available at: https://CRAN.R-project.org/package=forestplot.

Google Scholar

18. Smedley D, Haider S, Durinck S, Pandini L, Provero P, Allen J, et al. The BioMart community portal: An innovative alternative to large, centralized data repositories. Nucleic Acids Res (2015) 43(W1):W589–98. doi: 10.1093/nar/gkv350

CrossRef Full Text | Google Scholar

19. Durinck S, Moreau Y, Kasprzyk A, Davis S, De Moor B, Brazma A, et al. BioMart and bioconductor: A powerful link between biological databases and microarray data analysis. Bioinformatics (2005) 21(16):3439–40. doi: 10.1093/bioinformatics/bti525

CrossRef Full Text | Google Scholar

20. Yu G, Wang LG, Han Y, He QY. clusterProfiler: An r package for comparing biological themes among gene clusters. OMICS (2012) 16(5):284–7. doi: 10.1089/omi.2011.0118

CrossRef Full Text | Google Scholar

21. Wu T, Hu E, Xu S, Chen M, Guo P, Dai Z, et al. clusterProfiler 4.0: A universal enrichment tool for interpreting omics data. Innovation (N Y) (2021) 2(3):100141. doi: 10.1016/j.xinn.2021.100141

CrossRef Full Text | Google Scholar

22. Walter W, Sanchez-Cabo F, Ricote M. GOplot: an r package for visually combining expression data with functional analysis. Bioinformatics (2015) 31(17):2912–4. doi: 10.1093/bioinformatics/btv300

CrossRef Full Text | Google Scholar

23. Wickham H. ggplot2: Elegant graphics for data analysis. New York: Springer-Verlag (2016).

Google Scholar

24. Szklarczyk D, Morris JH, Cook H, Kuhn M, Wyder S, Simonovic M, et al. The STRING database in 2017: Quality-controlled protein-protein association networks, made broadly accessible. Nucleic Acids Res (2017) 45(D1):D362–D8. doi: 10.1093/nar/gkw937

CrossRef Full Text | Google Scholar

25. Cline MS, Smoot M, Cerami E, Kuchinsky A, Landys N, Workman C, et al. Integration of biological networks and gene expression data using cytoscape. Nat Protoc (2007) 2(10):2366–82. doi: 10.1038/nprot.2007.324

CrossRef Full Text | Google Scholar

26. Michael H, Tian L, Ghebremichael M. The ROC curve for regularly measured longitudinal biomarkers. Biostatistics (2019) 20(3):433–51. doi: 10.1093/biostatistics/kxy010

CrossRef Full Text | Google Scholar

27. Robin X, Turck N, Hainard A, Tiberti N, Lisacek F, Sanchez JC, et al. pROC: an open-source package for r and s+ to analyze and compare ROC curves. BMC Bioinf (2011) 12:77. doi: 10.1186/1471-2105-12-77

CrossRef Full Text | Google Scholar

28. Hu H, Miao YR, Jia LH, Yu QY, Zhang Q, Guo AY. AnimalTFDB 3.0: a comprehensive resource for annotation and prediction of animal transcription factors. Nucleic Acids Res (2019) 47(D1):D33–D8. doi: 10.1093/nar/gky822

CrossRef Full Text | Google Scholar

29. O'Leary NA, Wright MW, Brister JR, Ciufo S, Haddad D, McVeigh R, et al. Reference sequence (RefSeq) database at NCBI: Current status, taxonomic expansion, and functional annotation. Nucleic Acids Res (2016) 44(D1):D733–45. doi: 10.1093/nar/gkv1189

CrossRef Full Text | Google Scholar

30. Li JH, Liu S, Zhou H, Qu LH, Yang JH. starBase v2.0: decoding miRNA-ceRNA, miRNA-ncRNA and protein-RNA interaction networks from large-scale CLIP-seq data. Nucleic Acids Res (2014) 42(Database issue):D92–7. doi: 10.1093/nar/gkt1248

CrossRef Full Text | Google Scholar

31. Uchino S, Ishikawa H, Miyauchi A, Hirokawa M, Noguchi S, Ushiama M, et al. Age- and gender-specific risk of thyroid cancer in patients with familial adenomatous polyposis. J Clin Endocrinol Metab (2016) 101(12):4611–7. doi: 10.1210/jc.2016-2043

CrossRef Full Text | Google Scholar

32. Ito Y, Miyauchi A, Oda H. Low-risk papillary microcarcinoma of the thyroid: A review of active surveillance trials. Eur J Surg Oncol (2018) 44(3):307–15. doi: 10.1016/j.ejso.2017.03.004

CrossRef Full Text | Google Scholar

33. Ito Y, Miyauchi A, Kihara M, Higashiyama T, Kobayashi K, Miya A. Patient age is significantly related to the progression of papillary microcarcinoma of the thyroid under observation. Thyroid (2014) 24(1):27–34. doi: 10.1089/thy.2013.0367

CrossRef Full Text | Google Scholar

34. Kazaure HS, Roman SA, Sosa JA. The impact of age on thyroid cancer staging. Curr Opin Endocrinol Diabetes Obes (2018) 25(5):330–4. doi: 10.1097/MED.0000000000000430

CrossRef Full Text | Google Scholar

35. Smallridge RC, Chindris AM, Asmann YW, Casler JD, Serie DJ, Reddi HV, et al. RNA Sequencing identifies multiple fusion transcripts, differentially expressed genes, and reduced expression of immune function genes in BRAF (V600E) mutant vs BRAF wild-type papillary thyroid carcinoma. J Clin Endocrinol Metab (2014) 99(2):E338–47. doi: 10.1210/jc.2013-2792

CrossRef Full Text | Google Scholar

36. Jarzab B, Wiench M, Fujarewicz K, Simek K, Jarzab M, Oczko-Wojciechowska M, et al. Gene expression profile of papillary thyroid cancer: Sources of variability and diagnostic implications. Cancer Res (2005) 65(4):1587–97. doi: 10.1158/0008-5472.CAN-04-3078

CrossRef Full Text | Google Scholar

37. Bommarito A, Richiusa P, Carissimi E, Pizzolanti G, Rodolico V, Zito G, et al. BRAFV600E mutation, TIMP-1 upregulation, and NF-kappaB activation: closing the loop on the papillary thyroid cancer trilogy. Endocr Relat Cancer (2011) 18(6):669–85. doi: 10.1530/ERC-11-0076

CrossRef Full Text | Google Scholar

38. Ilie MI, Lassalle S, Long-Mira E, Hofman V, Zangari J, Benaim G, et al. In papillary thyroid carcinoma, TIMP-1 expression correlates with BRAF (V600E) mutation status and together with hypoxia-related proteins predicts aggressive behavior. Virchows Arch (2013) 463(3):437–44. doi: 10.1007/s00428-013-1453-x

CrossRef Full Text | Google Scholar

39. Boufraqech M, Zhang L, Nilubol N, Sadowski SM, Kotian S, Quezado M, et al. Lysyl oxidase (LOX) transcriptionally regulates SNAI2 expression and TIMP4 secretion in human cancers. Clin Cancer Res (2016) 22(17):4491–504. doi: 10.1158/1078-0432.CCR-15-2461

CrossRef Full Text | Google Scholar

40. Boufraqech M, Nilubol N, Zhang L, Gara SK, Sadowski SM, Mehta A, et al. miR30a inhibits LOX expression and anaplastic thyroid cancer progression. Cancer Res (2015) 75(2):367–77. doi: 10.1158/0008-5472.CAN-14-2304

CrossRef Full Text | Google Scholar

41. Jolly LA, Novitskiy S, Owens P, Massoll N, Cheng N, Fang W, et al. Fibroblast-mediated collagen remodeling within the tumor microenvironment facilitates progression of thyroid cancers driven by BrafV600E and pten loss. Cancer Res (2016) 76(7):1804–13. doi: 10.1158/0008-5472.CAN-15-2351

CrossRef Full Text | Google Scholar

42. Hebrant A, Floor S, Saiselet M, Antoniou A, Desbuleux A, Snyers B, et al. miRNA expression in anaplastic thyroid carcinomas. PLoS One (2014) 9(8):e103871. doi: 10.1371/journal.pone.0103871

CrossRef Full Text | Google Scholar

43. Yoo HK, Park H, Hwang HS, Kim HJ, Choi Y-H, Kook KH. Ganglioside GT1b increases hyaluronic acid synthase 2 via PI3K activation with TLR2 dependence in orbital fibroblasts from thyroid eye disease patients. BMB Rep (2021) 54(2):136–41. doi: 10.5483/BMBRep.2021.54.2.178

CrossRef Full Text | Google Scholar

44. Giusca SE, Amalinei C, Lozneanu L, Ciobanu Apostol D, Andriescu EC, Scripcariu A, et al. Heterogeneous periostin expression in different histological variants of papillary thyroid carcinoma. BioMed Res Int (2017) 2017:8701386. doi: 10.1155/2017/8701386

CrossRef Full Text | Google Scholar

45. Proenca MA, de Oliveira JG, Cadamuro AC, Succi M, Netinho JG, Goloni-Bertolo EM, et al. TLR2 and TLR4 polymorphisms influence mRNA and protein expression in colorectal cancer. World J Gastroenterol (2015) 21(25):7730–41. doi: 10.3748/wjg.v21.i25.7730

CrossRef Full Text | Google Scholar

46. Zhu L, Yuan H, Jiang T, Wang R, Ma H, Zhang S. Association of TLR2 and TLR4 polymorphisms with risk of cancer: A meta-analysis. PLoS One (2013) 8(12):e82858. doi: 10.1371/journal.pone.0082858

CrossRef Full Text | Google Scholar

47. Yang S, Liu L, Xu D, Li X. The relationship of the TLR9 and TLR2 genetic polymorphisms with cervical cancer risk: A meta-analysis of case-control studies. Pathol Oncol Res (2020) 26(1):307–15. doi: 10.1007/s12253-018-0465-x

CrossRef Full Text | Google Scholar

48. Hussaini HM, Parachuru VPB, Seymour GJ, Rich AM. FoxP3 and TLR2 in co-expression in oral cancer. Acta Histochem (2017) 119(7):768. doi: 10.1016/j.acthis.2017.06.006

CrossRef Full Text | Google Scholar

49. Seaman S, Zhu Z, Saha S, Zhang XM, Yang MY, Hilton MB, et al. Eradication of tumors through simultaneous ablation of CD276/B7-H3-Positive tumor cells and tumor vasculature. Cancer Cell (2017) 31(4):501–15.e8. doi: 10.1016/j.ccell.2017.03.005

CrossRef Full Text | Google Scholar

50. Luo Y, Yang YC, Shen CK, Ma B, Xu WB, Wang QF, et al. Immune checkpoint protein expression defines the prognosis of advanced thyroid carcinoma. Front Endocrinol (Lausanne) (2022) 13:859013. doi: 10.3389/fendo.2022.859013

CrossRef Full Text | Google Scholar

51. Zhao B, Huang Z, Zhu X, Cai H, Huang Y, Zhang X, et al. Clinical significance of the expression of Co-stimulatory molecule B7-H3 in papillary thyroid carcinoma. Front Cell Dev Biol (2022) 10:819236. doi: 10.3389/fcell.2022.819236

CrossRef Full Text | Google Scholar

52. Kusafuka K, Yamashita M, Iwasaki T, Tsuchiya C, Kubota A, Hirata K, et al. Periostin expression and its supposed roles in benign and malignant thyroid nodules: An immunohistochemical study of 105 cases. Diagn Pathol (2021) 16(1):86. doi: 10.1186/s13000-021-01146-8

CrossRef Full Text | Google Scholar

53. Bai Y, Kakudo K, Nakamura M, Ozaki T, Li Y, Liu Z, et al. Loss of cellular polarity/cohesiveness in the invasive front of papillary thyroid carcinoma and periostin expression. Cancer Lett (2009) 281(2):188–95. doi: 10.1016/j.canlet.2009.02.043

CrossRef Full Text | Google Scholar

54. Zhang J, Chen T, Yang X, Cheng H, Spath SS, Clavijo PE, et al. Attenuated TRAF3 fosters activation of alternative NF-kappaB and reduced expression of antiviral interferon, TP53, and RB to promote HPV-positive head and neck cancers. Cancer Res (2018) 78(16):4613–26. doi: 10.1158/0008-5472.CAN-17-0642

CrossRef Full Text | Google Scholar

55. Selzer E, Wilfing A, Sexl V, Freissmuth M. Effects of type I-interferons on human thyroid epithelial cells derived from normal and tumour tissue. Naunyn Schmiedebergs Arch Pharmacol (1994) 350(3):322–8. doi: 10.1007/BF00175039

CrossRef Full Text | Google Scholar

56. Lahat N, Sheinfeld M, Sobel E, Kinarty A, Kraiem Z. Divergent effects of cytokines on human leukocyte antigen-DR antigen expression of neoplastic and non-neoplastic human thyroid cells. Cancer (1992) 69(7):1799–807. doi: 10.1002/1097-0142(19920401)69:7<1799::AID-CNCR2820690723>3.0.CO;2-8

CrossRef Full Text | Google Scholar

57. Cunha LL, Domingues GAB, Morari EC, Soares FA, Vassallo J, Ward LS. The immune landscape of the microenvironment of thyroid cancer is closely related to differentiation status. Cancer Cell Int (2021) 21(1):387. doi: 10.1186/s12935-021-02084-7

CrossRef Full Text | Google Scholar