- Open Access
Identification of ferroptosis genes in immune infiltration and prognosis in thyroid papillary carcinoma using network analysis
BMC Genomics volume 22, Article number: 576 (2021)
Papillary thyroid carcinoma (PTC) is the most common thyroid cancer. While many patients survive, a portion of PTC cases display high aggressiveness and even develop into refractory differentiated thyroid carcinoma. This may be alleviated by developing a novel model to predict the risk of recurrence. Ferroptosis is an iron-dependent form of regulated cell death (RCD) driven by lethal accumulation of lipid peroxides, is regulated by a set of genes and shows a variety of metabolic changes. To elucidate whether ferroptosis occurs in PTC, we analyse the gene expression profiles of the disease and established a new model for the correlation.
The thyroid carcinoma (THCA) datasets were downloaded from The Cancer Genome Atlas (TCGA), UCSC Xena and MisgDB, and included 502 tumour samples and 56 normal samples. A total of 60 ferroptosis related genes were summarised from MisgDB database. Gene set enrichment analysis (GSEA) and Gene set variation analysis (GSVA) were used to analyse pathways potentially involving PTC subtypes. Single sample GSEA (ssGSEA) algorithm was used to analyse the proportion of 28 types of immune cells in the tumour immune infiltration microenvironment in THCA and the hclust algorithm was used to conduct immune typing according to the proportion of immune cells. Spearman correlation analysis was performed on the ferroptosis gene expression and the correlation between immune infiltrating cells proportion. We established the WGCNA to identify genes modules that are highly correlated with the microenvironment of immune invasion. DEseq2 algorithm was further used for differential analysis of sequencing data to analyse the functions and pathways potentially involving hub genes. GO and KEGG enrichment analysis was performed using Clusterprofiler to explore the clinical efficacy of hub genes. Univariate Cox analysis was performed for hub genes combined with clinical prognostic data, and the results was included for lasso regression and constructed the risk regression model. ROC curve and survival curve were used for evaluating the model. Univariate Cox analysis and multivariate Cox analysis were performed in combination with the clinical data of THCA and the risk score value, the clinical efficacy of the model was further evaluated.
We identify two subtypes in PTC based on the expression of ferroptosis related genes, with the proportion of cluster 1 significantly higher than cluster 2 in ferroptosis signature genes that are positively associated. The mutations of Braf and Nras are detected as the major mutations of cluster 1 and 2, respectively. Subsequent analyses of TME immune cells infiltration indicated cluster 1 is remarkably richer than cluster 2. The risk score of THCA is in good performance evaluated by ROC curve and survival curve, in conjunction with univariate Cox analysis and multivariate Cox analysis results based on the clinical data shows that the risk score of the proposed model could be used as an independent prognostic indicator to predict the prognosis of patients with papillary thyroid cancer.
Our study finds seven crucial genes, including Ac008063.2, Apoe, Bcl3, Acap3, Alox5ap, Atxn2l and B2m, and regulation of apoptosis by parathyroid hormone-related proteins significantly associated with ferroptosis and immune cells in PTC, and we construct the risk score model which can be used as an independent prognostic index to predict the prognosis of patients with PTC.
Thyroid cancer (TC) is the most common endocrine malignant tumour worldwide and its incidence has been increasing . In 1990, there were an estimated 95,030 incident cases of TC and 22,070 deaths; this increased to 255,490 incident cases and 41,240 deaths in 2017 . In 2020, TC has been found in 9th place for incidence among all cancers, with more than 586,000 cases diagnosed worldwide. In addition, there is a 3-fold higher incidence rate of TC in women, which is about 10.1 per 100,000 cases, representing one in every 20 cancers diagnosed among women . TC shows the highest incidence rates in Northern America, Australia/New Zealand, Eastern Asia, and Southern Europe. The incidence rate of TC is also increasing in China, and it has become one of the ten major cancers threatening the health of Chinese residents . The rapid increase of thyroid cancer has been largely attributed to the use of progressively sensitive diagnostic imaging modalities [5, 6]. A study also found that obesity positively correlated with 16% of TC cases and 63% of large-size tumors diagnosed from 2013 to 2015 in the USA . Differentiated thyroid cancer, which originates from thyroid follicular epithelial cells, accounts for more than 95% of TC, among which papillary thyroid carcinoma (PTC) constitutes more than 85% of cases [8, 9].
Traditional treatment methods of PTC include radical surgery, endocrine therapy and 131I therapy, and most treatments display relatively good efficacy, e.g., patients who are disease-free after one course of 131I ablation show a range of in 1–3% recurrence rate . However, a small proportion of PTC displays high aggressiveness and even develops into RAI-refractory differentiated thyroid cancer (RAIR-DTC). This causes roughly 15% of the patients to experience recurrence within 10 years after the initial treatment [11,12,13]. Accurate assessment of the prognosis of PTC is essential to ensure that high-risk and advanced patients receive appropriate treatment without over-treating low-risk patients. The American Thyroid Association (ATA) currently recommends the use of TNM staging to predict mortality and proposes a system to assess the risk of recurrence, which includes the size of the main tumour and whether it has grown into nearby areas (T), the extent of spread to nearby lymph nodes (N), and whether the cancer has metastasised to other organs of the body (M) [9, 14]. It has been addressed that tumour microenvironment (TME) possibly plays an important role in response to chemotherapy and antiangiogenic therapy [15, 16]. In the process of tumorigenesis, tumour cells interact with the surrounding microenvironment to promote immune tolerance, which then develops into immune escape that eventually results in tumour development and angiogenesis, invasion, metastasis, and chronic inflammation .
TME has been found to be a complex and continuously evolving matrix, containing such cells as stromal cells, fibroblasts, endothelial, innate and adaptive immune cells . The transcriptome of PTC has been well-characterised and the major molecular events associated with most PTC cases, including the aberrant expression of Brafv600e, Ras, point mutations of Tert and RET/PTC rearrangements were elucidated [18,19,20]. However, reports on the interaction between PTC and stromal tissue, lymphocytic infiltrate, or normal thyrocytes and expression profiles of associated marker genes are limited. A study of 23 PTC patients using oligonucleotide microarrays identified 19 genes differentially expressed with reference to 10 patients with other thyroid disease, such as Lrp4, Eva1, Tmprss4, Qpct, and Slc34a2 . A recent study comparing the gene expression profiles of PTC and normal thyroid in both micro-dissected cells and whole tissue slides revealed crosstalk between cancer cells and TME possibly involving the functions of Ptcsc, Ctgf, Tff3, Fn1, Mpped2 .
As a type of programmed cell death, ferroptosis is dependent on iron and induced by the accumulation of oxidatively damaged phospholipids, associated with the malfunction of glutathione-dependent antioxidant defences that are mediated by glutathione peroxidase 4 (Gpx4) via different pathways [23,24,25]. and, resulting in a large amount of ROS, which promotes ferroptosis The declines in metabolism of lipid peroxides catalysed by Gpx4 and glutathione (GSH) level intracellularly, lead to Fe2+ oxidising lipids in a Fenton-like manner, which enhances ferroptosis due to the elevation of lipid reactive oxygen species (ROS) in cells [24, 26]. In recent years, the induction of ferroptosis has been investigated as an alternative and/or joint therapeutic approach to trigger cancer cell death, with respect to other types of cell death, especially for the treatment of malignancies resistance issues in some cancers; the marker genes and modulator molecules of ferroptosis have been identified [13, 27, 28]. However, the role of ferroptosis in PTC currently remains elusive.
In this study, we integrated the THCA dataset from TCGA and the clinical databases of THCA from UCSC Xena to identify reliable differentially expressed genes (DEGs) relevant to ferroptosis in PTC. Then, univariate Cox survival analysis and lasso Cox regression analysis were performed to identify DEGs of PTC, and we proposed a prognostic prediction model using DEGs and clinical data from the TCGA-THCA and UCSC Xena datasets. We performed Multivariate Cox survival analysis on hub genes combined with clinical prognosis data and showed that the model can be used to predict PTC prognosis.
The papillary thyroid carcinoma (THCA) datasets were obtained from two platforms. THCA database level 3 count was searched in The Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov/repository). The datasets in TCGA were transformed into Transcripts Per Million (TPM) values and the mutation data of The Single Nucleotide Polymorphism Database (SNP) were downloaded for THCA. The clinical data of THCA were obtained using the UCSC Xena browser (https://xenabrowser.net/). Only patients with primary tumours who did not receive neoadjuvant therapy were included in this retrospective analysis . Their clinicopathological, genetic, epigenetic, and survival data were downloaded for a secondary analysis. In this study, the datasets of 502 tumour samples and 56 normal samples were used.
Identification of THCA subtypes based on related genes of ferroptosis
A set of 60 related genes of ferroptosis genes were summarised from the MisgDB database (https://www.gsea-msigdb.org/gsea/msigdb). The ConsensusClusterPlus package was used for consensus clustering and molecular subtype screening. In brief, k-means clustering was used, with 50 iterations (each using 80% of the samples) . The best cluster number was determined by the clustering score for the cumulative distribution function (CDF) curve, and the relative changes in the area under the CDF curve were evaluated. Principal component analysis (PCA) and t-distributed Stochastic Neighbor Embedding (t-SNE), commonly used for dimensionality reduction, were used to verify the reliability of the consensus clusters.
The ssGSEA score xi for each THCA sample i was converted to xi′ using the equation:
while xmax and xmin represent the maximum and minimum single-sample gene-set enrichment analysis (ssGSEA) scores for all samples in the THCA dataset, respectively. The relationship among subtypes, clinical parameters and ferroptosis genes gene expression were shown by thermography with R-package pheatmap.
Comparison of immune cell subgroups among molecular subtypes of THCA
We continued to use maftools to explore the difference between two subtypes and drew the waterfall of them. In order to investigate the biological difference and understand the different pathways involved by the subtypes, Gene Set Enrichment Analysis (GSEA) was used with the clusterProfiler package. Gene set variation analysis (GSVA) was used to further analyze the difference of pathways between the subtypes using GSVA package [30, 31].
Immune infiltration microenvironment in THCA
The difference between the subtypes may be due to the complexity of tumour microenvironment, so we used the ssGSEA algorithm to analyse the proportion of 28 kinds of immune cells in THCA, and used hclust algorithm to classify the immune cells according to the proportion of immune cells. Finally, a box-plot was constructed to contrast immune cell content between the subtypes. If the normal distribution and homogeneity of variance were met, T-test was used; if not, Wilcox test was used. The above operations were completed by using complexheatmap package, ssGSEA package, corrplot package and ggplot2 package.
Gene co-expression network analysis
To investigate characteristic immune gene distribution in each molecular subtype of THCA and identify genes or gene modules highly related to immune cell infiltration, the WGCNA R package was used to evaluate 502 immune-related genes comprising the expression matrix. WGCNA network was constructed and the significant modules were identified. Genes of the significant modules were selected as subtype related hub genes for follow-up analysis.
Data validation of the differentially expressed genes (DEGs)
To further analyse DEGs, DEseq2 algorithm was used for differential analysis of sequencing data (using adjust P-value < 0.05 and logFC > 2 as thresholds). Volcano plots and heatmap were applied to visualize the distribution of the overlapping DEGs between the training and test sets . The above operation was done using the DEseq2 package.
Functional enrichment, pathway analysis of hub genes
In order to identify the possible functions and pathway of hub genes, Gene Ontology (GO) analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway [33, 34] enrichment analyses were performed using the clusterProfiler package . The visualisation module of clusterProfiler was used for displaying analysis results. P < 0.05 was selected as the cut-off criterion.
Cox regression and survival analysis
To explore the clinical efficacy of hub genes, univariate Cox analysis was performed on hub genes combined with clinical prognosis data, and the results of univariate Cox regression (P < 0.05) were included in lasso regression. The lasso regression results were incorporated into multivariate Cox regression. The risk regression model was built according to the multivariate analysis results, and the formula of the model was as follows:
ROC curve and survival curve were used to evaluate the results of the model. Univariate Cox analysis and multivariate Cox analysis were performed by combining THCA clinical data and risk score values for further evaluation.
All the above analyses were completed by R software. Adjust P-value < 0.05, P-value < 0.05 and FC > 2 were used as statistical thresholds. The statistical methods and algorithms used are described in the corresponding steps.
Identification of THCA subtypes based on immune gene sets
The 502 tumour samples of THCA were divided into k subtypes (k = 2–9) using the R package ConsensusClusterPlus. Consensus distributions for each k were displayed, which dissect the optimal k value at which the sample distribution was stable (Fig. 1A). Based on the consensus score of the CDF curve, the optimal cluster number was determined as k = 2, and the relative changes in the area under the CDF curve were evaluated. These were finally divided into two distinct and nonoverlapping subtypes, i.e., cluster 1 and cluster 2. The consensus matrix heat map of these two clusters was shown in Fig. 1B. A principal component analysis (PCA) plot (Fig. 1C) and a two dimensional t-SNE analysis (Fig. 1D) further verified the LUAD cohort and the reliability of the consensus clusters.
Association among subtypes, clinical parameters and ferroptosis signature genes
To explore the association among subtypes, clinical parameters and ferroptosis signature genes, we analysed gene cluster, tumour stage, survival status, pathologic stage, pharmaceutical, radiation, race, gender, age, and ferroptosis signature genes in the TCGA cohort. The relationship among subtypes, clinical parameters and ferroptosis signature genes is derived by thermography (Fig. 2). On the left side of the heatmap, genes that were positively associated with ferroptosis were coded red (such as Csf2, Hsbp1, Pebp1, Nfs1, Hmgcr, Sqle, Mt1g, and Pgd), whereas blue colour was used to indicate genes that are negatively associated with ferroptosis (incl. Tp53, Keap1, Rpl8, Nox1, Fth1, Hmox1, Acsl4, Sat1, Gpx4, Ptgs2, Gls2, Fancd2, Alox5, Slc1a5 and Dpp4). Gene cluster differed significantly between the two subtypes. In particular, the proportion of cluster 2 was significantly higher than cluster 1 in ferroptosis signature genes that are positively associated, yet the expression of these genes generally appeared higher in cluster 2, such as Acsf2, Mt1g, Gclc and Aifm2. Moreover, cluster 1 had remarkably larger proportion of genes that are negatively associated with ferroptosis, for instance, Tp53, Dpp4, Slc1a5, Alox5 and Ptgs2.
In terms of the spread magnitude to nearby lymph nodes (stage_N), N1 was largely relevant to cluster 1, while N0 and NX were more present in cluster 2. The metastasis to distant sites (stage_M) showed an increased proportion of MX in cluster 2, while the extent of the tumour (stage_T) distributed across both cluster 1 and 2, with latter exhibited higher proportion of T2 (Fig. 2). There were more pathological stage III and IV cases associated with cluster 1, whereas cluster 2 subtype showed more genes associated with stage II, which represents non-invasive tumour with no spread to lymph nodes and no metastasis.
Somatic mutation landscape of the two subtypes were identified by GSEA
We continued to explore the differences between the two subtypes using the maftools package to explore the mutation profile and drawing the waterfall map of them (Fig. 3). It was found that the somatic mutation types of the two subtypes were missense mutations. The top 2 somatic mutation genes of subtype cluster 1 were Braf (88%) and Tnn (7%) (Fig. 3A), with the other 47 genes with low mutation frequency less than 5%. Cluster 2 showed somatic mutations of Nras, Hras and Tg with more than 10% of mutation frequency (Fig. 3B). This indicates biological differences between the two subtypes.
Gene set enrichment analysis
GSEA and GSVA were used to analyse the enriched pathways in each cluster to investigate the molecular differences in the phenotypes of these two subtypes. GSEA analysis suggested that subtype cluster 1 was mainly involved in bacterial invasion of epithelial cells, cell adhesion molecules, completion and coaggregation cascades and intestinal immune network for IgA production (Fig. 4A). In subtype cluster 2 phenylalanine, tyrosine and tryptophan biosynthesis, thyroid hormone synthesis, tyrosine metabolism, ubiquinone and other terpenoid−quinone biosynthesis were enriched (Fig. 4B). GSVA was used to further explore the differences in participating KEGG pathways between the two subtypes, and the results were shown in Fig. 4C (adjust P-value< 0.05; see Additional file 1: Table S1 for more details). In addition, there were more immune response relevant pathways present in cluster 1 than cluster 2, such as antigen processing and presentation (P-value = 1.27E-22), natural killer cell mediated cytotoxicity (P-value = 1.36E-37), Fc gamma R-mediated phagocytosis (P-value = 1.12E-23), cytokine-cytokine receptor interaction (P-value = 7.18E-45), and toll-like receptor signalling pathway (P-value = 1.17E-24). Cluster 2 was enriched with more metabolism pathways, for instance, metabolisms of glycine, serine and threonine (P-value = 1.04E-30), porphyrin and chlorophyll (P-value = 2.93E-33), and fatty acid (P-value = 4.05E-41).
Immune characteristics of the two subtypes
The difference between two immune subtypes might be formed because of the complexity of the TME. To explore the biological behaviours between the two clusters, we performed ssGSEA algorithm to analyse the proportion of 28 kinds of immune cells in immune infiltration microenvironment in THCA. The results showed that the degree of immune infiltration in cluster 1 was largely higher than in cluster 2, suggested by the greater activation of immune response relevant cells, among which activated B cells, CD8 and CD4 T cells appeared to be the top three cell types (Fig. 5A). Spearman correlation analysis was used to calculate the correlation between ferroptosis gene expression and the proportion of immune infiltrating cells (Fig. 5B). This showed positive correlation of Alox5, Fancd2, Dpp4, Hmox1, Ptgs2, Slc1a5, Sat1 and Acsl4 with most cell types with relatively high magnitudes. In contrast, the negative correlation of these cells with ferroptosis was indicated by the relative low expression of Acsf2, Pebp1, Nfs1, Hsbp1, Hmgcr, Acsl3 and Gss. The difference of immune cell contents between the two subtypes were compared (Fig. 5C and Additional file 2: Table S2). Notably, cluster1 was remarkably enriched (P-value < 0.0001) with nearly all innate immune cell types, except eosinophil and monocyte, compared to cluster 2.
Identification of DEGs distinct phenotypes derived from gene co-expression network
The distribution of genes expression profiles in the subtypes was investigated to identify genes or gene modules highly related to immune infiltration microenvironment. The WGCNA R-package was used to evaluate the expression matrix data from 502 tumour samples. Analysis of network topology for various soft-thresholding powers indicated that when the power value was equal to 5 (β = 5), the predicted gene co-expression network exhibited scale-free topology by the fit index greater than 0.8, with inherent modular features (Additional file 3: Fig. S1A). The adjacency function was used to generate the adjacency matrix based on the β and gene expression matrix. The hierarchical clustering was built based on the TOM dissimilarity measure (Fig. 6A). A total of 25 seven co-expression modules were detected. The module preservation statistics was employed to achieve reliable and preserved modules (Fig. S1B). The co-expression network was examined by the NULMS Stanford dataset, with genes assigned to modules based on the modules in the reference dataset. There are 16 modules strongly preserved (Z-summary more than 10), such as blue, green, turquoise yellow and red modules; while seven modules are moderately preserved (5 < Z-summary < 10).
The association between module eigengenes, subtypes and immunity were then computed through the Pearson’s correlation coefficient, which evaluates the P-value was calculated for any given correlation. The blue module was most significantly correlated with subtypes, and the yellow module had the highest correlation with immune typing (Additional file 3: Fig. S1C). Therefore, the genes of these two modules were selected as subtype related hub genes for subsequent analysis, which included a total of 1747 hub genes. The DEseq2 algorithm was used to further analyse the hub genes in the sequencing data, and significant differentially expression of 814 hub genes were found, among which 308 were up-regulated and 506 were down-regulated in tumour tissues (Fig. 6B). The heat map comparing the expression of these 814 genes in Log2FC values was displayed in Fig. 6C (See Additional file 4: Table S3 for the full gene list).
Gene ontology and KEGG pathway analysis showed functional enrichment in immune regulations
To further analyse the functions and pathways that 814 hub genes are potentially involved in, GO term and KEGG pathway enrichment analyses were performed using the clusterProfiler package. GO analysis showed that 814 hub genes were mainly involved in biological processes such as regulation of small GTPase mediated signal transduction, positive regulation of endocytosis and regulation of plasma lipoprotein particle levels. In terms of cellular components, these genes showed high relevance to membrane (Fig. 7A). KEGG pathway analysis showed that 814 hub genes were mainly enriched in Salmonella infection and endocytosis pathways (Fig. 7B). The enrichment analysis of Wikipathway revealed that regulation of apoptosis by parathyroid hormone-related protein, apoptosis modulation and signalling, and pathogenic Escherichia coli infection were the top 3 pathways (Fig. 7C). Besides, several cancer relevant pathways were also enriched in the hub genes (Additional file 4: Table S3).
Establishment and assessment of the predict model
To explore the clinical efficacy of 814 hub genes, Univariate Cox analysis was performed with clinical prognosis data and the results of P-value < 0.05 were recorded in Table S3. The results of Univariate Cox regression (P-value < 0.05) were included in lasso regression. Dimensionalisation was reduced according to the Lambda curve (Fig. 8A) and the proportional hazards model curve (Fig. 8B), which showed that the deviance was the smallest when the number of genes in the model was 19. The lasso regression results were incorporated into multivariate Cox regression, as seven genes were identified by LASSO regression (Additional file 4: Table S3), including Ac008063.2, Apoe, Bcl3, Acap3, Alox5ap, Atxn2l and B2m. The risk regression model diagram was displayed in Fig. 8C according to the multivariate analysis results. In order to evaluate the specificity and sensitivity of the prognostic model, ROC curve and survival curve were used for evaluation. The AUC value of the ROC curve was 0.748 (Fig. 8D) and the survival curve showed a significant difference between the high-risk and low-risk groups (P-value = 6.329E-4), indicating a good result of the model (Fig. 8E).
The risk score predicted prognosis
To further evaluate the clinical efficacy of the model, univariate Cox analysis (Fig. 9A) and multivariate Cox analysis (Fig. 9B) were performed by combining THCA clinical data and risk score values, which implied that the risk score value of the model could be used as an independent prognostic indicator to predict the prognosis of patients with thyroid papillary carcinoma (see Additional file 5: Table S4 for full statistical results).
The incidence of PTC is a heavy disease burden around the world, and has increased quite significantly over the past two decades [36,37,38]. Patients with recurrent PTC who undergo surgical treatment have a high risk of complications, including tracheal collapse, laryngeal edema and recurrent laryngeal nerve injury [9, 36, 39]. Furthermore, PTC is a highly heterogeneous disease and tumour progression involves a complex network comprised of multiple signalling pathways [40, 41]. Nomogram is widely used to evaluate clinical prognosis in oncology because it can integrate a variety of prognostic determinants, including molecular biological and clinicopathological parameters [8, 40, 42]. The combination of our prognostic gene signature with conventional clinical parameters may provide a better tool to predict PTC prognosis.
The discovery of distinct lethal subroutines of RCD has facilitated great progress in cancer treatment [43,44,45]. Ferroptosis was identified as a new form of RCD distinct from apoptosis, necrosis, autophagy, and other forms of cell death. It is an iron-dependent cell death process and caused by lipid peroxidation relying on reactive oxygen species (ROS) generation [24, 46, 47]. There has been growing interest in exploring the mechanisms underpinning ferroptosis in recent years, and several seminal discoveries have elucidated the process [44, 48, 49]. Ferroptosis is thought to only be inhibited by the phospholipid hydroperoxide-reducing enzyme glutathione peroxidase 4 and radical-trapping antioxidants [26, 44]. GSH-GPX4 and FSP1-CoQ10 are the two main ferroptosis resistance pathways [50, 51].
Growing evidence suggests that immune cells in TME play vital roles in tumorigenesis. These innate immune cells include macrophages, neutrophils, dendritic cells, innate lymphoid cells, myeloid-derived suppressor cells, and natural killer cells which potentially possess tumour-antagonising or tumour-promoting functions [29, 52, 53]. Cancer cells exhibited iron ion aggregation in the TME for active proliferation, thus from the perspective of iron homeostasis, regulating ferroptosis could effectively kill tumour cells. In addition, ferroptosis also plays an important immunological role, participating in the process of tumour surveillance by immune cells and tumour immunity [52, 54, 55]. The metabolism of immune cells affects their differentiation and function . Given the complex interaction of environmental factors in TME as a profound impact on the metabolic activities of immunity, matrix and tumour cell types and ferroptosis [48, 57]. Thus, exploring immunophenotype in TME could help reveal the significance of ferroptosis in cancer treatment. The cross study of ferroptosis and immune cells might lead to a discussion of clinical application.
In this study, we seek to investigate whether there was ferroptosis regulation via exploring the crucial genes and pathways of ferroptosis in PTC. A total of 502 PTC samples and 56 normal samples integrated from TCGA and UCSC Xena datasets were analysed, and 60 ferroptosis related genes were derived from MisgDB datasets. Two distinct subtypes (i.e. cluster 1, cluster 2) were identified in our study accordingly based on these related genes. The proportion of cluster 1 was significantly higher than cluster 2 in the signature genes that either positively or negatively associated with ferroptosis. The main mutation types of these two subtypes were identified as missense mutation. The missense mutations of Braf have been detected in circulating DNA in the serum of some patients with PTC . A close association of Braf mutations with extrathyroidal extension, lymph node metastasis, and stage III/IV of PTC has been strongly suggested . The significance of missense mutations of Nras and Hras in PTC has also been identified [60,61,62].
The distinct enriched KEGG pathways revealed GSEA analysis of cluster 1 and cluster 2 implied biological differences between the two genotypes. The GSVA analysis showed that there were more immune response relevant pathways enriched in cluster 1 (adjust P-value < 0.05, Additional file 1: Table S1). The difference between the two subtypes possibly resulted from the complexity of TME, since the degree of immune infiltration in cluster 1 was remarkably higher than that in cluster 2, which was reflected by the contents of many immune cell types. Subsequent analysis of TME cell infiltration indicated THCA cluster1 was remarkably rich in innate immune cell infiltration including T cells, natural killer cells, macrophage, eosinophil, mast cell, MDSC, plasmacytoid dendritic cell.
We subsequently analysed the related hub genes of the two modules with the highest correlation between molecular typing and immune typing in WGCNA network, and there were 814 DEGs including 308 up-regulated genes and 506 down-regulated genes in PTC tissues. Then, lasso regression was evaluated and incorporated into multivariate Cox regression, the crucial genes (i.e., Ac008063.2, Apoe, Bcl3, Acap3, Alox5ap, Atxn2l and B2m) were identified and the risk regression model was constructed according to the results of multivariate analysis. ROC curve and survival curve indicated a good result of the model. Based on these findings, the risk score of the model could be used as an independent prognostic indicator to predict the prognosis of patients with PTC.
All of these genes, with the exception of Ac008063.2, have been implicated in cancer cell proliferation. With the exception of Apeo, which is highly expressed in liver cancer, all of the genes have low cancer specificity, according to Human Protein Atlas (Ac008063.2 could not be identified). B2m has been identified as a potential biomarker for thyroid cancer, kidney disfunction and renal disease . Higher rates of B2m mutation are correlated with lower patient survival rates, which is speculated to be due to increased immune evasion [64, 65]. Atxn2l is a stress response molecule which has been identified as a potential cancer prognosis gene for cancer, such as adrenocortical carcinoma [66, 67]. Alox5ap has been identified as specifically upregulated in leukemia stem cells and unchanged in hematopoietic stem cells, indicating specificity in cancer cell self-renewal and differentiation [68, 69]. There is relatively little oncological literature on Acap3, however it has been identified as highly down-regulated in late-stage liver cancer patients . Bcl3 up-regulation has been associated with poorer prognostic outcomes [71, 72]. It is classified as an oncogene and has been implicated in breast cancer cell migration [73, 74]. Apoe is associated with tumorigenesis in many cancers, including lung, gastric and thyroid cancer, and a higher risk of metastasis [75,76,77].
To our knowledge, a prognostic model based on ferroptosis genes and immune cells in tumour immune infiltration microenvironment has not been reported up to now. Compared with whole genome sequencing, our prediction model based on a limited number of gene expression levels might be more economical and practical. Further, the risk score of the model could predict the prognosis of patients with PTC. However, the external validation of the seven-gene signature and prognostic nomogram is needed in more independent cohorts. In addition, to clarify the molecular mechanism underlying the seven genes in relation to PTC, the spatial expression pattern and quantity of these genes at protein level warrants further investigation.
We demonstrated that ferroptosis was associated with immune cell infiltration in TME of patients with PTC, based on which a seven-gene signature and a prognostic model to predict the prognosis was established. The seven genes appeared closely related to the progression and prognosis of PTC and thus could be potential therapeutic targets. The predicted model proved to be reliable in predicting the prognosis of patients with PTC and might thus be beneficial for individualised treatment and medical decision making.
Availability of data and materials
The datasets generated and/or analysed during the current study are available in The Cancer Genome Atlas database (https://portal.gdc.cancer.gov/projects/TCGA-THCA, accession number phs000178) and the UCSC Xena (https://xenabrowser.net/datapages/?cohort=TCGA%20Thyroid%20Cancer%20(THCA)&removeHub=https%3A%2F%2Fxena.treehouse.gi.ucsc.edu%3A443).
Papillary thyroid carcinoma
The Cancer Genome Atlas
Gene set enrichment analysis
Gene set variation analysis
Kyoto Encyclopedia of Genes and Genomes
Weighted correlation network analysis
Receiver operating characteristic
Area under the ROC Curve
Siegel RL, Miller KD, Jemal A. Cancer statistics, 2020. CA Cancer J Clin. 2020;70(1):7–30.
Deng Y, Li H, Wang M, Li N, Tian T, Wu Y, et al. Global burden of thyroid Cancer from 1990 to 2017. JAMA Netw Open. 2020;3(6):e208759. https://doi.org/10.1001/jamanetworkopen.2020.8759.
Ferlay J, Ervik M, Lam F, Colombet M, Mery L, Piñeros M, et al. Global cancer observatory: cancer today. Lyon: International Agency for Research on Cancer; 2020.
Du L, Zhao Z, Zheng R, Li H, Zhang S, Li R, et al. Epidemiology of thyroid Cancer: incidence and mortality in China, 2015. Front Oncol. 2020;10:1702. https://doi.org/10.3389/fonc.2020.01702.
Udelsman R, Zhang Y. The epidemic of thyroid cancer in the United States: the role of endocrinologists and ultrasounds. Thyroid. 2014;24(3):472–9. https://doi.org/10.1089/thy.2013.0257.
Brito JP, Morris JC, Montori VM. Thyroid cancer: zealous imaging has increased detection and treatment of low risk tumours. Bmj. 2013;347(aug27 4):f4706. https://doi.org/10.1136/bmj.f4706.
Kitahara CM, Pfeiffer RM, Sosa JA, Shiels MS. Impact of overweight and obesity on US papillary thyroid Cancer incidence trends (1995-2015). J Natl Cancer Inst. 2020;112(8):810–7. https://doi.org/10.1093/jnci/djz202.
Wu M, Yuan H, Li X, Liao Q, Liu Z. Identification of a five-gene signature and establishment of a prognostic nomogram to predict progression-free interval of papillary thyroid carcinoma. Front Endocrinol (Lausanne). 2019;10:790. https://doi.org/10.3389/fendo.2019.00790.
Carling T, Udelsman R. Thyroid cancer. Annu Rev Med. 2014;65(1):125–37. https://doi.org/10.1146/annurev-med-061512-105739.
Verburg FA, Van Santen HM, Luster M. Pediatric papillary thyroid cancer: current management challenges. Onco Targets Ther. 2017;10:165–75. https://doi.org/10.2147/OTT.S100512.
Rivera M, Ghossein RA, Schoder H, Gomez D, Larson SM, Tuttle RM. Histopathologic characterization of radioactive iodine-refractory fluorodeoxyglucose-positron emission tomography-positive thyroid carcinoma. Cancer. 2008;113(1):48–56. https://doi.org/10.1002/cncr.23515.
Riesco-Eizaguirre G, Galofre JC, Grande E, Zafon Llopis C, Ramon y Cajal Asensio T, Navarro Gonzalez E, et al. Spanish consensus for the management of patients with advanced radioactive iodine refractory differentiated thyroid cancer. Endocrinol Nutr. 2016;63(4):e17–24.
Luo J, Zhang B, Cui L, Liu T, Gu Y. FMO1 gene expression independently predicts favorable recurrence-free survival of classical papillary thyroid cancer. Future Oncol. 2019;15(12):1303–11. https://doi.org/10.2217/fon-2018-0885.
Haugen BR, Alexander EK, Bible KC, Doherty GM, Mandel SJ, Nikiforov YE, et al. 2015 American Thyroid Association management guidelines for adult patients with thyroid nodules and differentiated thyroid Cancer: the American Thyroid Association guidelines task force on thyroid nodules and differentiated thyroid Cancer. Thyroid. 2016;26(1):1–133. https://doi.org/10.1089/thy.2015.0020.
Roma-Rodrigues C, Mendes R, Baptista PV, Fernandes AR. Targeting tumor microenvironment for cancer therapy. Int J Mol Sci. 2019;20(4):840. https://doi.org/10.3390/ijms20040840.
De Palma M, Biziato D, Petrova TV. Microenvironmental regulation of tumour angiogenesis. Nat Rev Cancer. 2017;17(8):457–74. https://doi.org/10.1038/nrc.2017.51.
Denton AE, Roberts EW, Fearon DT. Stromal cells in the tumor microenvironment. Adv Exp Med Biol. 2018;1060:99–114. https://doi.org/10.1007/978-3-319-78127-3_6.
Cancer Genome Atlas Research Network. Integrated genomic characterization of papillary thyroid carcinoma. Cell. 2014;159(3):676–90. https://doi.org/10.1016/j.cell.2014.09.050.
Huang Y, Prasad M, Lemon WJ, Hampel H, Wright FA, Kornacker K, et al. Gene expression in papillary thyroid carcinoma reveals highly consistent profiles. Proc Natl Acad Sci U S A. 2001;98(26):15044–9. https://doi.org/10.1073/pnas.251547398.
Han J, Chen M, Wang Y, Gong B, Zhuang T, Liang L, et al. Identification of biomarkers based on differentially expressed genes in papillary thyroid carcinoma. Sci Rep. 2018;8(1):9912. https://doi.org/10.1038/s41598-018-28299-9.
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. https://doi.org/10.1158/0008-5472.CAN-04-3078.
Oczko-Wojciechowska M, Pfeifer A, Jarzab M, Swierniak M, Rusinek D, Tyszkiewicz T, et al. Impact of the tumor microenvironment on the gene expression profile in papillary thyroid Cancer. Pathobiology. 2020;87(2):143–54. https://doi.org/10.1159/000507223.
Stockwell BR, Friedmann Angeli JP, Bayir H, Bush AI, Conrad M, Dixon SJ, et al. Ferroptosis: a regulated cell death Nexus linking metabolism, redox biology, and disease. Cell. 2017;171(2):273–85. https://doi.org/10.1016/j.cell.2017.09.021.
Dixon SJ, Lemberg KM, Lamprecht MR, Skouta R, Zaitsev EM, Gleason CE, et al. Ferroptosis: an iron-dependent form of nonapoptotic cell death. Cell. 2012;149(5):1060–72. https://doi.org/10.1016/j.cell.2012.03.042.
Galluzzi L, Vitale I, Aaronson SA, Abrams JM, Adam D, Agostinis P, et al. Molecular mechanisms of cell death: recommendations of the nomenclature committee on cell death 2018. Cell Death Differ. 2018;25(3):486–541. https://doi.org/10.1038/s41418-017-0012-4.
Li J, Cao F, Yin HL, Huang ZJ, Lin ZT, Mao N, et al. Ferroptosis: past, present and future. Cell Death Dis. 2020;11(2):88. https://doi.org/10.1038/s41419-020-2298-2.
Elfarra AA, Krause RJ. Potential roles of flavin-containing monooxygenases in sulfoxidation reactions of l-methionine, N-acetyl-l-methionine and peptides containing l-methionine. Biochim Biophys Acta. 2005;1703(2):183–9. https://doi.org/10.1016/j.bbapap.2004.11.011.
Krueger SK, Williams DE. Mammalian flavin-containing monooxygenases: structure/function, genetic polymorphisms and role in drug metabolism. Pharmacol Ther. 2005;106(3):357–87. https://doi.org/10.1016/j.pharmthera.2005.01.001.
Tang C, Ma J, Liu X, Liu Z. Identification of four immune subtypes in bladder Cancer based on immune gene sets. Front Oncol. 2020;10. https://doi.org/10.3389/fonc.2020.544610.
Subramanian A, Tamayo P, Mootha VK, Mukherjee S, Ebert BL, Gillette MA, et al. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci U S A. 2005;102(43):15545–50. https://doi.org/10.1073/pnas.0506580102.
Lei Z, Yu S, Ding Y, Liang J, Halifu Y, Xiang F, et al. Identification of key genes and pathways involved in vitiligo development based on integrated analysis. Medicine (Baltimore). 2020;99(31):e21297. https://doi.org/10.1097/MD.0000000000021297.
Kang K, Xie F, Mao J, Bai Y, Wang X. Significance of tumor mutation burden in immune infiltration and prognosis in cutaneous melanoma. Front Oncol. 2020;10:573141. https://doi.org/10.3389/fonc.2020.573141.
Kanehisa M, Goto S. KEGG: Kyoto encyclopedia of genes and genomes. Nucleic Acids Res. 2000;28(1):27–30. https://doi.org/10.1093/nar/28.1.27.
Kanehisa M, Sato Y, Kawashima M, Furumichi M, Tanabe M. KEGG as a reference resource for gene and protein annotation. Nucleic Acids Res. 2016;44(D1):D457–62. https://doi.org/10.1093/nar/gkv1070.
Liao Y, Wang J, Jaehnig EJ, Shi Z, Zhang B. WebGestalt 2019: gene set analysis toolkit with revamped UIs and APIs. Nucleic Acids Res. 2019;47(W1):W199–205. https://doi.org/10.1093/nar/gkz401.
Cabanillas ME, McFadden DG, Durante C. Thyroid cancer. Lancet. 2016;388(10061):2783–95. https://doi.org/10.1016/S0140-6736(16)30172-6.
Bergdorf K, Ferguson DC, Mehrad M, Ely K, Stricker T, Weiss VL. Papillary thyroid carcinoma behavior: clues in the tumor microenvironment. Endocr Relat Cancer. 2019;26(6):601–14. https://doi.org/10.1530/ERC-19-0074.
Fisher SB, Perrier ND. The incidental thyroid nodule. CA Cancer J Clin. 2018;68(2):97–105. https://doi.org/10.3322/caac.21447.
Antonelli A, La Motta C. Novel therapeutic clues in thyroid carcinomas: the role of targeting cancer stem cells. Med Res Rev. 2017;37(6):1299–317. https://doi.org/10.1002/med.21448.
Balachandran VP, Gonen M, Smith JJ, DeMatteo RP. Nomograms in oncology: more than meets the eye. Lancet Oncol. 2015;16(4):e173–80. https://doi.org/10.1016/S1470-2045(14)71116-7.
Zeiger MA, Schneider EB. BRAF V600E mutation independently predicts central compartment lymph node metastasis in patients with papillary thyroid cancer. Ann Surg Oncol. 2013;20(1):3–4. https://doi.org/10.1245/s10434-012-2614-x.
Dong S, Huo H, Mao Y, Li X, Dong L. A risk score model for the prediction of osteosarcoma metastasis. FEBS Open Bio. 2019;9(3):519–26. https://doi.org/10.1002/2211-5463.12592.
Hassannia B, Vandenabeele P, Vanden Berghe T. Targeting Ferroptosis to Iron out Cancer. Cancer Cell. 2019;35(6):830–49. https://doi.org/10.1016/j.ccell.2019.04.002.
Doll S, Freitas FP, Shah R, Aldrovandi M, da Silva MC, Ingold I, et al. Xavier da Silva TN, Panzilius E, Scheel CH et al: FSP1 is a glutathione-independent ferroptosis suppressor. Nature. 2019;575(7784):693–8. https://doi.org/10.1038/s41586-019-1707-0.
Chen X, Kang R, Kroemer G, Tang D. Broadening horizons: the role of ferroptosis in cancer. Nat Rev Clin Oncol. 2021;18(5):280–96. https://doi.org/10.1038/s41571-020-00462-0.
Yan HF, Zou T, Tuo QZ, Xu S, Li H, Belaidi AA, et al. Ferroptosis: mechanisms and links with diseases. Signal Transduct Target Ther. 2021;6(1):49. https://doi.org/10.1038/s41392-020-00428-9.
Tang D, Chen X, Kang R, Kroemer G. Ferroptosis: molecular mechanisms and health implications. Cell Res. 2021;31(2):107–25. https://doi.org/10.1038/s41422-020-00441-1.
Angeli JPF, Shah R, Pratt DA, Conrad M. Ferroptosis inhibition: mechanisms and opportunities. Trends Pharmacol Sci. 2017;38(5):489–98. https://doi.org/10.1016/j.tips.2017.02.005.
Bersuker K, Hendricks JM, Li Z, Magtanong L, Ford B, Tang PH, et al. The CoQ oxidoreductase FSP1 acts parallel to GPX4 to inhibit ferroptosis. Nature. 2019;575(7784):688–92. https://doi.org/10.1038/s41586-019-1705-2.
Ubellacker JM, Tasdogan A, Ramesh V, Shen B, Mitchell EC, Martin-Sandoval MS, et al. Lymph protects metastasizing melanoma cells from ferroptosis. Nature. 2020;585(7823):113–8. https://doi.org/10.1038/s41586-020-2623-z.
Maher P, Currais A, Schubert D. Using the Oxytosis/Ferroptosis pathway to understand and treat age-associated neurodegenerative diseases. Cell Chem Biol. 2020;27(12):1456–71. https://doi.org/10.1016/j.chembiol.2020.10.010.
Lei X, Lei Y, Li JK, Du WX, Li RG, Yang J, et al. Immune cells within the tumor microenvironment: biological functions and roles in cancer immunotherapy. Cancer Lett. 2020;470:126–33. https://doi.org/10.1016/j.canlet.2019.11.009.
Wu F, Li F, Lin X, Xu F, Cui RR, Zhong JY, et al. Exosomes increased angiogenesis in papillary thyroid cancer microenvironment. Endocr Relat Cancer. 2019;26(5):525–38. https://doi.org/10.1530/ERC-19-0008.
La Rosa P, Petrillo S, Turchi R, Berardinelli F, Schirinzi T, Vasco G, et al. The Nrf2 induction prevents ferroptosis in Friedreich's ataxia. Redox Biol. 2021;38:101791. https://doi.org/10.1016/j.redox.2020.101791.
Qiu Y, Cao Y, Cao W, Jia Y, Lu N. The application of Ferroptosis in diseases. Pharmacol Res. 2020;159:104919. https://doi.org/10.1016/j.phrs.2020.104919.
Dong C. Cytokine regulation and function in T cells. Annu Rev Immunol. 2021;39(1):51–76. https://doi.org/10.1146/annurev-immunol-061020-053702.
Kaymak I, Williams KS, Cantor JR, Jones RG. Immunometabolic interplay in the tumor microenvironment. Cancer Cell. 2021;39(1):28–37. https://doi.org/10.1016/j.ccell.2020.09.004.
Chuang TC, Chuang AY, Poeta L, Koch WM, Califano JA, Tufano RP. Detectable BRAF mutation in serum DNA samples from patients with papillary thyroid carcinomas. Head Neck. 2010;32(2):229–34. https://doi.org/10.1002/hed.21178.
Xing M. BRAF mutation in papillary thyroid cancer: pathogenic role, molecular bases, and clinical implications. Endocr Rev. 2007;28(7):742–62. https://doi.org/10.1210/er.2007-0007.
Huang M, Yan C, Xiao J, Wang T, Ling R. Relevance and clinicopathologic relationship of BRAF V600E, TERT and NRAS mutations for papillary thyroid carcinoma patients in Northwest China. Diagn Pathol. 2019;14(1):74. https://doi.org/10.1186/s13000-019-0849-6.
Fakhruddin N, Jabbour M, Novy M, Tamim H, Bahmad H, Farhat F, et al. BRAF and NRAS mutations in papillary thyroid carcinoma and concordance in BRAF mutations between primary and corresponding lymph node metastases. Sci Rep. 2017;7(1):4666. https://doi.org/10.1038/s41598-017-04948-3.
Schulten HJ, Salama S, Al-Ahmadi A, Al-Mansouri Z, Mirza Z, Al-Ghamdi K, et al. Comprehensive survey of HRAS, KRAS, and NRAS mutations in proliferative thyroid lesions from an ethnically diverse population. Anticancer Res. 2013;33(11):4779–84.
Jasim A, Mohammed A, Ibrahim A. Beta-2-Microglobulin as a Marker in Patients with Thyroid Cancer. Iraqi Postgraduate Med J. 2019;18(1):6.
Middha S, Yaeger R, Shia J, Stadler ZK, King S, Guercio S, et al. Majority of B2M-mutant and -deficient colorectal carcinomas achieve clinical benefit from immune checkpoint inhibitor therapy and are microsatellite instability-high. JCO Precis Oncol. 2019;3:PO.18.00321. https://doi.org/10.1200/PO.18.00321.
Castro A, Ozturk K, Pyke RM, Xian S, Zanetti M, Carter H. Elevated neoantigen levels in tumors with somatic mutations in the HLA-A, HLA-B, HLA-C and B2M genes. BMC Med Genet. 2019;12(Suppl 6):107. https://doi.org/10.1186/s12920-019-0544-1.
Liang W, Sun F. Prognostic Alternative mRNA Splicing in Adrenocortical Carcinoma. Front Endocrinol (Lausanne). 2021;12:538364. https://doi.org/10.3389/fendo.2021.538364.
Lin L, Li X, Pan C, Lin W, Shao R, Liu Y, et al. ATXN2L upregulated by epidermal growth factor promotes gastric cancer cell invasiveness and oxaliplatin resistance. Cell Death Dis. 2019;10(3):173. https://doi.org/10.1038/s41419-019-1362-2.
Chen Y, Hu Y, Zhang H, Peng C, Li S. Loss of the Alox5 gene impairs leukemia stem cells and prevents chronic myeloid leukemia. Nat Genet. 2009;41(7):783–92. https://doi.org/10.1038/ng.389.
Chen Y, Li D, Li S. The Alox5 gene is a novel therapeutic target in cancer stem cells of chronic myeloid leukemia. Cell Cycle. 2009;8(21):3488–92. https://doi.org/10.4161/cc.8.21.9852.
Qi S, Su L, Li J, Zhang C, Ma Z, Liu G, et al. Arf6-driven endocytic recycling of CD147 determines HCC malignant phenotypes. J Exp Clin Cancer Res. 2019;38(1):471. https://doi.org/10.1186/s13046-019-1464-9.
Brenne AT, Fagerli UM, Shaughnessy JD Jr, Våtsveen TK, Rø TB, Hella H, et al. High expression of BCL3 in human myeloma cells is associated with increased proliferation and inferior prognosis. Eur J Haematol. 2009;82(5):354–63. https://doi.org/10.1111/j.1600-0609.2009.01225.x.
Legge DN, Chambers AC, Parker CT, Timms P, Collard TJ, Williams AC. The role of B-cell Lymphoma-3 (BCL-3) in enabling the hallmarks of cancer: implications for the treatment of colorectal carcinogenesis. Carcinogenesis. 2020;41(3):249–56. https://doi.org/10.1093/carcin/bgaa003.
Wakefield A, Soukupova J, Montagne A, Ranger J, French R, Muller WJ, et al. Bcl3 selectively promotes metastasis of ERBB2-driven mammary tumors. Cancer Res. 2013;73(2):745–55. https://doi.org/10.1158/0008-5472.CAN-12-1321.
Maldonado V, Melendez-Zajgla J. Role of Bcl-3 in solid tumors. Mol Cancer. 2011;10(1):152. https://doi.org/10.1186/1476-4598-10-152.
Ostendorf BN, Bilanovic J, Adaku N, Tafreshian KN, Tavora B, Vaughan RD, et al. Common germline variants of the human APOE gene modulate melanoma progression and survival. Nat Med. 2020;26(7):1048–53. https://doi.org/10.1038/s41591-020-0879-3.
Chang NW, Chen DR, Wu CT, Aouizerat BE, Chen FN, Hung SJ, et al. Influences of apolipoprotein E polymorphism on the risk for breast cancer and HER2/neu status in Taiwan. Breast Cancer Res Treat. 2005;90(3):257–61. https://doi.org/10.1007/s10549-004-4656-7.
Zhao Z, Zou S, Guan X, Wang M, Jiang Z, Liu Z, et al. Apolipoprotein E overexpression is associated with tumor progression and poor survival in colorectal Cancer. Front Genet. 2018;9:650. https://doi.org/10.3389/fgene.2018.00650.
Foshan First People’s Hospital Deng Feng project 2019A008, Foshan municipal Government (2015AG1003), Guangdong Provincial Government (2016A020213001) of China, MCRs project of Genecology Research Centre of University of the Sunshine Coast, and National Natural Science Foundation of China (31971355).
Ethics approval and consent to participate
No specific permits were required, all data were downloaded from online data repository.
Consent for publication
The authors declare that have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Table S1. Enrichment analysis of KEGG pathways by using GSVA.
Table S2. The comparison of immune cell contents identified in the two subtypes.
Fig. S1. Co-expression analysis of 1747 immune genes in the TCGA cohort. (A) Calculation of the scale-free fit index of various soft-thresholding powers and the analysis of the mean connectivity of various soft-thresholding powers. (B) Median rank and Z-summary statistics in the module preservation process. The plot shows the module position in the test dataset based on the Median rank(left). The plot illustrates the analysis of the Z-summary between different modules in the test dataset(right). (C) Heatmap of module trait relationships between gene modules and traits of THCA. The correlation coefficient in each cell represents the correlation between the modules and traits, which increases in size from blue to red. The corresponding P-value is annotated. The module–trait relationships were demonstrated by correlation values and P-values (In parenthesis) with a range of colors; the degree of correlation between modules and clinical features is shown. Rows are module eigengene (ME) regards to each module, and the columns indicate traits.
Table S3. The full list of significant differentially expressed genes.
Table S4. The statistical results of Unicox and Multicox analysis.
About this article
Cite this article
Lin, R., Fogarty, C.E., Ma, B. et al. Identification of ferroptosis genes in immune infiltration and prognosis in thyroid papillary carcinoma using network analysis. BMC Genomics 22, 576 (2021). https://doi.org/10.1186/s12864-021-07895-6
- Thyroid papillary carcinoma
- Immune-related genes
- Tumour microenvironment
- Ferroptosis genes
- Prognostic model