Construction of ceRNA Network Based on TCGA and UCSC Xena Ocial Website to Provide Certain Help for Immunotherapy of Thyroid Cancer.

Background: Thyroid cancer is the most common endocrine tumor, and papillary thyroid carcinoma is the most common. Commonly used treatment methods include surgery and radiotherapy. As an emerging treatment option in recent years, immunotherapy has played an essential role in the treatment of many tumors, but there are few studies in thyroid cancer. Method: We rst downloaded the mRNA and lncRNA data of thyroid cancer patients from the TCGA ocial website https://portal.gdc.cancer.gov/, and downloaded the miRNA of thyroid cancer patients from UCSC Xena https://xenabrowser.net/ data. Through the tumor microenvironment's immune score, tumor samples are divided into cold tumors and hot tumors. We obtained the differential genes of DEGs, lncRNA and miRNA, and jointly constructed a ceRNA network through differential analysis of the mRNA data of cold and hot tumours. Results: We rst assessed the immune inltration of the tumor microenvironment of patients with thyroid cancer, and divided the samples into cold tumors and hot tumors based on the immune score. A total of 568 up-regulated DEGs and 412 down-regulated DEGs were screened by analysing the difference between hot and cold tumours. Then we conducted a differential analysis of lncRNA and miRNA and obtained 629 differential genes of lncRNA, and 114 differential genes of miRNA. Finally, we constructed a ceRNA network for the acquired differential genes. Five miRNA-centric hubnets, hsa-mir-204, hsa-mir-128, hsa-mir-214, hsa-mir-150 and hsa-mir-338, have been identied. Conclusion: Our research has identied the immune-related mRNA, lncRNA and miRNA DEGs of thyroid cancer, and constructed a ceRNA network. This result is of great signicance for exploring immune-related mechanisms of thyroid cancer and provides certain help for immunotherapy of thyroid cancer.


Background
In 1863, Rudolf Virchow, the father of modern cytopathology, proposed the relationship between microin ammation and subsequent cancer development [8].In 1909, Paul Ehrlich proposed using the immune system to control cancer [9]. Tumor microenvironment, especially the immune system, plays an essential role in regulating tumor progression and tumor response to treatment. It mainly stimulates tumor speci c immune response by inducing the immunogenic death of tumor cells or participating in immune response mechanism [10]. According to the immune in ltration of tumor microenvironment, tumor is divided into "cold tumor" and "hot tumor". In the tumor microenvironment of "hot" tumor, the degree of immune invasion and immune effect is relatively higher, with strong antigen presentation ability and T cell activation. Such immune in ltration will lead to the production of tumor speci c CD8 + T cells, which can clear cancer cells and generate systemic tumor speci c immunity, forming a long-term antitumor memory response [11,12].While "cold" tumors have no immune cell in ltration in the tumor microenvironment or are mainly in ltrated by suppressive regulatory cell subtypes (including regulatory T cells (Tregs), regulatory B cells (bregs) and myeloid suppressor cells (MDSCs)) [13][14][15]. The growth of cancer is not inhibited in immunology. In thyroid cancer, we divided tumor samples into cold tumor group and hot tumor group by the degree of immune invasion, and further screened the differential genes of lncRNA, miRNA and mRNA between the two groups to screen out the immune mechanism related to thyroid cancer.
Cancer Genome Atlas (TCGA) database can be used for large-scale global gene expression pro le analysis and database mining to nd the potential correlation between genes and the overall survival rate of various malignant tumors [16].In this study, we rst from TCGA o cial website https://portal.gdc.cancer.gov/ The fpkm data of mRNA and lncRNA of thyroid cancer patients were downloaded and transformed into TPM from the o cial website of ucscxena https://xenabrowser.net/ . The clinical information corresponding to the miRNA data set of patients with thyroid cancer was downloaded. We used ssgsea, MCP counter, cibersort, xcell four packages to evaluate the immune cells of thyroid cancer tumor samples and normal samples. We used xcell's own algorithm and estimate package to evaluate the immune and matrix scores. In addition, thyroid cancer was divided into four subtypes by congruent clustering to understand the differences in immune cell types, immune related molecules, tumor size distribution and grading among the four subtypes. According to the difference of immune score and matrix score, cluster 3 and Cluster 4 were de ned as a hot tumor, and cluster 1 and cluster 2 were de ned as cold tumors. Then, differential analysis of cold and hot tumors, enrichment analysis of differential genes, and protein interaction network construction were carried out.
Similarly, we compared the expression of lncRNA and miRNA in cold and hot tumors. Finally, the differentially expressed lncRNA, miRNA and mRNA obtained by us were used to construct the immune related Cerna network. Finally, we identi ed ve immune connected CERN networks in thyroid cancer, which is of great signi cance to understand the mechanism of immune invasion and immunotherapy.
The purpose of this study was to construct a Cerna regulatory network using microarray data collected from a public database, and preliminarily identify the regulatory mechanism mediated by a novel lncRNA miRNA mRNA Cerna in thyroid cancer. This study may provide a target for the development of new therapeutic strategies for thyroid cancer.

Data download and immune in ltration assessment
We are from TCGA o cial website https://portal.gdc.cancer.gov/ The mRNA and lncRNA data of thyroid cancer patients were downloaded from ucscxena database https://xenabrowser.net/ To download miRNA data from. Among them, there were 57 normal samples and 511 tumor samples. In tumor microenvironment, there are four commonly used methods to evaluate immune cell in ltration: single sample gene set enrichment analysis (ssgsea), microenvironment cell populations (MCP) -counter, cibersort and xcell [17][18][19][20]. In order to reduce the error, we put four methods into the study. Inclusion criteria: P < 0.05 (the p value of cibersort was p < 0.1 to obtain su cient samples).

Immune and matrix evaluation of tumor and paracancerous samples
In the tumor microenvironment, the immune and matrix scores were calculated by xcell's own calculation method and estimate package, Immune and stromal cells are two major types of non tumor components and have been proposed to be valuable for diagnostic and diagnostic assessment of tumors. There were signi cant differences in the distribution of immune cells and immune scores between tumor tissues and adjacent tissues.

Correlation of immune cells
We used the results of four evaluation methods to calculate the correlation between different immune cells. Among them, blue represents positive correlation and red represents negative correlation. Through the correlation analysis of thyroid cancer cells and adjacent immune cells, we know that the interaction between adjacent and cancer cells is signi cantly different.

Cell consistent clustering
Using the method of cell consistent clustering, we divided thyroid cancer samples into different subtypes according to the results of immune in ltration and drew the cumulative distribution function (CDF) to identify the optimal number of subsets. Finally, we identi ed four different subtypes and attracted the heat map to compare the immune matrix score and the distribution of immune cells among different subtypes.
Differences of immune matrix score, immune related molecules and tumor size distribution and grading among different subtypes The progress and metastasis of cancer depend on the two-way interaction between cancer cells and their environment, forming the tumor microenvironment (TME) [21].Tumor microenvironment is different in different tumour progression stages, which can promote tumor formation and inhibit tumor formation.
We have known that immune cells can be activated to promote tumor formation and progression [22]. The immune matrix score can predict the immune invasion of tumor microenvironment. We compared the immune matrix score, immune-related molecules' expression, tumor size distribution, and grade of four subtypes of immune cell in ltration.

Analysis of the difference between cold and hot tumors
Cancer progression requires tumor cells to produce immune tolerance. Although the immune tolerance mechanism is involved, the tumor can be divided into two subtypes according to T lymphocyte in ltration [23,24]. The different proportion of effector T cells and regulatory T cells in thermal tumors re ects the different degree of immunosuppression and affects the progress of tumor [25]. We compared the results of immune matrix score among different subtypes. Cluster 3 and 4 with higher immune matrix scores were classi ed as hot tumors, and cluster 1 and 2 were cold tumors. The lists of differentially expressed genes (DEGs) between controls and patients with DTC were generated using the LIMMA method[26], where statistical signi cance was set as |log-fold change (logFC)|>1 and Benjamini and Hochberg-corrected false discovery rates (FDR)<0.05. A hierarchical cluster heatmap based on Euclidean distance was generated using the pheatmap R package (Version: 1.0.12) and represents DEGs' expression intensity and direction.

Enrichment analysis of differentially expressed genes and construction of protein interaction network of differentially expressed genes
We use R software package "clusterpro le" (version 3.7.0; http://cytoscape.org/ Gene Ontology (go) and the Kyoto Encyclopedia of genes and genomes (KEGG) pathway enrichment analysis were performed.
Objective to search for tumor related regulatory molecules and enrich them Based on the above research, we obtained the differential genes between tumor and adjacent tumor, and then crossed the genes up-regulated in cold tumor (ADJP < 0.05) with the genes over expressed in tumor. These genes are related to cancer and participate in negative immune regulation. Similarly, the upregulated genes in hot tumors are crossed with the low expressed genes in cancer, which are positive immune regulatory genes. Finally, we used go and KEGG enrichment analysis to obtain the immune regulatory molecules.
Differences between lncRNA and miRNA in cold and hot tumors We performed the differential analysis (up: | log2fc | > 1, down: | log2fc | > 1, adjusted p-value < 0.05) by comparing cold and hot tumors in the R computing environment using limma package.

Data download, evaluation of immune invasion and comparison of immune cells between cancer and adjacent tumor
In this study, we rst from TCGA o cial website https://portal.gdc.cancer.gov/ We downloaded the mRNA and lncRNA data of osteosarcoma patients, miRNA and corresponding clinical data are UCSC xenahttps://xenabrowser.net/ It is downloaded directly from the database. We used four methods: single sample gene set environment analysis (ssgsea), microenvironment cell populations (MCP) -counter, cibersort and xcell to compare the differences of immune cells between cancer and adjacent tissues. The four methods included different numbers of immune cells (cibersort: 22, ssgsea: 28, MCP counter: 10, xcell: 67). In order to reduce possible errors, we included all four evaluation methods into the study. We used four methods to compare the difference of immune cells between cancer and paracancerous. The results showed that the number of most immune cells in cancer was signi cantly lower than that in paracancerous ( Figure 1A-D), and the immune and matrix scores in cancer were lower than that in paracancerous ( Figure 1E).

Correlation of immune cells
We analyzed the correlation of immune cells in tumor and paracancerous samples (Figure 2A-H). The results showed that the correlation of immune cells in tumor and paracancerous samples was signi cantly different, which may be due to the synergistic in ltration of immune cells activated by cancer antigen. In tumor tissue, the synergistic effect of different immune cells constitutes the tumor immune microenvironment, and plays an important role in tumor invasion and development.

Cell consensus clustering
We classi ed thyroid cancer samples into different subtypes according to the results of different immune in ltration. The cumulative distribution function (CDF) is drawn to identify the number of optimal subgroups ( Figure 3A, B & C). Finally, we identi ed four different subtypes. We compared the distribution of cells and the expression of immune matrix score between different subtypes by mapping heat map. Among the four subtypes, we compared the differences of tumor purity, immune score and stromal score in innate and adaptive immune mechanisms. We observed that the tumorpurity of cluster3 and 4 was lower than that of the other two subtypes, while the immunescore and stromalscore were higher than those of the other two subtypes ( Figure 3D).
Immune matrix scores and immune-related molecules among different subtypes Immune in ltration in the tumor microenvironment can be assessed by the immune matrix score. We respectively compared the immune matrix scores and the expression levels of immune-related molecules in 4 different immune in ltrating subtypes ( Figure 4A-F). According to the ESTINMATE algorithm, the matrix score is between -600 and 1700, and the immune score is between -800 and 3500. Among them, the immune scores of clusters 3 and 4 are signi cantly higher than those of clusters 1 and 2. This means that immune and matrix scores are meaningful in the classi cation of subtypes.

Differences in tumor size distribution and grade
The American Joint Committee on Cancer/International Cancer Control (AJCC/UICC) released the 8th edition of the TNM staging system [36].Where T represents the size of the cancer, N represents the lymph node metastasis or not, and M represents whether there is distant metastasis. However, the TNM staging system alone is not enough to accurately re ect the disease stage of the tumor, and other tumor-related factors have also been paid more and more attention, so the different grades of tumors (stage1-4) are named [37].We draw a pie chart to understand the difference in tumor size distribution and grading between different subtypes of thyroid cancer ( Figure 5A-H). Among them, the ratio of cluster4 at T1 and stage1 is signi cantly higher than that of other clusters, which suggests that hot tumors have a lower pathological stage than cold tumors.

Analysis of differences between cold and hot tumors
We de ne clusters 3 and 4 with higher immune matrix scores as hot tumors, and clusters 1 and 2 with lower immune matrix scores as cold tumors. We drew a heat map and found that the immune cell types of thyroid cancer and adjacent samples were signi cantly different ( Figure 6A). We used limma package (|Log2FC| > 2, adjusted P value < 0.05) to identify differentially expressed genes, and drew a volcano map to compare the molecular differences in gene expression between cold and hot tumors ( Figure 6B). Enrichment analysis of differential genes and construction of protein interaction network of differential genes Through screening, we found 568 up-regulated differential genes and 412 down-regulated differential genes. We use the clusterPro ler package to perform GO enrichment and KEGG pathway enrichment analysis on differential genes related to cold and hot tumors [38]. The up-regulated differential genes are mainly enriched in regulation of leukocyte activation (GO:0002694), T cell activation (GO:0042110), regulation of lymphocyte activation (GO:0051249) and leukocyte cell-cell adhesion (GO:0007159); downregulated The genes are mainly enriched in hormone metabolic process (GO:0042445), thyroid hormone metabolic process (GO:0042403), hormone biosynthetic process (GO:0042446) and thyroid hormone generation (GO:0006590) ( Figure 7A-D). KEGG pathway enrichment analysis showed that the upregulated differential genes were mainly enriched in Hematopoietic cell lineage (hsa04640), Cytokinecytokine receptor interaction (hsa04060) and Viral protein interaction with cytokine and cytokine receptor (hsa04061), and the down-regulated differential genes were mainly enriched in Thyroid hormone synthesis (hsa04918), Rap1 signaling pathway (hsa04015) and Cortisol synthesis and secretion (hsa04927) (Figure 7E-H). We used the STRING online tool to construct a protein interaction network (PPI) of differential genes related to cold and hot tumors. Among the up-regulated differential genes, we chose the node connection greater than 20 as the hub gene, and in the down-regulated differential gene we chose the node connection Those greater than 10 are de ned as hub genes ( Figure 7I: down-regulated differential gene protein interaction network; Figure 7J: up-regulated differential gene protein interaction network). We have 19 up-regulated hub genes and 9 down-regulated hub genes. In the protein interaction network, we have observed a wide range of links between the markers related to the matrix in the differential gene and the markers related to the immune in the differential gene, which may be related to the higher scores of immunity and matrix in hot tumors .
Looking for tumor-related regulatory molecules and enrichment analysis of tumor-related molecules T cell-based cancer immunotherapy, such as checkpoint suppression or adoptive cell therapy, has greatly changed the way cancer is treated [39].We know that immunity plays a very important role in the occurrence and development of tumors. We rst calculated the differential genes between cancer and adjacent cancers, and then intersected the genes up-regulated in cold tumors (adjpvalue<0.05) with the highly expressed genes in the tumors. A total of 717 genes were negatively regulated immune genes related to cancer. So these genes may be therapeutic targets. Similarly, the genes that are up-regulated in hot tumors (adjpvalue<0.05) are crossed with genes that are low-expressed in cancer. There are 1246 genes in total. These genes are positively regulated immunity, so agonists need to be added during the treatment. We enrich these differential genes. In GO enrichment analysis, the up-regulated differential genes are mainly enriched in: T cell activation (GO:0042110), regulation of T cell activation (GO:0050863), regulation of lymphocyte activation (GO:0051249), external side of plasma membrane (GO:0009897), cytokine activity (GO:0005125), etc. ( Figure 8A&B). In the enrichment analysis of the KEGG pathway, the up-regulated differential genes were mainly enriched in: Cytokine-cytokine receptor interaction (hsa04060) and Viral protein interaction with cytokine and cytokine receptor (hsa04061) ( Figure 8C&D). The down-regulated differential genes are not enriched in related pathways.

Differences between lncRNA and miRNA in cold and hot tumors
We will compare the expression differences of the downloaded lncRNA and miRNA between the two groups of cold and hot tumors ( Figure 9A-D). We used limma package to identify differentially expressed genes, and drawn volcano maps to compare the differentially expressed genes between cold and hot tumors. A total of 629 differential lncRNAs (|Log2FC| > 1, adjusted P value < 0.05) and 114 differential miRNAs (|Log2FC| > 0.5, adjusted P value < 0.05) were obtained.
To better understand the effect of lncRNAs on mRNAs mediated by combination with miRNAs in DTC, we built a ceRNA network based on the above mentioned data and used the ggalluvial R package (Version: 0.9.1) to visualize the network. miRNA target gene prediction passed miRDB, Targetscan, and miRTarBase take the intersection. The circle represents lncRNA (adj<0.05), the diamond represents mRNA (adj<0.05, logFC>1, logFC<0.05), and the triangle represents miRNA (adj<0.05). Finally, the cytohub module calculates the rst 5 networks de ned as hubnet ( Figure 10). Five miRNAs, hsa-mir-204, hsa-mir-128, hsamir-214, hsa-mir-150 and hsa-mir-338, are located in the central area of the network, and are of great signi cance to the immunity of thyroid cancer.

Discussion
Thyroid cancer is one of the malignant tumors, the screening of RNAs transcripts was facilitated in the past 20 years, lncRNAs and miRNAs were tested to be an emerging factor which strongly associated with tumorigenesis and metastasis in thyroid carcinoma [40][41][42] .Interactions between tumor cells and various components of TME are signi cant and contribute to all the hallmarks of cancer [43]. TME can affect how a tumor grows and spreads. Therefore, identi cation of key genes in thyroid cancer microenvironment aids appropriate management and treatment of thyroid cancer. The immune in ltration of TME is very important for the immune-related treatment of thyroid cancer. Our research aims to identify immunerelated mRNAs, lncRNAs and miRNAs, and further explore the relationship between them.
More and more studies indicate that lncRNAs play a vital role in biological functions through multiple levels of regulation, which involve transcriptional, posttranscriptional, and epigenetic regulation [44,45].A large number of studies have shown that there is a complex and closely related regulatory network between miRNA and lncRNA. In triple negative breast cancer [46], the relationship between miRNA and lncRNA has been determined. The ceRNA hypothesis was proposed to explain the mechanism of tumorigenesis; this hypothesis provides a novel guiding theory and suggests valuable strategies and research directions for the diagnosis and treatment of malignancies [47], lncRNAs with sequences similar to their target miRNA can regulate mRNA expression by acting as a sponge of miRNA.
At present, there are few studies on ceRNA networks related to thyroid cancer. In this study, we downloaded mRNA, lncRNA and miRNA of patients with thyroid cancer from TCGA. We rst divide the sample into cold tumors and hot tumors through immune evaluation [48].In the tumor microenvironment of "hot" tumors, the degree of immune in ltration and immune effect is relatively higher, and it has strong antigen presentation ability and T cell activation. Such immune in ltration will lead to the production of tumor-speci c CD8 + T cells, which can eliminate cancer cells and generate systemic tumor-speci c immunity, forming a long-term anti-tumor memory response [11,12]. However, "cold" tumors have no immune cell in ltration in the tumor microenvironment or are mainly in ltrated by inhibitory regulatory cell subtypes (including regulatory T cells (Tregs), regulatory B cells (Bregs) and myeloid suppressor cells (MDSCs)) [13][14][15], the growth of cancer is immunologically uninhibited. We used DEGs to construct ceRNA, and nally identi ed a hubnet centered on ve miRNAs: hsa-mir-204, hsa-mir-128, hsa-mir-214, hsa-mir-150 and hsa-mir-338. miRNAs are approximately 22 nucleotides long RNA molecules that bind to the 3'-untranslated region (3'-UTR) of their respective target genes, and exert their in uence on gene expression by inhibiting protein translation or degrading mRNA [49].We already know that miRNA plays a very important role in the occurrence and development of cancer. hsa-mir-204 has been determined to be related to many cancer malignancies, including: melanoma [50] , breast cancer [51] and liver cancer [52].hsa-mir-128-3p can increase the chemical sensitivity of colorectal cancer to chemotherapy drugs [53].hsa-mir-128-3p has an important relationship with the drug resistance of chemotherapeutics, but there are few related reports in thyroid cancer. It is worth noting that hsa-mir-214 has been found to play a very important role in regulating the proliferation and metastasis of papillary thyroid cancer cells [54].hsamir-150 and hsa-mir-338 are also related to tumor proliferation and cell invasion, such as colorectal cancer [55],non-small cell carcinoma [56] and cervical cancer [57].In thyroid cancer, these related hub miRNAs and related lncRNA and mRNA jointly affect the tumor microenvironment. It may be closely related to immune cell in ltration and tumor invasion. Because the number of deaths from thyroid cancer samples is relatively small, in this study, hubnet was not analyzed for survival.
Inevitably, our research had some innate limitations which need to be addressed. although this bioinformatics study was designed well, the major drawback of our study was the lack of in vivo and vitro validation. Our research has identi ed 5 hubnets related to thyroid cancer, which may be related to cancer immune in ltration and tumor invasion. This result has important signi cance for understanding the tumor invasion and mechanism of action in patients with thyroid cancer, and provides possible markers for treatment.

Conclusions
In this study, we identi ed the hubnet centered on ve miRNAs: hsa-mir-204, hsa-mir-128, hsa-mir-214, hsa-mir-150 and hsa-mir-338. In the ceRNA network we established, the relationship between miRNA, lncRNA and immune-related mRNA was determined. The understanding of the molecular role of thyroid cancer tumor microenvironment is of great signi cance. LGQ performed the data analyses. LYS contributed signi cantly to process data and wrote the manuscript. All of the authors read and approved the nal manuscript.

Ethical Statement
The data of this study are from TCGA and UCSC, and do not involve animal experiments and human specimens, no ethics-related issues.  counter evaluation method shows that the number of CD8 T cells, monocytic lineage, myeloid dendritic cells, NK cells and T cells in cancer is less than that in adjacent cancer, and the number of neutrophils in cancer is higher than that in adjacent cancer. (C)ssgsea evaluation method shows that the number of most T cells and B cells in cancer is lower than that in adjacent cancer, while the number of mast cells, monocytic lineage and neutrophils in cancer is higher than that in adjacent cancer. (D)X-cell evaluation method shows that the number of B cells, T-cells and monocyte in cancer is less than that in adjacent cancer, while the number of dendritic cells, mast cells and natural killer cells is higher than that in adjacent cancer.(E)the immune and matrix scores in carcinoma were lower than those in adjacent carcinoma.     ratio of cluster4 in T1 and stage1 is signi cantly higher than that of other clusters.     Construction of ceRNA network. A ceRNA network consisting of 629 lncRNAs, 114 miRNAs and 980 mRNAs was constructed. We select the rst ve networks and de ne them as hubnet. The circle represents lncRNA (adj<0.05), the diamond represents mRNA (adj<0.05, logFC>1, logFC<0.05), and the triangle represents miRNA (adj<0.05).