Construction and validation of risk models of prognostic genes associated with parthanatos in papillary thyroid carcinoma based on bioinformatics

Construction and validation of risk models of prognostic genes associated with parthanatos in papillary thyroid carcinoma based on bioinformatics | Research Square window.SnipcartSettings = { analytics: { enabled: false } }; (function() { var accessVector = localStorage.getItem('access_vector') || ''; window.dataLayer = window.dataLayer || []; if (accessVector) { window.dataLayer.push({ user: { profile: { profileInfo: { snid: accessVector } } } }); } })(); (function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start':new Date().getTime(),event:'gtm.js'});var f=d.getElementsByTagName(s)[0],j=d.createElement(s),dl=l!='dataLayer'?'&l='+l:'';j.async=true;j.src='https://www.googletagmanager.com/gtm.js?id='+i+dl;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-K279D39R'); Browse Preprints In Review Journals COVID-19 Preprints AJE Video Bytes Research Tools Research Promotion AJE Professional Editing AJE Rubriq About Preprint Platform In Review Editorial Policies Our Team Advisory Board Help Center Sign In Submit a Preprint Cite Share Download PDF Research Article Construction and validation of risk models of prognostic genes associated with parthanatos in papillary thyroid carcinoma based on bioinformatics Rui Wang, Li Zhang, Shuxin Wen This is a preprint; it has not been peer reviewed by a journal. https://doi.org/ 10.21203/rs.3.rs-6673087/v1 This work is licensed under a CC BY 4.0 License Status: Under Review Version 1 posted 10 You are reading this latest preprint version Abstract Object: This study aimed to elucidate the role of parthanatos-related genes (PRGs) in papillary thyroid carcinoma (PTC) and construct a prognostic risk model to guide personalized treatment. Methods Using the GSE33630 dataset, differentially expressed PRGs were identified and analyzed via weighted gene co-expression network analysis (WGCNA) to pinpoint key module genes. Regression analysis selected seven prognostic genes for risk model construction. The model’s performance was validated, and a nomogram was developed for survival prediction. Further analyses included clinical feature correlations, immune infiltration, drug sensitivity, gene set enrichment analysis (GSEA), and experimental validation via RT-qPCR. Results Seven prognostic genes (TSHZ3, SERGEF, AKAP12, SGPP2, ASGR1, AK1, PELI2) were identified. The risk model demonstrated robust predictive accuracy, stratifying patients into high- and low-risk groups with significant survival differences. GSEA revealed 29 enriched pathways (e.g., ribosome, focal adhesion), while immune infiltration analysis highlighted CD56 + NK cells and AK1 as key immune correlates. Drug sensitivity screening identified 111 differential therapeutics. Functional analysis indicated AKAP12 had the strongest functional similarity among prognostic genes. Conclusion This study comprehensively mapped PRGs in PTC, established a validated risk model, and provided insights into immune-microenvironment interactions and therapeutic targets, advancing precision oncology for PTC. Papillary thyroid carcinoma Parthanatos Prognostic gene Risk model Immune checkpoint inhibitors Figures Figure 1 Figure 2 Figure 3 Figure 4 Figure 5 Figure 6 Figure 7 Figure 8 Figure 9 Figure 10 Introduction Thyroid cancer is a malignant tumor originating from the follicular epithelium or parafollicular epithelial cells of the thyroid gland, and is the most common malignant tumor in the head and neck. In recent years, the incidence of thyroid cancer in the world has increased rapidly, and thyroid cancer in China has continued to increase at a rate of 20% per year [ 1 ] . According to the origin and differentiation of tumors, thyroid cancer is further divided into: Papillary Thyroid Carcinoma (PTC), Follicular Thyroid Carcinoma (FTC), Medullary Thyroid Carcinoma (MTC), and poorly differentiated thyroid carcinoma (PDTC) and anaplastic thyroid cancer (ATC), of which PTC is the most common, accounting for about 90% of the entire thyroid cancer. PTC generally has a good prognosis, however, distant metastasis and poor prognosis can also be found in some aggressive subtypes of PTC [ 2 – 3 ] . Therefore, it is of great significance to study the development mechanism of PTC for the treatment of thyroid cancer. The prognostic of PTC is associated with age, sex, loss of histologic differentiation, tumor size, presence of carcinomatous lymphangitis, extrathyroid extension, and presence of metastasis [ 4 – 5 ] . The prognostic prediction studies of PTC are mostly conducted by molecular markers, among which BRAF [ 6 – 7 ] and TERT [ 8 – 9 ] are the main studies. However, the prognostic prediction of BRAF is still controversial [ 10 – 11 ] . Therefore, the exploration of other prognostic genes and their molecular mechanisms provides a new idea for prognosis assessment and individualized treatment of PTC patients. The biological homeostasis of an organism depends on the balance between cell proliferation and death. Parthanatos is a poly(ADP-ribose) polymerase 1 (PARP1) hyperactivation dependent and cysteine-independent cell death which is distinct from apoptosis, programmed necrosis, or other known forms of cell death [ 12 ] . Parthanatos is a multi-step pathway that plays a crucial role in tumorigenesis. Once the DNA is damaged, it leads to excessive activation of PARP, then PARP-1 produces free PAR by PARG-mediated hydrolyzation, inducing the release of apoptosis-inducing factor(AIF) from the outer mitochondrial membrane [ 12 – 14 ] . Once AIF is released, it can bind to macrophage migration inhibitory factor (MIF), translocates to the nucleus, and eventually lead to a large number of DNA fragments [ 15 – 16 ] . Parthanatos has been studied in hepatocellular carcinoma [ 17 ] , breast cancer [ 18 ] , prostate cancer [ 19 ] , cervical cancer [ 20 ] and other cancers, however, the function and mechanism of parthanatos in PTC remain unclear and require further extensive research. In recent years, it has become increasingly mature to screen key genes and study their molecular mechanisms through bioinformatics [ 21 – 23 ] . In this study, prognostic genes associated with parthanatos in PTC was identified with bioinformatics through two public databases, then, their potential molecular mechanisms was been explored, and a prognostic risk model was constructed to predict the survival rate of PTC patients, providing effective targets for personalized treatment and prognostic intervention for PTC patients. Materials and methods 2.1 Dataset selection and description We obtained genomic expression profiles and clinical feature data (including age, gender, subtype, grade, stage, lymphovascular invasion, lymph nodes positive by he, etc.) of the TCGA-THCA dataset from the Cancer Genome Atlas Program database (TCGA, https://portal.gdc.cancer.gov/ ), comprising 510 thyroid carcinoma (THCA) tissues and 58 normal tissue samples, which were utilized for prognostic analysis and validation. The samples from the TCGA-THCA dataset were apportioned randomly into two subsets in a ratio of 6 to 4, with the larger segment designated as the training set and the remainder allocated for validation purposes. The chip data of papillary thyroid carcinoma (PTC) patients and normal thyroid tissue samples (49/45) from the GSE33630 (GPL570 platform) dataset were procured from the Gene Expression Omnibus database (GEO, https://www.ncbi.nlm.nih.gov/geo/ ) database for differential expression analysis. Additionally, 9 parthanatos-related genes (PRGs) were identified from prior research [ 24 ] . 2.2 Differential expression analysis To delineate the landscape of differentially expressed genes (DEGs) in PTC patients relative to the control cohort, the expression matrix derived from the GSE33630 dataset was meticulously analyzed using the "limma " package (version 3.54.0) [ 25 ] . This analytical sieve employed a stringent criterion of adj.p.value 1 to discern DEGs. To gain an overall understanding of the distribution of DEGs, the "ggplot2" (version = 3.4.1) package was used to draw volcano plots for visualizing the differential genes and the "ComplexHeatmap" (version 2.14.0) [ 26 ] package was utilized to generate a heatmap illustrative of the top 10 genes demonstrating significant up-regulation and the top 10 genes exhibiting marked down-regulation, meticulously sorted by the respective fold changes in their expression levels. 2.3 Weighted correlation network analysis (WGCNA) In the TCGA-THCA disease samples, ssGSEA scores were calculated based on 9 PRGs using the "GSVA" (version 1.46.0) [ 27 ] package. Tumor samples were categorically segregated into high and low ssGSEA score cohorts, demarcated by an optimal threshold value derived from the ssGSEA scoring metric. Survival disparities between the high and low ssGSEA score were delineated through the construction of Kaplan-Meier (KM) survival curves. The difference in KM curves is considered significant when p < 0.02. Then, the expression matrix of disease samples in TCGA-THCA was analyzed utilizing WGCNA through the implementation of the "WGCNA" (version 1.71) [ 28 ] package. To commence the analytical process, an initial hierarchical clustering was performed on the expression data from all samples within the TCGA-THCA dataset, employing the Euclidean distance metric to assess the similarity of expression levels. This step was crucial for identifying and excluding any outliers among the samples, thereby refining the dataset for subsequent analysis. Subsequently, to ensure the constructed network adhered to a scale-free topology, an optimal soft threshold needed to be determined before constructing the co-expression network and set to r2 = 0.85. Ultimately, an unscaled network was crafted utilizing the chosen soft threshold parameter, effectively segmenting the genes into distinct modules, each adorned with a unique color identifier. Using the parthanatos related ssGSEA scores as the phenotype, associations between phenotypes and modules were established. This involved calculating the correlation coefficient between modules and phenotype traits, followed by generating a heatmap of correlations. The most positively and negatively associated modules with the parthanatos related ssGSEA scores (|cor| >0.3, p < 0.05) were identified in this study as key modules. Utilizing the "ggVenn" package (version 0.1.9), an overlap analysis was conducted between the DEGs and the key module genes, resulting in the identification of differential genes associated with parthanatos in PTC, referred to as DE-PRGs. 2.4 Gene function enrichment and protein-protein interaction (PPI) analysis We employed the "clusterProfiler" (version 4.7.1.3) [ 29 ] package to perform Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses on the DE-PRGs. These analyses were employed with a stringent cutoff of p-values < 0.05 to curate a concise list of significantly enriched biological processes and pathways. To establish a PPI network for DE-PRGs, genetic entities were submitted to the STRING database ( http://www.string-db.org/ ) to get results (interaction score > 0.15). Subsequently, the outcomes were imported into Cytoscape (version 3.7.2) [ 30 ] for graphical representation and analysis. 2.5 Construction of risk models In order to identify the genes associated with the prognosis of PTC patients, a univariate Cox regression analysis was performed on DE-PRGs within the TCGA-THCA training dataset using the "survival" (version 3.5-3) package. The data used consisted of both gene expression information of THCA samples and training set data from TCGA-THCA. Prior to conducting univariate Cox regression analysis, the assumption of proportional hazards (PH) for the risk model results was tested using the R software's cox.zph function (p > 0.05) for DE-PRGs. After excluding those that failed the PH test, the remaining DE-PRGs were used to construct univariate Cox regression model (p < 0.05 & HR ≠ 1). High HR values were excluded to avoid influencing subsequent analysis, and the remaining survival genes were then integrated into our subsequent analytical procedures. The Forestplot (version 2.0.1) package was employed to generate a graphical forest plot, elucidating the outcomes of the univariate Cox regression analysis in a visually intuitive manner. Thereafter, the remaining genes underwent a LASSO regression analytical process utilizing the "glmnet" (version 4.1-6) [ 31 ] package to delineate prognostic genes and to establish a risk prediction model. The analysis was conducted with the "family" parameter specified as Cox, thereby incorporating the Cox proportional hazards model, and prioritizing strongly correlated features for inclusion. The model was considered optimal when the model error was minimized, and the corresponding lambda value at this point was referred to as lambdamin. 2.6 Evaluation of risk models In the TCGA-THCA training set, individual patients with PTC were assigned a score derived from the expression of prognostic genes and the risk coefficients ascertained through LASSO regression analysis. This segmentation process resulted in the classification of patients into distinct high- and low-risk cohorts, determined by an empirically derived optimal cutoff for the risk scoring mechanism. The computation of the risk score was achieved through a linear combination of the gene expression levels and their respective weights within the predictive framework. The formulaic representation is as follows: Where, "risk score" represents the risk score, "coef" signifying the risk coefficients attributed to each specific gene, and "expr" represents the expression of the respective genes. The "survminer" (version 0.4.9) package was employed to delineate the survival curves for overall survival (OS) among the two risk cohorts. Furthermore, to gauge the precision of risk prognostication, receiver operating characteristic (ROC) curves were generated employing the "survivalROC" (version 0.4.9) package based on TCGA-THCA training set data, with survival time points at 3/5/7 years. The distinct characteristics separating the two risk cohorts were ascertained by assessing the transcriptomic levels of prognostic genes was analyzed based on the threshold criteria |cor|>0.3 and p<0.05. 2.7 Validation of the risk model To gauge the efficacy of the risk model, it underwent rigorous validation within the TCGA-THCA validation set and TCGA-THCA dataset. The risk model derived from the aforementioned training cohort was applied, and using the optimal cutoff value of risk score, to stratify the PTC patients within the TCGA-THCA validation subset and TCGA-THCA dataset into distinct high- and low-risk cohorts. Subsequently, the "survminer" package was deployed to illustrate the survival and ROC curves for overall survival, thereby comparing the prognostic outcomes between the two risk cohorts. Then, the differences between the two risk cohorts in the validation set were observed drawing on the expression patterns. Simultaneously, the "psych" (version = 2.2.9) package was employed to calculate the Spearman correlation coefficient between prognostic genes (|cor|>0.3 and p<0.05). 2.8 Construction and evaluation of nomogram Building upon the foundational analysis of prognostic gene expression, survival was considered as the outcome event. A predictive nomogram was then artfully constructed through the "rms" (version 6.5.0) package. For the purposes of validation, calibration curves were meticulously plotted with the "regplot" (version 1.1), while ROC curves were delineated by the "survivalROC" (AUC < 0.6 ), thereby assessing the nomogram's clinical utility in forecasting the progression of PTC. 2.9 Correlation analysis of survival differences and clinical features Within the TCGA-THCA training set, the Wilcoxon test was employed to elucidate the disparities in survival outcomes among the high- and low-risk cohorts across various subgroups categorized by distinct clinical feature (including age, gender, pathological stage, T stage, N stage, and M stage ), with a threshold set at p < 0.05. K-M curves were plotted using the "survminer" package for visualization. 2.10 Gene set enrichment analysis (GSEA) A differential expression analysis was performed on the respective cohort within the TCGA-THCA training set, utilizing the "DESeq2" (version 1.38.0) [ 32 ] package, to delineate the biological pathways underlying the high- and low-risk cohorts. According to the ranking of log 2 FC, the "clusterProfiler" package was employed to conduct GSEA utilizing the gene set “c2.cp.kegg.v2023.1.Hs.symbols.gm” derived from the latest version of the Molecular Signatures Database (MSigDB, https://www.gsea-msigdb.org/gsea/msigdb/index.jsp)(|NES| >1 and p < 0.05). 2.11 Analysis of the immune microenvironment In the TCGA-THCA training set, the PTC patients underwent immune infiltration analysis utilizing the ssGSEA algorithm, exploring the immune cell infiltration status in PTC patients. Firstly, we derived the enrichment scores for 28 immune cells in the different risk cohorts. Subsequently, the Wilcoxon test was utilized to ascertain the presence of statistically significant disparities in the distribution of immune cell populations between the two risk cohorts within the TCGA-THCA training dataset. Then, for all samples in the training set of TCGA-THCA, the "corrplot" (0.92) was utilized for visualize the correlation matrix to perform Spearman correlation analysis on differential immune cells. Simultaneously, Spearman correlation analysis was conduct on differential immune cells, risk scores, and prognostic genes, considering the results to be significant if |cor| >0.3 and p < 0.05. Ultimately, the quantification of stromal and immune cells in tumor tissue was achieved through the use of the "ESTIMATE" algorithm, thereby enabling the prediction of the immune score, stromal score, and ESTIMATE score. Additionally, it facilitated the estimation of tumor purity, as well as the levels of immune and stromal infiltration within the malignant tumor tissue. Differences in the distribution of three scores were analyzed between two risk cohorts using the Wilcoxon test. To discern the disparities in the distribution of the three scores between the high- and low-risk cohorts, the Wilcoxon test was employed. 2.12 Analysis of immune checkpoint inhibitors (ICIs) We obtained 14 ICIs genes (ASXL1, BCL2, CD33, CD47, CHEK1, PLK1, CTLA4, DOT1L, FLT3, IDH1, IDH2, MCL1, PDCD1, MDM2) from the pervoius studies [ 33 – 34 ] . Within the TCGA-THCA training dataset, the Wilcoxon test was invoked to delineate the disparities in ICIs between the high- and low-risk cohorts. Those ICIs exhibiting statistically significant variances (p < 0.05) across the two risk cohorts were chosen for further exploration. Following this, the "cor" function within the R programming environment was utilized to elucidate the Spearman correlation analysis, elucidating the relationship among the risk score and the differential expression of ICIs. 2.13 Drug sensitivity analysis To ascertain the disparities in drug responsiveness between the high- and low-risk cohorts, the Genomics of Drug Sensitivity in Cancer (GDSC, https://www.cancerrxgene.org/ ) database furnished the half-maximal inhibitory concentrations (IC 50 ) and PTC drug. The "pRRophetic" (version 0.5) [ 25 ] package was deployed to forecast the IC 50 for each tumor specimen within the TCGA-THCA training cohort. Subsequently, the divergences in drug sensitivity between the two risk cohorts were subjected to Wilcoxon test (p < 0.05). 2.14 Functional analysis of prognostic genes In order to discern other genes associated with prognostic gene function, the prediction of genes and gene function associated with prognosis was achieved through the use of the the GeneMANIA database ( http://genemania.org ). The functional similarity between prognosis genes was examined employing the "GOSemSim" (version 2.24.0) [ 35 ] package, with a score greater than 0.5 being determined as indicating high relevance. 2.15 Construction of the lncRNA-miRNA-mRNA network for prognostic genes The miRNAs that regulated prognostic genes were forecasted based on the TarBase v8.0( http://mirtarbase.mbc.nctu.edu.tw/ ) database and miRDB database ( https://mirdb.org/ ) using the Network Analyst platform. The "ggVenn" (version 0.1.9) package was utilized to intersect the miRNAs predicted from the two databases, obtaining key miRNAs. Subsequently, the miRNet database ( https://www.mirnet.ca/miRNet/home.xhtml ) was employed to forecast the lncRNAs and circRNAs of the key miRNAs. The top 3 lncRNAs and circRNAs were selected as the regulators of the miRNAs. The Cytoscape (version 3.7.2) [ 30 ] platform was utilized to build a comprehensive lncRNA-miRNA-mRNA regulatory network. 2.16 Validation of gene expression for prognostic genes The expression differences of prognostic genes in PTC cohort and normal cohort were analyzed employing the Wilcoxon test based on the entire sample of the TCGA-THCA dataset (p < 0.05), to further corroborate the expression of prognostic genes. Finally, immunohistochemical staining maps of prognostic genes were procured from the Human Protein Atlas database (HPA, http://www.proteinalas.org/ ) to visualize the differences in protein levels of prognostic genes between the PTC cohort and normal cohort. 2.17 Validation of reverse transcription-quantitative real-time Polymerase chain reaction (RT-qPCR) The expression of genes was validated via RT-qPCR. The tissue samples were obtained from Third Hospital of Shanxi Medical University, and it has been approved by the ethics committee (Approval notice number: YXLL-2024-173) and all the participants have signed the informed consent. The results were plotted using Graphpad Prism 10. The TRizol reagent kit was used to extract RNA from 5 pairs of frozen tissue samples, with samples 1–5 designated as the Control group and samples 6–10 as the PTC group. All experimental procedures were carried out according to the instructions, resulting in total RNA. 1 µl of the extracted RNA was taken for concentration detection using the Nano Photometer N50, and the purity/concentration was recorded to calculate the amount for subsequent reverse transcription steps. The remaining RNA was either immediately subjected to reverse transcription or stored at -80°C. The SweScript First Strand cDNA Synthesis SuperMix for PCR from Servicebio Company was used to reverse transcribe the RNA into cDNA, following the instructions. The cDNA was diluted 5–20 times with ddH 2 O (RNase/ARase free), with 3ul of cDNA, 5ul of 2xUniversal Blue SYBR Green qPCR Master Mix, 1ul of forward primer (10µM), 1ul of reverse primer (10µM). Subsequently, 40 cycles of reaction were carried out using the CFX96 real-time quantitative PCR instrument, with the program detailed in Supplementary Table 1 . The primer sequences for the TSHZ3, SERGEF, AKAP12, ASGR1, AK1, PELI2 and SGPP2 genes can be found in Supplementary Table 2 , and GAPDH was used as the reference gene. The 2-ΔΔCT method was used to calculate the relative gene expression levels. 2.18 Statistical analysis All analyses were conducted using the R programming language (version 4.2.3). The p-values and r-values for correlation analysis were calculated using non-parametric Spearman correlation coefficient. The disparities between groups were scrutinized utilizing the Wilcoxon test, with statistical significance being attributed to findings where p < 0.005. Results 3.1 Indentification of 82 DE-PRGs Through differential expression analysis of the PTC and the control cohort, 1,007 DEGs were identified (557 up-regulated genes; 450 down-regulated genes) ( Fig. 1 A ) . Additionally, the top 10 up- and down-regulated genes were displayed based on log 2 FC sorting ( Fig. 1 B ) . In the disease samples of TCGA-THCA, the parthanatos scores of the disease cohort were computed for the 9 PRGs obtained from previous studies. The tumor specimens were categorically segregated into distinct cohorts, distinguished by high and low score brackets, based on the optimally determined cutoff threshold derived from the ssGSEA scoring metrics. The significant difference in the high- and low-risk cohorts was observed through the K-M curve (p < 0.02) (Fig. 1 C), indicating that it was meaningful for us to use DE-PRGs genes as the background gene set for scoring. Subsequently, WGCNA was undertaken utilizing the expression matrix of disease samples within the TCGA-THCA dataset. Firstly, all samples were clustered, and no outlier samples were found ( Fig. 1 D ) . Then, sample traits (ssGSEA scores) were introduced, and hierarchical clustering was performed again using the Euclidean distance of sample expression levels ( Fig. 1 E ) . In the pursuit of constructing a scale-free gene co-expression network, the key step of finding the optimal soft threshold was carried out before building the co-expression network. Upon achieving a scale-free fit index of 0.85, the soft threshold with a connectivity close to 0 and a β value of 10 is chosen for screening ( Fig. 1 F ) . Utilizing the chosen soft threshold, the scale-free network was constructed, and these genes were divided into several modules, each marked with a different color. A total of 11 co-expression modules were identified through WGCNA analysis ( Fig. 1 G ) (excluding the grey module to which unclassified genes belong). Using the parthanatos related ssGSEA score as a phenotype, the module eigengenes and the correlation coefficient matrix between phenotype and traits were obtained ( Fig. 1 H ) . The blue module (cor = 0.51,p = 1e − 35 ) and the turquoise module (cor=-0.38, p = 9e − 19 ) with the highest significant correlation with the parthanatos related ssGSEA score, denoted as the key modules, were chosen. By selecting all module genes, 6,628 module genes were obtained. Finally, the obtained DEGs were intersected with key module genes, resulting in 82 intersecting genes, referred to as DE-PRGs ( Fig. 1 I ) . 3.2 Functional enrichment and protein interaction analysis of DE-PRGs In order to explore the biological functions and signaling pathways, we conducted GO and KEGG enrichment analyses on the 82 DE-PRGs. In the GO enrichment analysis, a total of 177 GO biological functions were enriched (p-value < 0.05 )(144 in biological processes (BPs), 16 in molecular functions (MFs), and 17 in cellular components (CCs)) ( Fig. 2 A ) . The enriched GO biological functions mainly included cellular cation homeostasis, divalent inorganic cation homeostasis, and cellular metal ion homeostasis. In the KEGG enrichment analysis, based on the screening criterion of p-value < 0.05, a total of 8 KEGG pathways were enriched ( Fig. 2 B ) . The enriched pathways included regulation of actin cytoskeleton, calcium signaling pathway, thyroid hormone synthesis, prostate cancer, complement and coagulation cascades, virion - ebolavirus and lyssavirus, sulfur metabolism, and thiamine metabolism. Simultaneously, to investigate the protein-level interaction of the DE-PRGs, we deposited the 82 DE-PRGs into the STRING to forge a DE-PRGs protein interaction network (interaction score > 0.15) ( Fig. 2 C ). In the PPI network, after removing isolated proteins, a network of interactions among 73 proteins was obtained, containing 73 nodes and 117 edges. Among them, AUTS2 exhibited the strongest interaction. 3.3 Construction of risk model and identification of prognostic genes In order to discern genes linked to prognosis in PTC patients the univariate Cox regression was conducted on 82 DE-PRGs within the TCGA-THCA training set. Initially, four DE-PRGs that did not pass the PH assumption test (p > 0.05) were excluded, and then the univariate Cox regression model was performed, uncovering a relationship among the expression of 11 genes and OS, with MIPOL1 exhibiting an exceptionally high HR value, potentially influenced by confounding factors and genetic mutations, and was therefore excluded to avoid biasing subsequent analyses. The remaining 10 genes associated with survival were included for further analysis, revealing that among these 10 survival-related genes, SERGEF and AK1 were protective factors, while the others were risk factors ( Fig. 3 A, Table 3 − 1) . Therefore, this study constructed a risk model through LASSO regression analysis, which includes seven genes: TSHZ3, SERGEF, AKAP12, ASGR1, AK1, PELI2 and SGPP2. These genes, were defined as prognostic genes ( Fig. 3 B, 3 C, Table 3 − 2) . When lambdamin was set to 0.007, the corresponding number of genes was 7. Utilizing the transcriptome data of the identified prognostic genes in conjunction with the risk coefficients ascertained through LASSO regression analysis, a customized risk evaluation was allocated for every PTC sample. The algorithmic formula for computing this risk score is as follows: riskscore = TSHZ3 * (1.099) + SERGEF * (-0.450) + AKAP12 * (0.204) + SGPP2 * (0.089) + ASGR1 * (0.211) + AK1 * (-0.005) + PELI2 * (0.510) Using the optimal cutoff value of risk score, PTC cohort were segregated into two distinct groups: a high-risk cohort and a low-risk cohort ( Fig. 3 D, 3 E ) , revealing a notable contrast in survival between the two risk cohorts (p = 0.0001), with lower survival rates among the high-risk cohort patients ( Fig. 3 F ) . Furthermore, the precision of the risk model was assessed utilizing TCGA-THCA training set data. The AUC of the ROC curve consistently exceeded 0.6 in the training set, signifying commendable performance of the risk model ( Fig. 3 G ) . In terms of the expression levels of prognostic genes, no notable disparities were discerned between the two cohorts ( Fig. 3 H ) . However, ASGR1 and PELI2 showed lower expression levels in the training set, while SERGEF and AK1 exhibited higher expression levels, with SGPP2 and SERGEF showing the highest negative correlation ( Fig. 3 I ) . Table 3 − 1 . PH test of 10 prognostic genes ID AKAP12 TSHZ3 F2R ASGR1 PELI2 AK1 SGPP2 CRYBG3 N4BP2 SERGEF PH 0.911 0.708 0.467 0.642 0.795 0.281 0.401 0.975 0.955 0.499 Table 3 − 2 . PRisk coefficient of prognostic genes Gene Coefficient TSHZ3 1.099478866 SERGEF -0.449517532 AKAP12 0.20384777 SGPP2 0.088714538 ASGR1 0.210755007 AK1 -0.004941586 PELI2 0.509768893 3.4 Validation of the risk assessment model and correlation of prognostic genes Within the TCGA-THCA validation set, a marked discrepancy in patient survivability was discerned between the high-risk and low-risk cohorts (p = 0.012), with lower survival rates observed in the high risk cohort ( Fig. 4 A, 4 B ) . The AUC values at 3/5/7 years in the ROC curve were all greater than 0.6, indicating excellent validation results ( Fig. 4 C ) . Similarly, in the TCGA-THCA dataset, a notable contrast in patient survival emerged between the two risk cohorts (p = 0.0001), with lower survival rates observed in the high-risk cohort ( Fig. 4 D, 4 E ) . The AUC values at 3/5/7 years in the ROC curve were all greater than 0.6, indicating excellent validation results ( Fig. 4 F ) . Furthermore, through the observation of the disparate expression of prognostic genes between the two risk cohorts in both datasets, notable distinctions in the expression levels were discerned. In the TCGA-THCA validation set, ASGR1, SGPP2, PELI2 and TSHZ3 showed higher expression in the high-risk cohort contrast to the low-risk cohort ( Fig. 4 G ) ; while in the TCGA-THCA dataset, AK1 and SERGEF showed lower expression in the high-risk cohort contrast to the low-risk cohort ( Fig. 4 H ) . Additionally, analyzing the correlation between prognostic genes revealed that SERGEF and AK1 (cor = 0.53 and P < 0.001) exhibited the highest positive correlation in both datasets ( Fig. 4 I, 4 J ) . 3.5 Establishment and verification of nomogram Utilizing the expression levels of prognostic genes, a nomogram model was constructed with survival as the outcome event ( Fig. 5 A ) . For the nomogram predicting survival risk, the individual points in the graph represents the score corresponding to the different values of each variable. The total points was derived from the sum of the individual points corresponding to all variable values. The higher the total points, the higher the patient's survival rate. Then, the survival probabilities for 3, 5, and 7 years were estimated for PTC patients using the total points. The calibration curve (p > 0.05, the closer the slope is to 1, the better the accuracy of the model's predictions) demonstrated that the nomogram model exhibited commendable prognostic power in predicting the clinical outcomes of PTC ( Fig. 5 B ) . The results from the ROC curve (AUC > 0.6 is considered to indicate model accuracy) revealed that the nomogram model exhibited commendable precision in forecasting the clinical outcomes of PTC ( Fig. 5 C ) . Therefore, it was believed that the nomogram model is effective in predicting PTC clinical outcomes. 3.6 Clinical characterization analysis in high- and low-risk cohorts Based on the TCGA-THCA training set, the survival differences between two risk cohorts were demonstrated in different clinical feature subgroups, including age, gender, pathological stage, T stage, N stage, and M stage (p < 0.05). It was discernible that there existed significant disparities in survival outcomes among two risk cohorts for patients aged over 47, gender, stage3, T3, N, and M0, with too few samples for M1 resulting in no outcome ( Fig. 6 A- 6 F ). It was manifest that the high-risk cohorts exhibit a survival rate that is inferior to that observed among the low-risk cohorts. 3.7 Enrichment of biological pathways and drug sensitivity analyses By conducting GSEA on the two risk cohorts, 29 pathways were ultimately enriched (|NES| >1, p < 0.05), with main enrichments in various signaling pathways such as ribosome, focal adhesion, arachidonic acid metabolism, ECM receptor interaction and oxidative phosphorylation ( Fig. 7 A ) . Through the analysis of differences in drug sensitivity between high- and low-risk cohorts, we found that 111 drugs showed significant differences, such as cisplatin, EHT.1864, embelin, obatoclax. mesylate and VX.680 ( Fig. 7 B ) . The top10 drugs (ranking by p-value) were selected to display the differences. And we found these drugs all had higher IC50 in the low risk cohort contrast to the high risk cohort. 3.8 Analysis of immune infiltration and ICIs First, an enrichment analysis was undertake on a diverse array of 28 immune cell populations, spanning two risk cohorts within the TCGA-THCA training dataset ( Fig. 8 A ) . Subsequently, the Wilcoxon test was employed to discern any discrepancies amongst the 28 immune cell types between the two distinct risk cohorts, revealing differences in 15 immune cells, defined as differential expression immune cells ( Fig. 8 B ) . Following this, Spearman correlation analysis was conducted on the differential expression immune cells, revealing an overall high positive correlation, such as central memory CD4 T cell and effector memeory CD8 T cell (r = 0.9), central memory CD4 T cell and MDSC (r = 0.93) ( Fig. 8 C ) . Further exploration of the correlation amongst the variant expression of immune cells, risk scores, and prognostic genes was conducted using Spearman correlation analysis. It was found that the highest positive and negative correlations were a pair each, namely TSHZ3 with CD56 natural killer cells (cor=-0.44, p = 6.26e − 16 ) and AK1 with CD56 bright natural killer cell (cor = 0.24, p = 1.64e − 05 ) ( Fig. 8 D ) . Finally, through the analysis of the predicted and distributed differences in immune score, stromal score, and estimate score, statistically significant variances in immune scores were observed among the two risk cohorts (p < 0.05), with the latter exhibiting a notably higher score ( Fig. 8 E ) . Currently, ICIs have emerged as a promising approach in cancer therapy. By scrutinizing the distinctions in ICIs between the high- and low-risk cohorts, 10 differential expression immune checkpoints were identified, defined as different ICIs ( Fig. 8 F ) . Lastly, correlation analysis revealed the highest positive correlation between risk scores and BCL2 ( Fig. 8 G ) . 3.9 Functional analysis, lncRNA-miRNA-mRNA regulatory network, and expression validation of prognostic genes We predicted through GeneMANIA that the related genes were all co-expression genes, without predicting their functions ( Fig. 9 A ) . Subsequently, using Friends analysis, we identified the functional similarity among prognostic genes. The highest score was obtained by AKAP12, indicating the strongest functional similarity with other prognostic genes ( Fig. 9 B ) . To elucidate the regulatory mechanisms in which prognostic genes are involved in PTC, we intersected the predicted results from the TarBase v8.0 and miRDB databases to obtain 52 crucial miRNAs. Then, using the miRNet database, we predicted the lncRNAs and circRNAs that regulate these crucial miRNA. Finally, we assembled a complete lncRNA-miRNA-mRNA regulatory network by Cytoscape, which consisted of 141 nodes and 259 edges ( Fig. 9 C ) . During the analysis of prognostic gene expression disparities between PTC cohort and control cohort, we observed heightened expression of AK1 and SERGEF in PTC cohort relative to control cohort, while the other 5 genes showed the opposite pattern ( Fig. 9 D ) . Furthermore, we downloaded the immunohistochemical staining patterns of prognostic genes from the HPA database to validate the differences in protein levels among the PTC cohort and the normal cohort. The findings illuminated discernible variances in the protein levels of prognostic genes between PTC tissues and normal tissues ( Fig. 9 E ) . 3.10 Key genes validation through RT-qPCR We checked the concentration of the extracted RNA and found it to be within the standard range ( Supplementary Table 3 ). Next, we performed RT-qPCR on all samples. By comparing the expression levels of key genes between the PTC group and the control group, we found significant differences in the expression of TSHZ3, SERGEF, AKAP12, SGPP2, and ASGR1 ( Fig. 10 A- 10 E ) , with ASGR1 showing the most significant difference (p < 0.0001). However, the expression differences of AK1 and PELI2 between the two groups were not significant ( Fig. 10 F, 10 G ) , but this may have been due to the limited sample size. Discussion In this study, 1007 DEGs were obtained by GEO database analysis, and 6628 module genes were obtained by ssGSEA and WGCNA after the expansion of 9 parthanatos related genes. 82 intersection genes were obtained by the intersection of DEGs and module genes. Then GO and KEGG enrichment analysis and PPI protein interaction analysis were performed. The biological functions of GO enriched mainly included cellular cation homeostasis, divalent inorganic cation homeostasis, cellular metal ion homeostasis, synapse organization, negative regulation of transport, cell cortex, clathrin-coated vesicle membrane, coated vesicle membrane, clathrin-coated vesicle, platelet dense tubular network, protein kinase A binding, SH3 domain binding, clathrin adaptor activity, cargo adaptor activity and nucleobase-containing compound kinase activity. We performed univariate Cox regression and LASSO regression based on these 82 genes, and finally used 7 genes to construct the risk model, namely TSHZ3, SERGEF, AKAP12, SGPP2, ASGR1, AK1, and PELI2. Among the 7 genes, SERGEF, AKAP12 and AK1 had been reported in thyroid cancer, especially AK1 had been covered in PTC. In an analysis of the gene expression and CNA profiles of Chernobyl Pediatric Patients, AK1 was found only in children who had been exposed to radiation, which means radiation altered AK1 gene expression [ 36 ] . In addition, AK1 is associated with prognosis in renal cell carcinoma [ 37 ] , myeloid leukemia [ 38 ] , and lung adenocarcinoma [ 39 ] . Hugo Gomez-Rueda et al. used the differences in gene expression values and a multivariate search tool coupled with a nearest centroid classifier to identify SERGEF and other gene as the biomarker which could improve the diagnosis of cytologically indeterminate thyroid cancers [ 40 ] . Similar to our study, Yin Tian et al . also used bioinformatics methods to study the potential regulatory mechanism of thyroid cancer prognosis, and ultimately constructed lasso risk factors using 9 genes including AKAP12 [ 41 ] . Although the other prognostic genes have not been reported in thyroid cancer, they had been studied in other cancers. TSHZ3 had been constructed a prognostic model in Bladder Cancer [ 42 ] , ASGR1 was studied in Hepatocellular carcinoma [ 43 ] , PELI2 was studied in myeloma [ 44 ] and gastric cancer [ 45 ] , and SGPP2 was studied in prostate cancer [ 46 ] . The most extensively studied prognostic genes are still used to construct prognostic related models, which is similar to our research. This indicates that using bioinformatics to construct risk models is a direction for predicting disease prognosis. Afterwards, we conducted in-depth explore on prognostic genes and risk scores, including correlation analysis with clinical features. The results of clinical feature manifest that the survival rate of the high-risk cohorts is significant inferior to the low-risk cohorts in aged, gender, stage3, T3, N, and M0. The prognostic genes TSHZ3 and ASGR1 have also been found in the different stages of the studies with other cancers. TSHZ3 and other 5 gene were selected as differentially expressed homeobox genes and were used to construct prognostic model in Bladder cancer [ 42 ] . In the analysis of clinical status, the expression of TSHZ3 was positively correlated with T stage, N stage, and BLCA stage. Yinjiang Zhang et al . found that the expression of ASGR1 was significantly correlated with T stage in Hepatocellular Carcinoma through public database, and the higher expression level exhibited a longer disease-free survival rate [ 47 ] . Although TSHZ3 and ASGR1 have been studied in different types of cancer, these findings suggest that TSHZ3 and ASGR1 have prognostic roles in different clinical stages or TNM stages. Studies has shown that immune system is associated with the occurrence, invasion, and metastasis of thyroid cancer, and this may affect the treatment and prognosis of patients [ 48 ] . In our study, analysis of immune infiltration closely followed analysis of clinical features. 28 types of immune infiltrating cells show differences between high and low-risk cohorts, and the expression of 28 immune cells revealed a positive correlation, which demonstrated that this signature was significantly correlated with immune infiltration. In the following analysis of the correlation between differential immune cells and prognostic genes, we found that CD56 and TSHZ3 showed the highest negative correlation and AK1 showed the highest positive correlation. TSHZ3 and AK1 were associated with immune cells were also found in some non-thyroid cancer studies. In a study of Bladder Cance, there were differences in the infiltration of most immune cells between high-risk cohorts and low-risk cohorts, and TSHZ3 was related to almost all immune cell types [ 42 ] . AK1 was associated with CD8+ T cells, and it might be involved in immune-related mechanisms were coverd in dilated cardiomyopathy with fibrosis [ 49 ] . In addition to the above mentioned contents, we also analyzed immune checkpoint analysis, GSEA functional enrichment analysis, drug sensitivity, functional similarity analysis, the construction of GeneMANIA network and ceRNA network for prognostic genes, and the mRNA and protein expression of prognostic genes between PTC group and normal group. Since the data of this study were from TCGA and GEO, its reliability cannot be guaranteed. In addition, although we verified the mRNA expression levels of prognostic genes, our sample size was insufficient (n = 10) and a large number of additional samples were needed. In order to further study the molecular mechanism of these 7 prognostic genes, we can not only detect their mRNA expression, but also carry out further studies in the laboratory through cell or animal experiments. For the accuracy of the risk model, we had validated with the data from TCGA, but But the amount of data needs to be replenished. In this study, 1007 DEGs were obtained by GEO database analysis, and 6628 module genes were obtained by ssGSEA and WGCNA after the expansion of 9 parthanatos related genes. 82 intersection genes were obtained by the intersection of DEGs and module genes. Then GO and KEGG enrichment analysis and PPI protein interaction analysis were performed. The biological functions of GO enriched mainly included cellular cation homeostasis, divalent inorganic cation homeostasis, cellular metal ion homeostasis, synapse organization, negative regulation of transport, cell cortex, clathrin-coated vesicle membrane, coated vesicle membrane, clathrin-coated vesicle, platelet dense tubular network, protein kinase A binding, SH3 domain binding, clathrin adaptor activity, cargo adaptor activity and nucleobase-containing compound kinase activity. We performed univariate Cox regression and LASSO regression based on these 82 genes, and finally used 7 genes to construct the risk model, namely TSHZ3, SERGEF, AKAP12, SGPP2, ASGR1, AK1, and PELI2. Among the 7 genes, SERGEF, AKAP12 and AK1 had been reported in thyroid cancer, especially AK1 had been covered in PTC. In an analysis of the gene expression and CNA profiles of Chernobyl Pediatric Patients, AK1 was found only in children who had been exposed to radiation, which means radiation altered AK1 gene expression [ 36 ] . In addition, AK1 is associated with prognosis in renal cell carcinoma [ 37 ] , myeloid leukemia [ 38 ] , and lung adenocarcinoma [ 39 ] . Hugo Gomez-Rueda et al. used the differences in gene expression values and a multivariate search tool coupled with a nearest centroid classifier to identify SERGEF and other gene as the biomarker which could improve the diagnosis of cytologically indeterminate thyroid cancers [ 40 ] . Similar to our study, Yin Tian et al . also used bioinformatics methods to study the potential regulatory mechanism of thyroid cancer prognosis, and ultimately constructed lasso risk factors using 9 genes including AKAP12 [ 41 ] . Although the other prognostic genes have not been reported in thyroid cancer, they had been studied in other cancers. TSHZ3 had been constructed a prognostic model in Bladder Cancer [ 42 ] , ASGR1 was studied in Hepatocellular carcinoma [ 43 ] , PELI2 was studied in myeloma [ 44 ] and gastric cancer [ 45 ] , and SGPP2 was studied in prostate cancer [ 46 ] . The most extensively studied prognostic genes are still used to construct prognostic related models, which is similar to our research. This indicates that using bioinformatics to construct risk models is a direction for predicting disease prognosis. Declarations Funding Not applicable. Conflict of Interest All authors declare that they have none of the conflicts of interest. Ethics The protocol was approved by the ethics committee of Third Hospital of Shanxi Medical University(Approval notice number: YXLL-2024-173) in accordance with the Declaration of Helsinki. Informed consents (Consent to Participate and Consent to Publish) were obtained from all participants or, if participants are under 18, from a parent and/or legal guardian. Author contributions Rui Wang and Shuxin Wen conceived and designed the project. Rui Wang and Li Zhang conducted the experiments. Rui Wang analyzed and interpreted the data. Rui Wang, Li Zhang and Shuxin Wen wrote the manuscript. Acknowledgment Not applicable. Clinical trial number Not applicable. Consent to Publish declaration Not applicable. Data Availability Statement The GSE33630 data are available as open data via the Gene Expression Omnibus database：https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE33630. The TCGA-THCA dataset are available as open data via the Cancer Genome Atlas Program database: https://portal.gdc.cancer.gov/projects/TCGA-THCA. References None. Guidelines for the diagnosis and treatment of thyroid cancer (2022)[J]. Chin J Practical Surg. 2022;42(12):16. Ito Y, Miyauchi A, Fujishima M, Masuoka H, Higashiyama T, Kihara M, Onoda N, Miya A. Prognostic significance of patient age in papillary thyroid carcinoma with no high-risk features. Endocr J. 2022;69(9):1131–6. La Vecchia C, Malvezzi M, Bosetti C, Garavello W, Bertuccio P, Levi F, Negri E. Thyroid cancer mortality and incidence: a global overview. Int J Cancer. 2015;136(9):2187–95. Schindler AM, van Melle G, Evequoz B, Scazziga B. Prognostic factors in papillary carcinoma of the thyroid. Cancer. 1991;68(2):324–30. Shaha AR. Implications of prognostic factors and risk groups in the management of differentiated thyroid cancer. Laryngoscope. 2004;114(3):393–402. Wei X, Wang X, Xiong J, Li C, Liao Y, Zhu Y, Mao J. Risk and Prognostic Factors for BRAFV600E Mutations in Papillary Thyroid Carcinoma. Biomed Res Int. 2022;2022:9959649. Xing M. Prognostic utility of BRAF mutation in papillary thyroid cancer. Mol Cell Endocrinol. 2010;321(1):86–93. Chen B, Shi Y, Xu Y, Zhang J. The predictive value of coexisting BRAFV600E and TERT promoter mutations on poor outcomes and high tumour aggressiveness in papillary thyroid carcinoma: A systematic review and meta-analysis. Clin Endocrinol (Oxf). 2021;94(5):731–42. Liu R, Xing M. TERT promoter mutations in thyroid cancer. Endocr Relat Cancer. 2016;23(3):R143–55. Scheffel RS, Maia AL. The long and still uncertain journey of BRAF as a prognostic tool in patients with papillary thyroid cancer. Arch Endocrinol Metab. 2019;63(2):95–6. Nasirden A, Saito T, Fukumura Y, Hara K, Akaike K, Kurisaki-Arakawa A, Asahina M, Yamashita A, Tomomasa R, Hayashi T, Arakawa A, Yao T. In Japanese patients with papillary thyroid carcinoma, TERT promoter mutation is associated with poor prognosis, in contrast to BRAF V600E mutation. Virchows Arch. 2016;469(6):687–96. Huang P, Chen G, Jin W, Mao K, Wan H, He Y. Molecular Mechanisms of Parthanatos and Its Role in Diverse Diseases. Int J Mol Sci. 2022;23(13):7292. Koehler RC, Dawson VL, Dawson TM. Targeting Parthanatos in Ischemic Stroke. Front Neurol. 2021;12:662034. Xu X, Sun B, Zhao C. Poly (ADP-Ribose) polymerase 1 and parthanatos in neurological diseases: From pathogenesis to therapeutic opportunities. Neurobiol Dis. 2023;187:106314. Wang Y, An R, Umanah GK, Park H, Nambiar K, Eacker SM, Kim B, Bao L, Harraz MM, Chang C, Chen R, Wang JE, Kam TI, Jeong JS, Xie Z, Neifert S, Qian J, Andrabi SA, Blackshaw S, Zhu H, Song H, Ming GL, Dawson VL, Dawson TM. A nuclease that mediates cell death induced by DNA damage and poly(ADP-ribose) polymerase-1. Science. 2016;354(6308):aad6872. Wang X, Ge P. Parthanatos in the pathogenesis of nervous system diseases. Neuroscience. 2020;449:241–50. Jiang X, Deng W, Tao S, Tang Z, Chen Y, Tian M, Wang T, Tao C, Li Y, Fang Y, Pu C, Gao J, Wang X, Qu W, Gai X, Ding Z, Fu Y, Zheng Y, Cao S, Zhou J, Huang M, Liu W, Xu J, Fan J, Shi Y. A RIPK3-independent role of MLKL in suppressing parthanatos promotes immune evasion in hepatocellular carcinoma. Cell Discov. 2023;9(1):7. Zou Y, Xie J, Zheng S, Liu W, Tang Y, Tian W, Deng X, Wu L, Zhang Y, Wong CW, Tan D, Liu Q, Xie X. Leveraging diverse cell-death patterns to predict the prognosis and drug sensitivity of triple-negative breast cancer patients after surgery. Int J Surg. 2022;107:106936. Li C, Zhang J, Wu Q, Kumar A, Pan G, Kelvin DJ. Nifuroxazide Activates the Parthanatos to Overcome TMPRSS2:ERG Fusion-Positive Prostate Cancer. Mol Cancer Ther. 2023;22(3):306–16. Li C, Xue Y, Wu J, Zhang L, Yang T, Ai M, Han J, Zheng X, Wang R, Boldogh I, Ba X. MTH1 inhibition synergizes with ROS-inducing agents to trigger cervical cancer cells undergoing parthanatos. Biochim Biophys Acta Mol Basis Dis. 2024;1870(5):167190. Chen DL, Cai JH, Wang CCN. Identification of Key Prognostic Genes of Triple Negative Breast Cancer by LASSO-Based Machine Learning and Bioinformatics Analysis. Genes (Basel). 2022;13(5):902. Wu B, Xi S. Bioinformatics analysis of differentially expressed genes and pathways in the development of cervical cancer. BMC Cancer. 2021;21(1):733. Zhang Y, Shan L, Li D, Tang Y, Qian W, Dai J, Du M, Sun X, Zhu Y, Wang Q, Zhou L. Identification of key biomarkers associated with immune cells infiltration for myocardial injury in dermatomyositis by integrated bioinformatics analysis. Arthritis Res Ther. 2023;25(1):69. Qin H, Abulaiti A, Maimaiti A, Abulaiti Z, Fan G, Aili Y, Ji W, Wang Z, Wang Y. Integrated machine learning survival framework develops a prognostic model based on inter-crosstalk definition of mitochondrial function and cell death patterns in a large multicenter cohort for lower-grade glioma. J Transl Med. 2023;21(1):588. Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W, Smyth GK. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. 2015;43(7):e47. Gu Z, Eils R, Schlesner M. Complex heatmaps reveal patterns and correlations in multidimensional genomic data. Bioinformatics. 2016;32(18):2847–9. Hänzelmann S, Castelo R, Guinney J. GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinformatics. 2013;14:7. Langfelder P, Horvath S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics. 2008;9:559. Wu T, Hu E, Xu S, Chen M, Guo P, Dai Z, Feng T, Zhou L, Tang W, Zhan L, Fu X, Liu S, Bo X, Yu G. clusterProfiler 4.0: A universal enrichment tool for interpreting omics data. Innov (Camb). 2021;2(3):100141. Shannon P, Markiel A, Ozier O, Baliga NS, Wang JT, Ramage D, Amin N, Schwikowski B, Ideker T. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 2003;13(11):2498–504. Friedman J, Hastie T, Tibshirani R. Regularization Paths for Generalized Linear Models via Coordinate Descent. J Stat Softw. 2010;33(1):1–22. Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15(12):550. Fu D, Zhang B, Wu S, Zhang Y, Xie J, Ning W, Jiang H. Prognosis and Characterization of Immune Microenvironment in Acute Myeloid Leukemia Through Identification of an Autophagy-Related Signature. Front Immunol. 2021;12:695865. Mas MT, Chen CY, Hitzeman RA, Riggs AD. Active human-yeast chimeric phosphoglycerate kinases engineered by domain interchange. Science. 1986;233(4765):788–90. Yu G. Gene Ontology Semantic Similarity Analysis Using GOSemSim. Methods Mol Biol. 2020;2117:207–15. Stein L, Rothschild J, Luce J, Cowell JK, Thomas G, Bogdanova TI, Tronko MD, Hawthorn L. Copy number and gene expression alterations in radiation-induced papillary thyroid carcinoma from chernobyl pediatric patients. Thyroid. 2010;20(5):475–87. Liu J, Liu B, Guo Y, Chen Z, Sun W, Gao W, Wu H, Wang Y. Key miRNAs and target genes played roles in the development of clear cell renal cell carcinoma. Cancer Biomark. 2018;23(2):279–90. Qin T, Zhao H, Shao Y, Hu N, Shi J, Fu L, Zhang Y. High expression of AK1 predicts inferior prognosis in acute myeloid leukemia patients undergoing chemotherapy. Biosci Rep. 2020;40(6):BSR20200097. Jan YH, Lai TC, Yang CJ, Huang MS, Hsiao M. A co-expressed gene status of adenylate kinase 1/4 reveals prognostic gene signature associated with prognosis and sensitivity to EGFR targeted therapy in lung adenocarcinoma. Sci Rep. 2019;9(1):12329. Gomez-Rueda H, Palacios-Corona R, Gutiérrez-Hermosillo H, Trevino V. A robust biomarker of differential correlations improves the diagnosis of cytologically indeterminate thyroid cancers. Int J Mol Med. 2016;37(5):1355–62. Tian Y, Xie T, Sun X. Analysis of the regulatory mechanisms of prognostic immune factors in thyroid cancer. Front Oncol. 2022;12:1059591. Dong B, Liang J, Li D, Song W, Song J, Zhu M, Zhao S, Ma Y, Yang T. Identification of a Prognostic Signature Associated With the Homeobox Gene Family for Bladder Cancer. Front Mol Biosci. 2021;8:688298. Roa-Colomo A, López Garrido MÁ, Molina-Vallejo P, Rojas A, Sanchez MG, Aranda-García V, Salmeron J, Romero-Gomez M, Muntane J, Padillo J, Alamo JM, Lorente JA, Serrano MJ, Garrido-Navas MC. Hepatocellular carcinoma risk-stratification based on ASGR1 in circulating epithelial cells for cancer interception. Front Mol Biosci. 2022;9:1074277. Masuda T, Haji S, Nakashima Y, Tsuda M, Kimura D, Takamatsu A, Iwahashi N, Umakoshi H, Shiratsuchi M, Kikutake C, Suyama M, Ohkawa Y, Ogawa Y. Identification of a drug-response gene in multiple myeloma through longitudinal single-cell transcriptome sequencing. iScience. 2022;25(8):104781. Zhang C, Jing LW, Li ZT, Chang ZW, Liu H, Zhang QM, Zhang QY. Identification of a prognostic 28-gene expression signature for gastric cancer with lymphatic metastasis. Biosci Rep. 2019;39(5):BSR20182179. Zhang Y, Kong X, Xin S, Bi L, Sun X. Discovery of Lipid Metabolism-Related Genes for Predicting Tumor Immune Microenvironment Status and Prognosis in Prostate Cancer. J Oncol. 2022;2022:8227806. Zhang Y, Wei H, Fan L, Fang M, He X, Lu B, Pang Z. CLEC4s as Potential Therapeutic Targets in Hepatocellular Carcinoma Microenvironment. Front Cell Dev Biol. 2021;9:681372. Yin H, Tang Y, Guo Y, Wen S. Immune Microenvironment of Thyroid Cancer. J Cancer. 2020;11(16):4884–96. Li W, Liu P, Liu H, Zhang F, Fu Y. Integrative analysis of genes reveals endoplasmic reticulum stress-related immune responses involved in dilated cardiomyopathy with fibrosis. Apoptosis. 2023;28(9–10):1406–21. Supplementary. Additional Declarations No competing interests reported. Supplementary Files SupplementaryTable13.docx Cite Share Download PDF Status: Under Review Version 1 posted Editorial decision: Revision requested 22 Sep, 2025 Reviews received at journal 06 Sep, 2025 Reviews received at journal 12 Aug, 2025 Reviewers agreed at journal 01 Aug, 2025 Reviewers agreed at journal 27 Jul, 2025 Reviewers invited by journal 15 Jul, 2025 Editor assigned by journal 15 Jul, 2025 Editor invited by journal 30 Jun, 2025 Submission checks completed at journal 29 Jun, 2025 First submitted to journal 28 Jun, 2025 You are reading this latest preprint version Research Square lets you share your work early, gain feedback from the community, and start making changes to your manuscript prior to peer review in a journal. As a division of Research Square Company, we’re committed to making research communication faster, fairer, and more useful. We do this by developing innovative software and high quality services for the global research community. Our growing team is made up of researchers and industry professionals working together to solve the most critical problems facing scientific publishing. Also discoverable on Platform About Our Team In Review Editorial Policies Advisory Board Help Center Resources Author Services Accessibility API Access RSS feed Manage Cookie Preferences © Research Square 2026 | ISSN 2693-5015 (online) Privacy Policy Terms of Service Do Not Sell My Personal Information {"props":{"pageProps":{"initialData":{"identity":"rs-6673087","acceptedTermsAndConditions":true,"allowDirectSubmit":false,"archivedVersions":[],"articleType":"Research Article","associatedPublications":[],"authors":[{"id":487028225,"identity":"dcaf976b-19d0-46fb-b7e8-150647981230","order_by":0,"name":"Rui Wang","email":"","orcid":"","institution":"1．Third Hospital of Shanxi Medical University, Shanxi Academy of Medical Sciences, Tongji Shanxi hospital","correspondingAuthor":false,"prefix":"","firstName":"Rui","middleName":"","lastName":"Wang","suffix":""},{"id":487028226,"identity":"63279d99-096f-407f-a077-11680486d281","order_by":1,"name":"Li Zhang","email":"","orcid":"","institution":"1．Third Hospital of Shanxi Medical University, Shanxi Academy of Medical Sciences, Tongji Shanxi hospital","correspondingAuthor":false,"prefix":"","firstName":"Li","middleName":"","lastName":"Zhang","suffix":""},{"id":487028227,"identity":"700b4269-bf41-4828-9c53-b3c574c773ad","order_by":2,"name":"Shuxin Wen","email":"data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAZAAAAAyAQMAAABI0h/eAAAABlBMVEX///8AAABVwtN+AAAACXBIWXMAAA7EAAAOxAGVKw4bAAAAqElEQVRIiWNgGAWjYBADOTb29gOkaTHm4zmTQJqWxHkSDgbEKTU43nv4M8+fw+ltEgwJDD8qthGh5cy5NMmZbWm5bdKNBxh7ztwmrMXsRo4Zw8cGm9w2mQMJzIxtxGi5/8b4Q8IfiXQ2iQQDIrXc4DGQ+MBmk0C8FvszOWYgvxi2AQP5IFF+kWw/YwwKMXn59vaDD35UEKEFBRwgUf0oGAWjYBSMAlwAAPdYOuFedU65AAAAAElFTkSuQmCC","orcid":"","institution":"1．Third Hospital of Shanxi Medical University, Shanxi Academy of Medical Sciences, Tongji Shanxi hospital","correspondingAuthor":true,"prefix":"","firstName":"Shuxin","middleName":"","lastName":"Wen","suffix":""}],"badges":[],"createdAt":"2025-05-15 13:23:28","currentVersionCode":1,"declarations":"","doi":"10.21203/rs.3.rs-6673087/v1","doiUrl":"https://doi.org/10.21203/rs.3.rs-6673087/v1","draftVersion":[],"editorialEvents":[],"editorialNote":"","failedWorkflow":false,"files":[{"id":87323720,"identity":"c762da31-78d1-4491-8279-0175443deec0","added_by":"auto","created_at":"2025-07-22 17:07:03","extension":"jpg","order_by":1,"title":"Figure 1","display":"","copyAsset":false,"role":"figure","size":940146,"visible":true,"origin":"","legend":"\u003cp\u003eIndentification of parthanatos-related differentially expressed genes in PTC(DE-PRGs). A: Volcano plot of DEGs. Expression levels of the DEGs are represented by the different colors, Red, up‑regulated; green, down‑regulated; grey, normal expression. B: The top half showed a heatmap of the expression density of DEGs, the bottom half showed a heatmap of the expression of DEGs. C: \u003cem\u003eKaplan-Meier\u003c/em\u003e curve was established by disease samples of TCGA-THCA. D: Sample clustering showed no outlier samples. Each branch represents a sample, and the vertical axis represents the Euclidean distance of sample expression. E: Sample clustering and heatmap of sample traits. F: Determination of soft threshold. G: Identification of co-expression modules. The top half showed the hierarchical clustering tree of genes, and the bottom half showed the gene module. H: The correlation between module eigengenes and traits was shown using heat map. I: 82 DE-PRGs were identified based on the Venn diagram. DEGs, differentially expressed genes; PTC, papillary thyroid carcinoma.\u003c/p\u003e","description":"","filename":"Fig1.jpg","url":"https://assets-eu.researchsquare.com/files/rs-6673087/v1/48b214b00cdf613f12e27047.jpg"},{"id":87322502,"identity":"0baed3f9-c2b4-46cd-8db6-53c846111986","added_by":"auto","created_at":"2025-07-22 16:51:03","extension":"jpg","order_by":2,"title":"Figure 2","display":"","copyAsset":false,"role":"figure","size":587143,"visible":true,"origin":"","legend":"\u003cp\u003eFunctional enrichment and protein interaction analysis of DE-PRGS. A: GO enrichment analysis enriched 177 GO biological functions. The number above the blue layer indicates the number of genes enriched. The higher the inner column, the more significant. B: KEGG enrichment analysis enriched 8 KEGG pathways. The number above the blue layer indicates the number of genes enriched. The higher the inner column, the more significant. C: PPI network was constructed based on 82 DE-PRGs. The stronger the interaction, the darker the color. GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; PPI, protein-protein interaction.\u003c/p\u003e","description":"","filename":"Fig2.jpg","url":"https://assets-eu.researchsquare.com/files/rs-6673087/v1/14db63d6931df4642b5b0880.jpg"},{"id":87323395,"identity":"8fb9fbe7-7161-48d8-b51f-328beee7d364","added_by":"auto","created_at":"2025-07-22 16:59:03","extension":"jpg","order_by":3,"title":"Figure 3","display":"","copyAsset":false,"role":"figure","size":709562,"visible":true,"origin":"","legend":"\u003cp\u003eConstruction of risk model and identification of prognostic genes. A: The forest figure for Univariate Cox regression analysis of the 10 survival-related genes. B-C: LASSO regression analysis of the 10 survival-related genes. D: Risk score distribution of patients. The dotted line is the median risk score and the corresponding number of patients. E: Overall survival (days) distribution of patients. The dotted line is the median risk score and the corresponding number of patients. F: \u003cem\u003eKaplan-Meier\u003c/em\u003e curve of high- and low- risk cohorts. G: ROC survival curves at 3, 5 and 7 years. H: Heatmap of prognostic genes expression in high- and low- risk cohorts. I: The correlation of prognostic genes. LASSO, Least absolute shrinkage and selection operator; ROC, receiver operating characteristic.\u003c/p\u003e","description":"","filename":"Fig3.jpg","url":"https://assets-eu.researchsquare.com/files/rs-6673087/v1/3b2c7bcb1d602e8911fa465c.jpg"},{"id":87321115,"identity":"47f6a7df-dc30-4696-a8a9-0e8a96d1fe37","added_by":"auto","created_at":"2025-07-22 16:35:03","extension":"jpg","order_by":4,"title":"Figure 4","display":"","copyAsset":false,"role":"figure","size":939134,"visible":true,"origin":"","legend":"\u003cp\u003eValidation of the risk model and correlation of porognostic genes. A: Risk score and overall survival (days) distribution of patients. The dotted line is the median risk score and the corresponding number of patients. B: \u003cem\u003eKaplan-Meier\u003c/em\u003e curve of high- and low- risk cohorts. C: ROC survival curves at 3, 5 and 7 years. D: Risk score and overall survival (days) distribution of patients. The dotted line is the median risk score and the corresponding number of patients. E: \u003cem\u003eKaplan-Meier\u003c/em\u003e curve of high- and low- risk cohorts. F: ROC survival curves at 3, 5 and 7 years. G: Heatmap of prognostic genes expression in high- and low-risk cohorts. H: The correlation of prognostic genes. I: Heatmap of prognostic genes expression in high- and low-risk cohorts. J: The correlation of prognostic genes. ROC, receiver operating characteristic.\u003c/p\u003e","description":"","filename":"Fig4.jpg","url":"https://assets-eu.researchsquare.com/files/rs-6673087/v1/1a5332f49d49c77fe050308d.jpg"},{"id":87321124,"identity":"343ec92e-9bd6-4bc4-ac6d-d1bb58466f66","added_by":"auto","created_at":"2025-07-22 16:35:03","extension":"jpg","order_by":5,"title":"Figure 5","display":"","copyAsset":false,"role":"figure","size":446657,"visible":true,"origin":"","legend":"\u003cp\u003eEstablishment and verification of nomogram. A: Nomogram model was constructed with survival as the outcome event by prognostic genes. B: The calibration curve of 3/5/7 year. C: ROC survival curves at 3, 5 and 7 years. ROC, receiver operating characteristic.\u003c/p\u003e","description":"","filename":"Fig5.jpg","url":"https://assets-eu.researchsquare.com/files/rs-6673087/v1/62462cfbdb9efa14574915ba.jpg"},{"id":87323721,"identity":"f63fdde4-16f4-473b-9824-c60065944c28","added_by":"auto","created_at":"2025-07-22 17:07:03","extension":"jpg","order_by":6,"title":"Figure 6","display":"","copyAsset":false,"role":"figure","size":433918,"visible":true,"origin":"","legend":"\u003cp\u003eClinical characterization analysis in high- and low-risk cohorts. A: Age. B: Gender. C: Pathological stage. D: T stage. E: N stage. F: M stage.\u003c/p\u003e","description":"","filename":"Fig6.jpg","url":"https://assets-eu.researchsquare.com/files/rs-6673087/v1/8481c73fbeaa1ae44e8a58dd.jpg"},{"id":87321121,"identity":"63757ea0-3930-4231-9ff1-f66aec5f232b","added_by":"auto","created_at":"2025-07-22 16:35:03","extension":"jpg","order_by":7,"title":"Figure 7","display":"","copyAsset":false,"role":"figure","size":226163,"visible":true,"origin":"","legend":"\u003cp\u003eGSEA enrichment analyses and drug sensitivity analyses. A: GSEA enrichment in high- and low- risk cohorts. B: drug sensitivity analyses in high- and low- risk cohorts. GSEA, Gene Set Enrichment Analysis.\u003c/p\u003e","description":"","filename":"Fig7.jpg","url":"https://assets-eu.researchsquare.com/files/rs-6673087/v1/745a96c8966b4c5d23137110.jpg"},{"id":87322123,"identity":"82b87a82-0090-4a6e-beeb-1def17d97df8","added_by":"auto","created_at":"2025-07-22 16:43:03","extension":"jpg","order_by":8,"title":"Figure 8","display":"","copyAsset":false,"role":"figure","size":993228,"visible":true,"origin":"","legend":"\u003cp\u003eAnalysis of immune infiltration and immune-checkpoint inhibitors. A: Heat map of 28 immune cells expressed in high- and low- risk cohort. B: Differences of 15 types of immune cells in high- and low- risk cohort. C: Heat map in differential immune cell correlation. D: Heat map between prognostic genes and differential immune cell correlation. E: ESTIMATE score in high- and low- risk cohort. F: Differences in the expression of ICIs in high- and low- risk cohorts. G: Heat map between risk score and differential ICIs correlation. ICIs, immune-checkpoint inhibitors.\u003c/p\u003e","description":"","filename":"Fig8.jpg","url":"https://assets-eu.researchsquare.com/files/rs-6673087/v1/c4e802f7fd3831a1680460bd.jpg"},{"id":87323722,"identity":"03b57d72-7149-4dd1-bb60-120a1e828220","added_by":"auto","created_at":"2025-07-22 17:07:03","extension":"jpg","order_by":9,"title":"Figure 9","display":"","copyAsset":false,"role":"figure","size":808374,"visible":true,"origin":"","legend":"\u003cp\u003eFunctional analysis, lncRNA-miRNA-mRNA regulatory network, and expression validation of prognostic genes. A: Gene and function related to prognostic genes function were predicted by GeneMANIA. B: Friends analyses of prognostic genes. C: CeRNA network of prognostic genes. D: Expression of prognostic genes. E: The protein levels of prognostic genes in PTC cohort and normal control cohort.\u003c/p\u003e","description":"","filename":"Fig9.jpg","url":"https://assets-eu.researchsquare.com/files/rs-6673087/v1/9ba32b8315672e7aecea1a48.jpg"},{"id":87322126,"identity":"b064e0fd-eaeb-43bd-9ce5-f631eb81fcb5","added_by":"auto","created_at":"2025-07-22 16:43:03","extension":"jpg","order_by":10,"title":"Figure 10","display":"","copyAsset":false,"role":"figure","size":317816,"visible":true,"origin":"","legend":"\u003cp\u003eThe expression of 7 prognostic genes genes. A: ASGR1. B: AKAP12. C: SERGEF. D: TSHZ3. E: SGPP2. F: PELI2. G: AK1.\u003c/p\u003e","description":"","filename":"Fig10.jpg","url":"https://assets-eu.researchsquare.com/files/rs-6673087/v1/d61496eac9634dce239c5a89.jpg"},{"id":87467095,"identity":"0d236bf2-baf1-4377-851e-2bdcbb1dde28","added_by":"auto","created_at":"2025-07-24 07:59:04","extension":"pdf","order_by":0,"title":"","display":"","copyAsset":false,"role":"manuscript-pdf","size":7814861,"visible":true,"origin":"","legend":"","description":"","filename":"manuscript.pdf","url":"https://assets-eu.researchsquare.com/files/rs-6673087/v1/10b62c11-249e-44af-9c43-7094fe94b900.pdf"},{"id":87323397,"identity":"7a391208-b7f5-4c85-a6c1-a1b019b84076","added_by":"auto","created_at":"2025-07-22 16:59:03","extension":"docx","order_by":1,"title":"","display":"","copyAsset":false,"role":"supplement","size":17119,"visible":true,"origin":"","legend":"","description":"","filename":"SupplementaryTable13.docx","url":"https://assets-eu.researchsquare.com/files/rs-6673087/v1/c3b18ffec122cdeea8275fdc.docx"}],"financialInterests":"No competing interests reported.","formattedTitle":"Construction and validation of risk models of prognostic genes associated with parthanatos in papillary thyroid carcinoma based on bioinformatics","fulltext":[{"header":"Introduction","content":"\u003cp\u003eThyroid cancer is a malignant tumor originating from the follicular epithelium or parafollicular epithelial cells of the thyroid gland, and is the most common malignant tumor in the head and neck. In recent years, the incidence of thyroid cancer in the world has increased rapidly, and thyroid cancer in China has continued to increase at a rate of 20% per year\u003csup\u003e[\u003cspan citationid=\"CR1\" class=\"CitationRef\"\u003e1\u003c/span\u003e]\u003c/sup\u003e. According to the origin and differentiation of tumors, thyroid cancer is further divided into: Papillary Thyroid Carcinoma (PTC), Follicular Thyroid Carcinoma (FTC), Medullary Thyroid Carcinoma (MTC), and poorly differentiated thyroid carcinoma (PDTC) and anaplastic thyroid cancer (ATC), of which PTC is the most common, accounting for about 90% of the entire thyroid cancer. PTC generally has a good prognosis, however, distant metastasis and poor prognosis can also be found in some aggressive subtypes of PTC\u003csup\u003e[\u003cspan citationid=\"CR2\" class=\"CitationRef\"\u003e2\u003c/span\u003e\u0026ndash;\u003cspan citationid=\"CR3\" class=\"CitationRef\"\u003e3\u003c/span\u003e]\u003c/sup\u003e. Therefore, it is of great significance to study the development mechanism of PTC for the treatment of thyroid cancer. The prognostic of PTC is associated with age, sex, loss of histologic differentiation, tumor size, presence of carcinomatous lymphangitis, extrathyroid extension, and presence of metastasis\u003csup\u003e[\u003cspan citationid=\"CR4\" class=\"CitationRef\"\u003e4\u003c/span\u003e\u0026ndash;\u003cspan citationid=\"CR5\" class=\"CitationRef\"\u003e5\u003c/span\u003e]\u003c/sup\u003e. The prognostic prediction studies of PTC are mostly conducted by molecular markers, among which BRAF\u003csup\u003e[\u003cspan citationid=\"CR6\" class=\"CitationRef\"\u003e6\u003c/span\u003e\u0026ndash;\u003cspan citationid=\"CR7\" class=\"CitationRef\"\u003e7\u003c/span\u003e]\u003c/sup\u003e and TERT\u003csup\u003e[\u003cspan citationid=\"CR8\" class=\"CitationRef\"\u003e8\u003c/span\u003e\u0026ndash;\u003cspan citationid=\"CR9\" class=\"CitationRef\"\u003e9\u003c/span\u003e]\u003c/sup\u003e are the main studies. However, the prognostic prediction of BRAF is still controversial\u003csup\u003e[\u003cspan citationid=\"CR10\" class=\"CitationRef\"\u003e10\u003c/span\u003e\u0026ndash;\u003cspan citationid=\"CR11\" class=\"CitationRef\"\u003e11\u003c/span\u003e]\u003c/sup\u003e. Therefore, the exploration of other prognostic genes and their molecular mechanisms provides a new idea for prognosis assessment and individualized treatment of PTC patients.\u003c/p\u003e\u003cp\u003eThe biological homeostasis of an organism depends on the balance between cell proliferation and death. Parthanatos is a poly(ADP-ribose) polymerase 1 (PARP1) hyperactivation dependent and cysteine-independent cell death which is distinct from apoptosis, programmed necrosis, or other known forms of cell death\u003csup\u003e[\u003cspan citationid=\"CR12\" class=\"CitationRef\"\u003e12\u003c/span\u003e]\u003c/sup\u003e. Parthanatos is a multi-step pathway that plays a crucial role in tumorigenesis. Once the DNA is damaged, it leads to excessive activation of PARP, then PARP-1 produces free PAR by PARG-mediated hydrolyzation, inducing the release of apoptosis-inducing factor(AIF) from the outer mitochondrial membrane\u003csup\u003e[\u003cspan additionalcitationids=\"CR13\" citationid=\"CR12\" class=\"CitationRef\"\u003e12\u003c/span\u003e\u0026ndash;\u003cspan citationid=\"CR14\" class=\"CitationRef\"\u003e14\u003c/span\u003e]\u003c/sup\u003e. Once AIF is released, it can bind to macrophage migration inhibitory factor (MIF), translocates to the nucleus, and eventually lead to a large number of DNA fragments\u003csup\u003e[\u003cspan citationid=\"CR15\" class=\"CitationRef\"\u003e15\u003c/span\u003e\u0026ndash;\u003cspan citationid=\"CR16\" class=\"CitationRef\"\u003e16\u003c/span\u003e]\u003c/sup\u003e. Parthanatos has been studied in hepatocellular carcinoma\u003csup\u003e[\u003cspan citationid=\"CR17\" class=\"CitationRef\"\u003e17\u003c/span\u003e]\u003c/sup\u003e, breast cancer\u003csup\u003e[\u003cspan citationid=\"CR18\" class=\"CitationRef\"\u003e18\u003c/span\u003e]\u003c/sup\u003e, prostate cancer\u003csup\u003e[\u003cspan citationid=\"CR19\" class=\"CitationRef\"\u003e19\u003c/span\u003e]\u003c/sup\u003e, cervical cancer\u003csup\u003e[\u003cspan citationid=\"CR20\" class=\"CitationRef\"\u003e20\u003c/span\u003e]\u003c/sup\u003e and other cancers, however, the function and mechanism of parthanatos in PTC remain unclear and require further extensive research.\u003c/p\u003e\u003cp\u003eIn recent years, it has become increasingly mature to screen key genes and study their molecular mechanisms through bioinformatics\u003csup\u003e[\u003cspan additionalcitationids=\"CR22\" citationid=\"CR21\" class=\"CitationRef\"\u003e21\u003c/span\u003e\u0026ndash;\u003cspan citationid=\"CR23\" class=\"CitationRef\"\u003e23\u003c/span\u003e]\u003c/sup\u003e. In this study, prognostic genes associated with parthanatos in PTC was identified with bioinformatics through two public databases, then, their potential molecular mechanisms was been explored, and a prognostic risk model was constructed to predict the survival rate of PTC patients, providing effective targets for personalized treatment and prognostic intervention for PTC patients.\u003c/p\u003e"},{"header":"Materials and methods","content":"\u003cdiv id=\"Sec3\" class=\"Section2\"\u003e\u003ch2\u003e2.1 Dataset selection and description\u003c/h2\u003e\u003cp\u003eWe obtained genomic expression profiles and clinical feature data (including age, gender, subtype, grade, stage, lymphovascular invasion, lymph nodes positive by he, etc.) of the TCGA-THCA dataset from the Cancer Genome Atlas Program database (TCGA, \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://portal.gdc.cancer.gov/\u003c/span\u003e\u003cspan address=\"https://portal.gdc.cancer.gov/\" targettype=\"URL\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e), comprising 510 thyroid carcinoma (THCA) tissues and 58 normal tissue samples, which were utilized for prognostic analysis and validation. The samples from the TCGA-THCA dataset were apportioned randomly into two subsets in a ratio of 6 to 4, with the larger segment designated as the training set and the remainder allocated for validation purposes. The chip data of papillary thyroid carcinoma (PTC) patients and normal thyroid tissue samples (49/45) from the GSE33630 (GPL570 platform) dataset were procured from the Gene Expression Omnibus database (GEO, \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://www.ncbi.nlm.nih.gov/geo/\u003c/span\u003e\u003cspan address=\"https://www.ncbi.nlm.nih.gov/geo/\" targettype=\"URL\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e) database for differential expression analysis. Additionally, 9 parthanatos-related genes (PRGs) were identified from prior research\u003csup\u003e[\u003cspan citationid=\"CR24\" class=\"CitationRef\"\u003e24\u003c/span\u003e]\u003c/sup\u003e.\u003c/p\u003e\u003c/div\u003e\n\u003ch3\u003e2.2 Differential expression analysis\u003c/h3\u003e\n\u003cp\u003eTo delineate the landscape of differentially expressed genes (DEGs) in PTC patients relative to the control cohort, the expression matrix derived from the GSE33630 dataset was meticulously analyzed using the \"limma \" package (version 3.54.0)\u003csup\u003e[\u003cspan citationid=\"CR25\" class=\"CitationRef\"\u003e25\u003c/span\u003e]\u003c/sup\u003e. This analytical sieve employed a stringent criterion of adj.p.value\u0026thinsp;\u0026lt;\u0026thinsp;0.05 \u0026amp; |log\u003csub\u003e2\u003c/sub\u003eFC| \u0026gt;1 to discern DEGs. To gain an overall understanding of the distribution of DEGs, the \"ggplot2\" (version\u0026thinsp;=\u0026thinsp;3.4.1) package was used to draw volcano plots for visualizing the differential genes and the \"ComplexHeatmap\" (version 2.14.0)\u003csup\u003e[\u003cspan citationid=\"CR26\" class=\"CitationRef\"\u003e26\u003c/span\u003e]\u003c/sup\u003e package was utilized to generate a heatmap illustrative of the top 10 genes demonstrating significant up-regulation and the top 10 genes exhibiting marked down-regulation, meticulously sorted by the respective fold changes in their expression levels.\u003c/p\u003e\n\u003ch3\u003e2.3 Weighted correlation network analysis (WGCNA)\u003c/h3\u003e\n\u003cp\u003eIn the TCGA-THCA disease samples, ssGSEA scores were calculated based on 9 PRGs using the \"GSVA\" (version 1.46.0)\u003csup\u003e[\u003cspan citationid=\"CR27\" class=\"CitationRef\"\u003e27\u003c/span\u003e]\u003c/sup\u003e package. Tumor samples were categorically segregated into high and low ssGSEA score cohorts, demarcated by an optimal threshold value derived from the ssGSEA scoring metric. Survival disparities between the high and low ssGSEA score were delineated through the construction of Kaplan-Meier (KM) survival curves. The difference in KM curves is considered significant when p\u0026thinsp;\u0026lt;\u0026thinsp;0.02.\u003c/p\u003e\u003cp\u003eThen, the expression matrix of disease samples in TCGA-THCA was analyzed utilizing WGCNA through the implementation of the \"WGCNA\" (version 1.71)\u003csup\u003e[\u003cspan citationid=\"CR28\" class=\"CitationRef\"\u003e28\u003c/span\u003e]\u003c/sup\u003e package. To commence the analytical process, an initial hierarchical clustering was performed on the expression data from all samples within the TCGA-THCA dataset, employing the Euclidean distance metric to assess the similarity of expression levels. This step was crucial for identifying and excluding any outliers among the samples, thereby refining the dataset for subsequent analysis. Subsequently, to ensure the constructed network adhered to a scale-free topology, an optimal soft threshold needed to be determined before constructing the co-expression network and set to r2\u0026thinsp;=\u0026thinsp;0.85. Ultimately, an unscaled network was crafted utilizing the chosen soft threshold parameter, effectively segmenting the genes into distinct modules, each adorned with a unique color identifier. Using the parthanatos related ssGSEA scores as the phenotype, associations between phenotypes and modules were established. This involved calculating the correlation coefficient between modules and phenotype traits, followed by generating a heatmap of correlations. The most positively and negatively associated modules with the parthanatos related ssGSEA scores (|cor| \u0026gt;0.3, p\u0026thinsp;\u0026lt;\u0026thinsp;0.05) were identified in this study as key modules.\u003c/p\u003e\u003cp\u003eUtilizing the \"ggVenn\" package (version 0.1.9), an overlap analysis was conducted between the DEGs and the key module genes, resulting in the identification of differential genes associated with parthanatos in PTC, referred to as DE-PRGs.\u003c/p\u003e\n\u003ch3\u003e2.4 Gene function enrichment and protein-protein interaction (PPI) analysis\u003c/h3\u003e\n\u003cp\u003eWe employed the \"clusterProfiler\" (version 4.7.1.3)\u003csup\u003e[\u003cspan citationid=\"CR29\" class=\"CitationRef\"\u003e29\u003c/span\u003e]\u003c/sup\u003e package to perform Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses on the DE-PRGs. These analyses were employed with a stringent cutoff of p-values\u0026thinsp;\u0026lt;\u0026thinsp;0.05 to curate a concise list of significantly enriched biological processes and pathways.\u003c/p\u003e\u003cp\u003eTo establish a PPI network for DE-PRGs, genetic entities were submitted to the STRING database (\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttp://www.string-db.org/\u003c/span\u003e\u003cspan address=\"http://www.string-db.org/\" targettype=\"URL\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e) to get results (interaction score\u0026thinsp;\u0026gt;\u0026thinsp;0.15). Subsequently, the outcomes were imported into Cytoscape (version 3.7.2)\u003csup\u003e[\u003cspan citationid=\"CR30\" class=\"CitationRef\"\u003e30\u003c/span\u003e]\u003c/sup\u003e for graphical representation and analysis.\u003c/p\u003e\n\u003ch3\u003e2.5 Construction of risk models\u003c/h3\u003e\n\u003cp\u003eIn order to identify the genes associated with the prognosis of PTC patients, a univariate Cox regression analysis was performed on DE-PRGs within the TCGA-THCA training dataset using the \"survival\" (version 3.5-3) package. The data used consisted of both gene expression information of THCA samples and training set data from TCGA-THCA. Prior to conducting univariate Cox regression analysis, the assumption of proportional hazards (PH) for the risk model results was tested using the R software's cox.zph function (p\u0026thinsp;\u0026gt;\u0026thinsp;0.05) for DE-PRGs. After excluding those that failed the PH test, the remaining DE-PRGs were used to construct univariate Cox regression model (p\u0026thinsp;\u0026lt;\u0026thinsp;0.05 \u0026amp; HR\u0026thinsp;\u0026ne;\u0026thinsp;1). High HR values were excluded to avoid influencing subsequent analysis, and the remaining survival genes were then integrated into our subsequent analytical procedures. The Forestplot (version 2.0.1) package was employed to generate a graphical forest plot, elucidating the outcomes of the univariate Cox regression analysis in a visually intuitive manner.\u003c/p\u003e\u003cp\u003eThereafter, the remaining genes underwent a LASSO regression analytical process utilizing the \"glmnet\" (version 4.1-6)\u003csup\u003e[\u003cspan citationid=\"CR31\" class=\"CitationRef\"\u003e31\u003c/span\u003e]\u003c/sup\u003e package to delineate prognostic genes and to establish a risk prediction model. The analysis was conducted with the \"family\" parameter specified as Cox, thereby incorporating the Cox proportional hazards model, and prioritizing strongly correlated features for inclusion. The model was considered optimal when the model error was minimized, and the corresponding lambda value at this point was referred to as lambdamin.\u003c/p\u003e\u003cdiv id=\"Sec8\" class=\"Section2\"\u003e\u003ch2\u003e2.6 Evaluation of risk models\u003c/h2\u003e\u003cp\u003eIn the TCGA-THCA training set, individual patients with PTC were assigned a score derived from the expression of prognostic genes and the risk coefficients ascertained through LASSO regression analysis. This segmentation process resulted in the classification of patients into distinct high- and low-risk cohorts, determined by an empirically derived optimal cutoff for the risk scoring mechanism. The computation of the risk score was achieved through a linear combination of the gene expression levels and their respective weights within the predictive framework. The formulaic representation is as follows:\u003c/p\u003e\u003cp\u003e\u003cspan class=\"InlineEquation\"\u003e\u003c/span\u003e\u003c/p\u003e\u003cp\u003eWhere, \"risk score\" represents the risk score, \"coef\" signifying the risk coefficients attributed to each specific gene, and \"expr\" represents the expression of the respective genes. The \"survminer\" (version 0.4.9) package was employed to delineate the survival curves for overall survival (OS) among the two risk cohorts. Furthermore, to gauge the precision of risk prognostication, receiver operating characteristic (ROC) curves were generated employing the \"survivalROC\" (version 0.4.9) package based on TCGA-THCA training set data, with survival time points at 3/5/7 years.\u003c/p\u003e\u003cp\u003eThe distinct characteristics separating the two risk cohorts were ascertained by assessing the transcriptomic levels of prognostic genes was analyzed based on the threshold criteria |cor|\u0026gt;0.3 and p\u0026lt;0.05.\u003c/p\u003e\u003c/div\u003e\n\u003ch3\u003e2.7 Validation of the risk model\u003c/h3\u003e\n\u003cp\u003eTo gauge the efficacy of the risk model, it underwent rigorous validation within the TCGA-THCA validation set and TCGA-THCA dataset. The risk model derived from the aforementioned training cohort was applied, and using the optimal cutoff value of risk score, to stratify the PTC patients within the TCGA-THCA validation subset and TCGA-THCA dataset into distinct high- and low-risk cohorts. Subsequently, the \"survminer\" package was deployed to illustrate the survival and ROC curves for overall survival, thereby comparing the prognostic outcomes between the two risk cohorts. Then, the differences between the two risk cohorts in the validation set were observed drawing on the expression patterns. Simultaneously, the \"psych\" (version\u0026thinsp;=\u0026thinsp;2.2.9) package was employed to calculate the Spearman correlation coefficient between prognostic genes (|cor|\u0026gt;0.3 and p\u0026lt;0.05).\u003c/p\u003e\n\u003ch3\u003e2.8 Construction and evaluation of nomogram\u003c/h3\u003e\n\u003cp\u003eBuilding upon the foundational analysis of prognostic gene expression, survival was considered as the outcome event. A predictive nomogram was then artfully constructed through the \"rms\" (version 6.5.0) package. For the purposes of validation, calibration curves were meticulously plotted with the \"regplot\" (version 1.1), while ROC curves were delineated by the \"survivalROC\" (AUC\u0026thinsp;\u0026lt;\u0026thinsp;0.6 ), thereby assessing the nomogram's clinical utility in forecasting the progression of PTC.\u003c/p\u003e\u003cdiv id=\"Sec11\" class=\"Section2\"\u003e\u003ch2\u003e2.9 Correlation analysis of survival differences and clinical features\u003c/h2\u003e\u003cp\u003eWithin the TCGA-THCA training set, the Wilcoxon test was employed to elucidate the disparities in survival outcomes among the high- and low-risk cohorts across various subgroups categorized by distinct clinical feature (including age, gender, pathological stage, T stage, N stage, and M stage ), with a threshold set at p\u0026thinsp;\u0026lt;\u0026thinsp;0.05. K-M curves were plotted using the \"survminer\" package for visualization.\u003c/p\u003e\u003c/div\u003e\u003cdiv id=\"Sec12\" class=\"Section2\"\u003e\u003ch2\u003e2.10 Gene set enrichment analysis (GSEA)\u003c/h2\u003e\u003cp\u003eA differential expression analysis was performed on the respective cohort within the TCGA-THCA training set, utilizing the \"DESeq2\" (version 1.38.0)\u003csup\u003e[\u003cspan citationid=\"CR32\" class=\"CitationRef\"\u003e32\u003c/span\u003e]\u003c/sup\u003e package, to delineate the biological pathways underlying the high- and low-risk cohorts. According to the ranking of log\u003csub\u003e2\u003c/sub\u003eFC, the \"clusterProfiler\" package was employed to conduct GSEA utilizing the gene set \u0026ldquo;c2.cp.kegg.v2023.1.Hs.symbols.gm\u0026rdquo; derived from the latest version of the Molecular Signatures Database (MSigDB, \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://www.gsea-msigdb.org/gsea/msigdb/index.jsp)(|NES|\u003c/span\u003e\u003cspan address=\"https://www.gsea-msigdb.org/gsea/msigdb/index.jsp)(|NES|\" targettype=\"URL\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e \u0026gt;1 and p\u0026thinsp;\u0026lt;\u0026thinsp;0.05).\u003c/p\u003e\u003c/div\u003e\u003cdiv id=\"Sec13\" class=\"Section2\"\u003e\u003ch2\u003e2.11 Analysis of the immune microenvironment\u003c/h2\u003e\u003cp\u003eIn the TCGA-THCA training set, the PTC patients underwent immune infiltration analysis utilizing the ssGSEA algorithm, exploring the immune cell infiltration status in PTC patients. Firstly, we derived the enrichment scores for 28 immune cells in the different risk cohorts. Subsequently, the Wilcoxon test was utilized to ascertain the presence of statistically significant disparities in the distribution of immune cell populations between the two risk cohorts within the TCGA-THCA training dataset. Then, for all samples in the training set of TCGA-THCA, the \"corrplot\" (0.92) was utilized for visualize the correlation matrix to perform Spearman correlation analysis on differential immune cells. Simultaneously, Spearman correlation analysis was conduct on differential immune cells, risk scores, and prognostic genes, considering the results to be significant if |cor| \u0026gt;0.3 and p\u0026thinsp;\u0026lt;\u0026thinsp;0.05. Ultimately, the quantification of stromal and immune cells in tumor tissue was achieved through the use of the \"ESTIMATE\" algorithm, thereby enabling the prediction of the immune score, stromal score, and ESTIMATE score. Additionally, it facilitated the estimation of tumor purity, as well as the levels of immune and stromal infiltration within the malignant tumor tissue. Differences in the distribution of three scores were analyzed between two risk cohorts using the Wilcoxon test. To discern the disparities in the distribution of the three scores between the high- and low-risk cohorts, the Wilcoxon test was employed.\u003c/p\u003e\u003c/div\u003e\u003cdiv id=\"Sec14\" class=\"Section2\"\u003e\u003ch2\u003e2.12 Analysis of immune checkpoint inhibitors (ICIs)\u003c/h2\u003e\u003cp\u003eWe obtained 14 ICIs genes (ASXL1, BCL2, CD33, CD47, CHEK1, PLK1, CTLA4, DOT1L, FLT3, IDH1, IDH2, MCL1, PDCD1, MDM2) from the pervoius studies \u003csup\u003e[\u003cspan citationid=\"CR33\" class=\"CitationRef\"\u003e33\u003c/span\u003e\u0026ndash;\u003cspan citationid=\"CR34\" class=\"CitationRef\"\u003e34\u003c/span\u003e]\u003c/sup\u003e. Within the TCGA-THCA training dataset, the Wilcoxon test was invoked to delineate the disparities in ICIs between the high- and low-risk cohorts. Those ICIs exhibiting statistically significant variances (p\u0026thinsp;\u0026lt;\u0026thinsp;0.05) across the two risk cohorts were chosen for further exploration. Following this, the \"cor\" function within the R programming environment was utilized to elucidate the Spearman correlation analysis, elucidating the relationship among the risk score and the differential expression of ICIs.\u003c/p\u003e\u003c/div\u003e\u003cdiv id=\"Sec15\" class=\"Section2\"\u003e\u003ch2\u003e2.13 Drug sensitivity analysis\u003c/h2\u003e\u003cp\u003eTo ascertain the disparities in drug responsiveness between the high- and low-risk cohorts, the Genomics of Drug Sensitivity in Cancer (GDSC, \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://www.cancerrxgene.org/\u003c/span\u003e\u003cspan address=\"https://www.cancerrxgene.org/\" targettype=\"URL\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e) database furnished the half-maximal inhibitory concentrations (IC\u003csub\u003e50\u003c/sub\u003e) and PTC drug. The \"pRRophetic\" (version 0.5)\u003csup\u003e[\u003cspan citationid=\"CR25\" class=\"CitationRef\"\u003e25\u003c/span\u003e]\u003c/sup\u003e package was deployed to forecast the IC\u003csub\u003e50\u003c/sub\u003e for each tumor specimen within the TCGA-THCA training cohort. Subsequently, the divergences in drug sensitivity between the two risk cohorts were subjected to Wilcoxon test (p\u0026thinsp;\u0026lt;\u0026thinsp;0.05).\u003c/p\u003e\u003c/div\u003e\u003cdiv id=\"Sec16\" class=\"Section2\"\u003e\u003ch2\u003e2.14 Functional analysis of prognostic genes\u003c/h2\u003e\u003cp\u003eIn order to discern other genes associated with prognostic gene function, the prediction of genes and gene function associated with prognosis was achieved through the use of the the GeneMANIA database (\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttp://genemania.org\u003c/span\u003e\u003cspan address=\"http://genemania.org\" targettype=\"URL\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e). The functional similarity between prognosis genes was examined employing the \"GOSemSim\" (version 2.24.0)\u003csup\u003e[\u003cspan citationid=\"CR35\" class=\"CitationRef\"\u003e35\u003c/span\u003e]\u003c/sup\u003e package, with a score greater than 0.5 being determined as indicating high relevance.\u003c/p\u003e\u003c/div\u003e\u003cdiv id=\"Sec17\" class=\"Section2\"\u003e\u003ch2\u003e2.15 Construction of the lncRNA-miRNA-mRNA network for prognostic genes\u003c/h2\u003e\u003cp\u003eThe miRNAs that regulated prognostic genes were forecasted based on the TarBase v8.0(\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttp://mirtarbase.mbc.nctu.edu.tw/\u003c/span\u003e\u003cspan address=\"http://mirtarbase.mbc.nctu.edu.tw/\" targettype=\"URL\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e) database and miRDB database (\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://mirdb.org/\u003c/span\u003e\u003cspan address=\"https://mirdb.org/\" targettype=\"URL\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e) using the Network Analyst platform. The \"ggVenn\" (version 0.1.9) package was utilized to intersect the miRNAs predicted from the two databases, obtaining key miRNAs. Subsequently, the miRNet database (\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://www.mirnet.ca/miRNet/home.xhtml\u003c/span\u003e\u003cspan address=\"https://www.mirnet.ca/miRNet/home.xhtml\" targettype=\"URL\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e) was employed to forecast the lncRNAs and circRNAs of the key miRNAs. The top 3 lncRNAs and circRNAs were selected as the regulators of the miRNAs. The Cytoscape (version 3.7.2)\u003csup\u003e[\u003cspan citationid=\"CR30\" class=\"CitationRef\"\u003e30\u003c/span\u003e]\u003c/sup\u003e platform was utilized to build a comprehensive lncRNA-miRNA-mRNA regulatory network.\u003c/p\u003e\u003c/div\u003e\u003cdiv id=\"Sec18\" class=\"Section2\"\u003e\u003ch2\u003e2.16 Validation of gene expression for prognostic genes\u003c/h2\u003e\u003cp\u003eThe expression differences of prognostic genes in PTC cohort and normal cohort were analyzed employing the Wilcoxon test based on the entire sample of the TCGA-THCA dataset (p\u0026thinsp;\u0026lt;\u0026thinsp;0.05), to further corroborate the expression of prognostic genes. Finally, immunohistochemical staining maps of prognostic genes were procured from the Human Protein Atlas database (HPA, \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttp://www.proteinalas.org/\u003c/span\u003e\u003cspan address=\"http://www.proteinalas.org/\" targettype=\"URL\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e) to visualize the differences in protein levels of prognostic genes between the PTC cohort and normal cohort.\u003c/p\u003e\u003c/div\u003e\u003cdiv id=\"Sec19\" class=\"Section2\"\u003e\u003ch2\u003e2.17 Validation of reverse transcription-quantitative real-time Polymerase chain reaction (RT-qPCR)\u003c/h2\u003e\u003cp\u003eThe expression of genes was validated via RT-qPCR. The tissue samples were obtained from Third Hospital of Shanxi Medical University, and it has been approved by the ethics committee (Approval notice number: YXLL-2024-173) and all the participants have signed the informed consent. The results were plotted using Graphpad Prism 10.\u003c/p\u003e\u003cp\u003eThe TRizol reagent kit was used to extract RNA from 5 pairs of frozen tissue samples, with samples 1\u0026ndash;5 designated as the Control group and samples 6\u0026ndash;10 as the PTC group. All experimental procedures were carried out according to the instructions, resulting in total RNA. 1 \u0026micro;l of the extracted RNA was taken for concentration detection using the Nano Photometer N50, and the purity/concentration was recorded to calculate the amount for subsequent reverse transcription steps. The remaining RNA was either immediately subjected to reverse transcription or stored at -80\u0026deg;C. The SweScript First Strand cDNA Synthesis SuperMix for PCR from Servicebio Company was used to reverse transcribe the RNA into cDNA, following the instructions. The cDNA was diluted 5\u0026ndash;20 times with ddH\u003csub\u003e2\u003c/sub\u003eO (RNase/ARase free), with 3ul of cDNA, 5ul of 2xUniversal Blue SYBR Green qPCR Master Mix, 1ul of forward primer (10\u0026micro;M), 1ul of reverse primer (10\u0026micro;M). Subsequently, 40 cycles of reaction were carried out using the CFX96 real-time quantitative PCR instrument, with the program detailed in \u003cb\u003eSupplementary Table\u0026nbsp;1\u003c/b\u003e. The primer sequences for the TSHZ3, SERGEF, AKAP12, ASGR1, AK1, PELI2 and SGPP2 genes can be found in \u003cb\u003eSupplementary Table\u0026nbsp;2\u003c/b\u003e, and GAPDH was used as the reference gene. The 2-ΔΔCT method was used to calculate the relative gene expression levels.\u003c/p\u003e\u003c/div\u003e\u003cdiv id=\"Sec20\" class=\"Section2\"\u003e\u003ch2\u003e2.18 Statistical analysis\u003c/h2\u003e\u003cp\u003eAll analyses were conducted using the R programming language (version 4.2.3). The p-values and r-values for correlation analysis were calculated using non-parametric Spearman correlation coefficient. The disparities between groups were scrutinized utilizing the Wilcoxon test, with statistical significance being attributed to findings where p\u0026thinsp;\u0026lt;\u0026thinsp;0.005.\u003c/p\u003e\u003c/div\u003e"},{"header":"Results","content":"\u003cdiv id=\"Sec22\" class=\"Section2\"\u003e\n \u003ch2\u003e3.1 Indentification of 82 DE-PRGs\u003c/h2\u003e\n \u003cp\u003eThrough differential expression analysis of the PTC and the control cohort, 1,007 DEGs were identified (557 up-regulated genes; 450 down-regulated genes) \u003cstrong\u003e(\u003c/strong\u003eFig. \u003cspan class=\"InternalRef\"\u003e1\u003c/span\u003eA\u003cstrong\u003e)\u003c/strong\u003e. Additionally, the top 10 up- and down-regulated genes were displayed based on log\u003csub\u003e2\u003c/sub\u003eFC sorting \u003cstrong\u003e(\u003c/strong\u003eFig. \u003cspan class=\"InternalRef\"\u003e1\u003c/span\u003eB\u003cstrong\u003e)\u003c/strong\u003e.\u003c/p\u003e\n \u003cp\u003eIn the disease samples of TCGA-THCA, the parthanatos scores of the disease cohort were computed for the 9 PRGs obtained from previous studies. The tumor specimens were categorically segregated into distinct cohorts, distinguished by high and low score brackets, based on the optimally determined cutoff threshold derived from the ssGSEA scoring metrics. The significant difference in the high- and low-risk cohorts was observed through the K-M curve (p\u0026thinsp;\u0026lt;\u0026thinsp;0.02) (Fig. \u003cspan class=\"InternalRef\"\u003e1\u003c/span\u003eC), indicating that it was meaningful for us to use DE-PRGs genes as the background gene set for scoring. Subsequently, WGCNA was undertaken utilizing the expression matrix of disease samples within the TCGA-THCA dataset. Firstly, all samples were clustered, and no outlier samples were found \u003cstrong\u003e(\u003c/strong\u003eFig. \u003cspan class=\"InternalRef\"\u003e1\u003c/span\u003eD\u003cstrong\u003e)\u003c/strong\u003e. Then, sample traits (ssGSEA scores) were introduced, and hierarchical clustering was performed again using the Euclidean distance of sample expression levels \u003cstrong\u003e(\u003c/strong\u003eFig. \u003cspan class=\"InternalRef\"\u003e1\u003c/span\u003eE\u003cstrong\u003e)\u003c/strong\u003e. In the pursuit of constructing a scale-free gene co-expression network, the key step of finding the optimal soft threshold was carried out before building the co-expression network. Upon achieving a scale-free fit index of 0.85, the soft threshold with a connectivity close to 0 and a \u0026beta; value of 10 is chosen for screening \u003cstrong\u003e(\u003c/strong\u003eFig. \u003cspan class=\"InternalRef\"\u003e1\u003c/span\u003eF\u003cstrong\u003e)\u003c/strong\u003e.\u003c/p\u003e\n \u003cp\u003eUtilizing the chosen soft threshold, the scale-free network was constructed, and these genes were divided into several modules, each marked with a different color. A total of 11 co-expression modules were identified through WGCNA analysis \u003cstrong\u003e(\u003c/strong\u003eFig. \u003cspan class=\"InternalRef\"\u003e1\u003c/span\u003eG\u003cstrong\u003e)\u003c/strong\u003e (excluding the grey module to which unclassified genes belong). Using the parthanatos related ssGSEA score as a phenotype, the module eigengenes and the correlation coefficient matrix between phenotype and traits were obtained \u003cstrong\u003e(\u003c/strong\u003eFig. \u003cspan class=\"InternalRef\"\u003e1\u003c/span\u003eH\u003cstrong\u003e)\u003c/strong\u003e. The blue module (cor\u0026thinsp;=\u0026thinsp;0.51,p\u0026thinsp;=\u0026thinsp;1e\u003csup\u003e\u0026minus;\u0026thinsp;35\u003c/sup\u003e) and the turquoise module (cor=-0.38, p\u0026thinsp;=\u0026thinsp;9e\u003csup\u003e\u0026minus;\u0026thinsp;19\u003c/sup\u003e) with the highest significant correlation with the parthanatos related ssGSEA score, denoted as the key modules, were chosen. By selecting all module genes, 6,628 module genes were obtained. Finally, the obtained DEGs were intersected with key module genes, resulting in 82 intersecting genes, referred to as DE-PRGs \u003cstrong\u003e(\u003c/strong\u003eFig. \u003cspan class=\"InternalRef\"\u003e1\u003c/span\u003eI\u003cstrong\u003e)\u003c/strong\u003e.\u003c/p\u003e\n \u003cdiv id=\"Sec23\" class=\"Section3\"\u003e\n \u003ch2\u003e3.2 Functional enrichment and protein interaction analysis of DE-PRGs\u003c/h2\u003e\n \u003cp\u003eIn order to explore the biological functions and signaling pathways, we conducted GO and KEGG enrichment analyses on the 82 DE-PRGs. In the GO enrichment analysis, a total of 177 GO biological functions were enriched (p-value\u0026thinsp;\u0026lt;\u0026thinsp;0.05 )(144 in biological processes (BPs), 16 in molecular functions (MFs), and 17 in cellular components (CCs)) \u003cstrong\u003e(\u003c/strong\u003eFig. \u003cspan class=\"InternalRef\"\u003e2\u003c/span\u003eA\u003cstrong\u003e)\u003c/strong\u003e. The enriched GO biological functions mainly included cellular cation homeostasis, divalent inorganic cation homeostasis, and cellular metal ion homeostasis. In the KEGG enrichment analysis, based on the screening criterion of p-value\u0026thinsp;\u0026lt;\u0026thinsp;0.05, a total of 8 KEGG pathways were enriched \u003cstrong\u003e(\u003c/strong\u003eFig. \u003cspan class=\"InternalRef\"\u003e2\u003c/span\u003eB\u003cstrong\u003e)\u003c/strong\u003e. The enriched pathways included regulation of actin cytoskeleton, calcium signaling pathway, thyroid hormone synthesis, prostate cancer, complement and coagulation cascades, virion - ebolavirus and lyssavirus, sulfur metabolism, and thiamine metabolism.\u003c/p\u003e\n \u003cp\u003eSimultaneously, to investigate the protein-level interaction of the DE-PRGs, we deposited the 82 DE-PRGs into the STRING to forge a DE-PRGs protein interaction network (interaction score\u0026thinsp;\u0026gt;\u0026thinsp;0.15) \u003cstrong\u003e(\u003c/strong\u003eFig. \u003cspan class=\"InternalRef\"\u003e2\u003c/span\u003eC\u003cstrong\u003e).\u003c/strong\u003e In the PPI network, after removing isolated proteins, a network of interactions among 73 proteins was obtained, containing 73 nodes and 117 edges. Among them, AUTS2 exhibited the strongest interaction.\u003c/p\u003e\n \u003c/div\u003e\n\u003c/div\u003e\n\u003cdiv id=\"Sec24\" class=\"Section2\"\u003e\n \u003ch2\u003e3.3 Construction of risk model and identification of prognostic genes\u003c/h2\u003e\n \u003cp\u003eIn order to discern genes linked to prognosis in PTC patients the univariate Cox regression was conducted on 82 DE-PRGs within the TCGA-THCA training set. Initially, four DE-PRGs that did not pass the PH assumption test (p\u0026thinsp;\u0026gt;\u0026thinsp;0.05) were excluded, and then the univariate Cox regression model was performed, uncovering a relationship among the expression of 11 genes and OS, with MIPOL1 exhibiting an exceptionally high HR value, potentially influenced by confounding factors and genetic mutations, and was therefore excluded to avoid biasing subsequent analyses. The remaining 10 genes associated with survival were included for further analysis, revealing that among these 10 survival-related genes, SERGEF and AK1 were protective factors, while the others were risk factors \u003cstrong\u003e(\u003c/strong\u003eFig. \u003cspan class=\"InternalRef\"\u003e3\u003c/span\u003eA, Table \u003cspan class=\"InternalRef\"\u003e3\u003c/span\u003e\u0026thinsp;\u003cstrong\u003e\u0026minus;\u0026thinsp;1)\u003c/strong\u003e. Therefore, this study constructed a risk model through LASSO regression analysis, which includes seven genes: TSHZ3, SERGEF, AKAP12, ASGR1, AK1, PELI2 and SGPP2. These genes, were defined as prognostic genes \u003cstrong\u003e(\u003c/strong\u003eFig. \u003cspan class=\"InternalRef\"\u003e3\u003c/span\u003eB, \u003cspan class=\"InternalRef\"\u003e3\u003c/span\u003eC, Table \u003cspan class=\"InternalRef\"\u003e3\u003c/span\u003e\u0026thinsp;\u003cstrong\u003e\u0026minus;\u0026thinsp;2)\u003c/strong\u003e. When lambdamin was set to 0.007, the corresponding number of genes was 7.\u003c/p\u003e\n \u003cdiv class=\"gridtable\"\u003e\u003c/div\u003e\n \u003cdiv class=\"gridtable\"\u003e\n \u003cdiv align=\"left\" class=\"colspec\"\u003e\u003cbr\u003e\u003c/div\u003e\n \u003c/div\u003e\n \u003cdiv class=\"gridtable\"\u003e\n \u003cdiv align=\"char\" class=\"colspec\"\u003e\u003cbr\u003e\u003c/div\u003e\n \u003c/div\u003e\n \u003cp\u003eUtilizing the transcriptome data of the identified prognostic genes in conjunction with the risk coefficients ascertained through LASSO regression analysis, a customized risk evaluation was allocated for every PTC sample. The algorithmic formula for computing this risk score is as follows:\u003c/p\u003e\n \u003cp\u003eriskscore\u0026thinsp;=\u0026thinsp;TSHZ3 * (1.099)\u0026thinsp;+\u0026thinsp;SERGEF * (-0.450)\u0026thinsp;+\u0026thinsp;AKAP12 * (0.204)\u0026thinsp;+\u0026thinsp;SGPP2 * (0.089)\u0026thinsp;+\u0026thinsp;ASGR1 * (0.211)\u0026thinsp;+\u0026thinsp;AK1 * (-0.005)\u0026thinsp;+\u0026thinsp;PELI2 * (0.510)\u003c/p\u003e\n \u003cp\u003eUsing the optimal cutoff value of risk score, PTC cohort were segregated into two distinct groups: a high-risk cohort and a low-risk cohort \u003cstrong\u003e(\u003c/strong\u003eFig. \u003cspan class=\"InternalRef\"\u003e3\u003c/span\u003eD, \u003cspan class=\"InternalRef\"\u003e3\u003c/span\u003eE\u003cstrong\u003e)\u003c/strong\u003e, revealing a notable contrast in survival between the two risk cohorts (p\u0026thinsp;=\u0026thinsp;0.0001), with lower survival rates among the high-risk cohort patients \u003cstrong\u003e(\u003c/strong\u003eFig. \u003cspan class=\"InternalRef\"\u003e3\u003c/span\u003eF\u003cstrong\u003e)\u003c/strong\u003e.\u003c/p\u003e\n \u003cp\u003eFurthermore, the precision of the risk model was assessed utilizing TCGA-THCA training set data. The AUC of the ROC curve consistently exceeded 0.6 in the training set, signifying commendable performance of the risk model \u003cstrong\u003e(\u003c/strong\u003eFig. \u003cspan class=\"InternalRef\"\u003e3\u003c/span\u003eG\u003cstrong\u003e)\u003c/strong\u003e. In terms of the expression levels of prognostic genes, no notable disparities were discerned between the two cohorts \u003cstrong\u003e(\u003c/strong\u003eFig. \u003cspan class=\"InternalRef\"\u003e3\u003c/span\u003eH\u003cstrong\u003e)\u003c/strong\u003e. However, ASGR1 and PELI2 showed lower expression levels in the training set, while SERGEF and AK1 exhibited higher expression levels, with SGPP2 and SERGEF showing the highest negative correlation \u003cstrong\u003e(\u003c/strong\u003eFig. \u003cspan class=\"InternalRef\"\u003e3\u003c/span\u003eI\u003cstrong\u003e)\u003c/strong\u003e.\u003c/p\u003e\n \u003cdiv class=\"gridtable\"\u003e\n \u003ctable id=\"Tab4\" border=\"1\"\u003e\n \u003ccaption language=\"En\"\u003e\n \u003cdiv class=\"CaptionNumber\"\u003eTable 3\u003c/div\u003e\n \u003cdiv class=\"CaptionContent\"\u003e\n \u003cp\u003e\u0026thinsp;\u003cstrong\u003e\u0026minus;\u0026thinsp;1\u003c/strong\u003e. PH test of 10 prognostic genes\u003c/p\u003e\n \u003c/div\u003e\n \u003c/caption\u003e\n \u003cthead\u003e\n \u003ctr\u003e\n \u003cth align=\"left\"\u003e\n \u003cp\u003eID\u003c/p\u003e\n \u003c/th\u003e\n \u003cth align=\"left\"\u003e\n \u003cp\u003eAKAP12\u003c/p\u003e\n \u003c/th\u003e\n \u003cth align=\"left\"\u003e\n \u003cp\u003eTSHZ3\u003c/p\u003e\n \u003c/th\u003e\n \u003cth align=\"left\"\u003e\n \u003cp\u003eF2R\u003c/p\u003e\n \u003c/th\u003e\n \u003cth align=\"left\"\u003e\n \u003cp\u003eASGR1\u003c/p\u003e\n \u003c/th\u003e\n \u003cth align=\"left\"\u003e\n \u003cp\u003ePELI2\u003c/p\u003e\n \u003c/th\u003e\n \u003cth align=\"left\"\u003e\n \u003cp\u003eAK1\u003c/p\u003e\n \u003c/th\u003e\n \u003cth align=\"left\"\u003e\n \u003cp\u003eSGPP2\u003c/p\u003e\n \u003c/th\u003e\n \u003cth align=\"left\"\u003e\n \u003cp\u003eCRYBG3\u003c/p\u003e\n \u003c/th\u003e\n \u003cth align=\"left\"\u003e\n \u003cp\u003eN4BP2\u003c/p\u003e\n \u003c/th\u003e\n \u003cth align=\"left\"\u003e\n \u003cp\u003eSERGEF\u003c/p\u003e\n \u003c/th\u003e\n \u003c/tr\u003e\n \u003c/thead\u003e\n \u003ctbody\u003e\n \u003ctr\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003ePH\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e0.911\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e0.708\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e0.467\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e0.642\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e0.795\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e0.281\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e0.401\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e0.975\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e0.955\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd align=\"left\"\u003e\n \u003cp\u003e0.499\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003c/tbody\u003e\n \u003c/table\u003e\n \u003c/div\u003e\n \u003cdiv class=\"gridtable\"\u003e\n \u003cdiv class=\"colspec\"\u003e\u003cbr\u003e\u003c/div\u003e\n \u003ctable id=\"Tab5\" border=\"1\"\u003e\n \u003ccaption language=\"En\"\u003e\n \u003cdiv class=\"CaptionNumber\"\u003eTable 3\u003c/div\u003e\n \u003cdiv class=\"CaptionContent\"\u003e\n \u003cp\u003e\u0026thinsp;\u003cstrong\u003e\u0026minus;\u0026thinsp;2\u003c/strong\u003e. PRisk coefficient of prognostic genes\u003c/p\u003e\n \u003c/div\u003e\n \u003c/caption\u003e\n \u003ctbody\u003e\n \u003ctr\u003e\n \u003ctd style=\"width: 79px;\"\u003e\n \u003cp\u003eGene\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd style=\"width: 110px;\"\u003e\n \u003cp\u003eCoefficient\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd style=\"width: 79px;\"\u003e\n \u003cp\u003eTSHZ3\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd style=\"width: 110px;\"\u003e\n \u003cp\u003e1.099478866\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd style=\"width: 79px;\"\u003e\n \u003cp\u003eSERGEF\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd style=\"width: 110px;\"\u003e\n \u003cp\u003e-0.449517532\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd style=\"width: 79px;\"\u003e\n \u003cp\u003eAKAP12\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd style=\"width: 110px;\"\u003e\n \u003cp\u003e0.20384777\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd style=\"width: 79px;\"\u003e\n \u003cp\u003eSGPP2\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd style=\"width: 110px;\"\u003e\n \u003cp\u003e0.088714538\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd style=\"width: 79px;\"\u003e\n \u003cp\u003eASGR1\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd style=\"width: 110px;\"\u003e\n \u003cp\u003e0.210755007\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd style=\"width: 79px;\"\u003e\n \u003cp\u003eAK1\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd style=\"width: 110px;\"\u003e\n \u003cp\u003e-0.004941586\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd style=\"width: 79px;\"\u003e\n \u003cp\u003ePELI2\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd style=\"width: 110px;\"\u003e\n \u003cp\u003e0.509768893\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003c/tbody\u003e\n \u003c/table\u003e\n \u003c/div\u003e\n \u003cdiv id=\"Sec25\" class=\"Section3\"\u003e\n \u003ch2\u003e3.4 Validation of the risk assessment model and correlation of prognostic genes\u003c/h2\u003e\n \u003cp\u003eWithin the TCGA-THCA validation set, a marked discrepancy in patient survivability was discerned between the high-risk and low-risk cohorts (p\u0026thinsp;=\u0026thinsp;0.012), with lower survival rates observed in the high risk cohort \u003cstrong\u003e(\u003c/strong\u003eFig. \u003cspan class=\"InternalRef\"\u003e4\u003c/span\u003eA, \u003cspan class=\"InternalRef\"\u003e4\u003c/span\u003eB\u003cstrong\u003e)\u003c/strong\u003e. The AUC values at 3/5/7 years in the ROC curve were all greater than 0.6, indicating excellent validation results \u003cstrong\u003e(\u003c/strong\u003eFig. \u003cspan class=\"InternalRef\"\u003e4\u003c/span\u003eC\u003cstrong\u003e)\u003c/strong\u003e. Similarly, in the TCGA-THCA dataset, a notable contrast in patient survival emerged between the two risk cohorts (p\u0026thinsp;=\u0026thinsp;0.0001), with lower survival rates observed in the high-risk cohort \u003cstrong\u003e(\u003c/strong\u003eFig. \u003cspan class=\"InternalRef\"\u003e4\u003c/span\u003eD, \u003cspan class=\"InternalRef\"\u003e4\u003c/span\u003eE\u003cstrong\u003e)\u003c/strong\u003e. The AUC values at 3/5/7 years in the ROC curve were all greater than 0.6, indicating excellent validation results \u003cstrong\u003e(\u003c/strong\u003eFig. \u003cspan class=\"InternalRef\"\u003e4\u003c/span\u003eF\u003cstrong\u003e)\u003c/strong\u003e.\u003c/p\u003e\n \u003cp\u003eFurthermore, through the observation of the disparate expression of prognostic genes between the two risk cohorts in both datasets, notable distinctions in the expression levels were discerned. In the TCGA-THCA validation set, ASGR1, SGPP2, PELI2 and TSHZ3 showed higher expression in the high-risk cohort contrast to the low-risk cohort \u003cstrong\u003e(\u003c/strong\u003eFig. \u003cspan class=\"InternalRef\"\u003e4\u003c/span\u003eG\u003cstrong\u003e)\u003c/strong\u003e; while in the TCGA-THCA dataset, AK1 and SERGEF showed lower expression in the high-risk cohort contrast to the low-risk cohort \u003cstrong\u003e(\u003c/strong\u003eFig. \u003cspan class=\"InternalRef\"\u003e4\u003c/span\u003eH\u003cstrong\u003e)\u003c/strong\u003e. Additionally, analyzing the correlation between prognostic genes revealed that SERGEF and AK1 (cor\u0026thinsp;=\u0026thinsp;0.53 and P\u0026thinsp;\u0026lt;\u0026thinsp;0.001) exhibited the highest positive correlation in both datasets \u003cstrong\u003e(\u003c/strong\u003eFig. \u003cspan class=\"InternalRef\"\u003e4\u003c/span\u003eI, \u003cspan class=\"InternalRef\"\u003e4\u003c/span\u003eJ\u003cstrong\u003e)\u003c/strong\u003e.\u003c/p\u003e\n \u003c/div\u003e\n \u003cdiv id=\"Sec26\" class=\"Section3\"\u003e\n \u003ch2\u003e3.5 Establishment and verification of nomogram\u003c/h2\u003e\n \u003cp\u003eUtilizing the expression levels of prognostic genes, a nomogram model was constructed with survival as the outcome event \u003cstrong\u003e(\u003c/strong\u003eFig. \u003cspan class=\"InternalRef\"\u003e5\u003c/span\u003eA\u003cstrong\u003e)\u003c/strong\u003e. For the nomogram predicting survival risk, the individual points in the graph represents the score corresponding to the different values of each variable. The total points was derived from the sum of the individual points corresponding to all variable values. The higher the total points, the higher the patient\u0026apos;s survival rate. Then, the survival probabilities for 3, 5, and 7 years were estimated for PTC patients using the total points. The calibration curve (p\u0026thinsp;\u0026gt;\u0026thinsp;0.05, the closer the slope is to 1, the better the accuracy of the model\u0026apos;s predictions) demonstrated that the nomogram model exhibited commendable prognostic power in predicting the clinical outcomes of PTC \u003cstrong\u003e(\u003c/strong\u003eFig. \u003cspan class=\"InternalRef\"\u003e5\u003c/span\u003eB\u003cstrong\u003e)\u003c/strong\u003e. The results from the ROC curve (AUC\u0026thinsp;\u0026gt;\u0026thinsp;0.6 is considered to indicate model accuracy) revealed that the nomogram model exhibited commendable precision in forecasting the clinical outcomes of PTC \u003cstrong\u003e(\u003c/strong\u003eFig. \u003cspan class=\"InternalRef\"\u003e5\u003c/span\u003eC\u003cstrong\u003e)\u003c/strong\u003e. Therefore, it was believed that the nomogram model is effective in predicting PTC clinical outcomes.\u003c/p\u003e\n \u003c/div\u003e\n \u003cdiv id=\"Sec27\" class=\"Section3\"\u003e\n \u003ch2\u003e3.6 Clinical characterization analysis in high- and low-risk cohorts\u003c/h2\u003e\n \u003cp\u003eBased on the TCGA-THCA training set, the survival differences between two risk cohorts were demonstrated in different clinical feature subgroups, including age, gender, pathological stage, T stage, N stage, and M stage (p\u0026thinsp;\u0026lt;\u0026thinsp;0.05). It was discernible that there existed significant disparities in survival outcomes among two risk cohorts for patients aged over 47, gender, stage3, T3, N, and M0, with too few samples for M1 resulting in no outcome \u003cstrong\u003e(\u003c/strong\u003eFig. \u003cspan class=\"InternalRef\"\u003e6\u003c/span\u003eA-\u003cspan class=\"InternalRef\"\u003e6\u003c/span\u003eF\u003cstrong\u003e).\u003c/strong\u003e It was manifest that the high-risk cohorts exhibit a survival rate that is inferior to that observed among the low-risk cohorts.\u003c/p\u003e\n \u003c/div\u003e\n\u003c/div\u003e\n\u003cdiv id=\"Sec28\" class=\"Section2\"\u003e\n \u003ch2\u003e3.7 Enrichment of biological pathways and drug sensitivity analyses\u003c/h2\u003e\n \u003cp\u003eBy conducting GSEA on the two risk cohorts, 29 pathways were ultimately enriched (|NES| \u0026gt;1, p\u0026thinsp;\u0026lt;\u0026thinsp;0.05), with main enrichments in various signaling pathways such as ribosome, focal adhesion, arachidonic acid metabolism, ECM receptor interaction and oxidative phosphorylation \u003cstrong\u003e(\u003c/strong\u003eFig. \u003cspan class=\"InternalRef\"\u003e7\u003c/span\u003eA\u003cstrong\u003e)\u003c/strong\u003e.\u003c/p\u003e\n \u003cp\u003eThrough the analysis of differences in drug sensitivity between high- and low-risk cohorts, we found that 111 drugs showed significant differences, such as cisplatin, EHT.1864, embelin, obatoclax. mesylate and VX.680 \u003cstrong\u003e(\u003c/strong\u003eFig. \u003cspan class=\"InternalRef\"\u003e7\u003c/span\u003eB\u003cstrong\u003e)\u003c/strong\u003e. The top10 drugs (ranking by p-value) were selected to display the differences. And we found these drugs all had higher IC50 in the low risk cohort contrast to the high risk cohort.\u003c/p\u003e\n\u003c/div\u003e\n\u003cdiv id=\"Sec29\" class=\"Section2\"\u003e\n \u003ch2\u003e3.8 Analysis of immune infiltration and ICIs\u003c/h2\u003e\n \u003cp\u003eFirst, an enrichment analysis was undertake on a diverse array of 28 immune cell populations, spanning two risk cohorts within the TCGA-THCA training dataset \u003cstrong\u003e(\u003c/strong\u003eFig. \u003cspan class=\"InternalRef\"\u003e8\u003c/span\u003eA\u003cstrong\u003e)\u003c/strong\u003e. Subsequently, the Wilcoxon test was employed to discern any discrepancies amongst the 28 immune cell types between the two distinct risk cohorts, revealing differences in 15 immune cells, defined as differential expression immune cells \u003cstrong\u003e(\u003c/strong\u003eFig. \u003cspan class=\"InternalRef\"\u003e8\u003c/span\u003eB\u003cstrong\u003e)\u003c/strong\u003e. Following this, Spearman correlation analysis was conducted on the differential expression immune cells, revealing an overall high positive correlation, such as central memory CD4 T cell and effector memeory CD8 T cell (r\u0026thinsp;=\u0026thinsp;0.9), central memory CD4 T cell and MDSC (r\u0026thinsp;=\u0026thinsp;0.93) \u003cstrong\u003e(\u003c/strong\u003eFig. \u003cspan class=\"InternalRef\"\u003e8\u003c/span\u003eC\u003cstrong\u003e)\u003c/strong\u003e. Further exploration of the correlation amongst the variant expression of immune cells, risk scores, and prognostic genes was conducted using Spearman correlation analysis. It was found that the highest positive and negative correlations were a pair each, namely TSHZ3 with CD56 natural killer cells (cor=-0.44, p\u0026thinsp;=\u0026thinsp;6.26e\u003csup\u003e\u0026minus;\u0026thinsp;16\u003c/sup\u003e) and AK1 with CD56 bright natural killer cell (cor\u0026thinsp;=\u0026thinsp;0.24, p\u0026thinsp;=\u0026thinsp;1.64e\u003csup\u003e\u0026minus;\u0026thinsp;05\u003c/sup\u003e) \u003cstrong\u003e(\u003c/strong\u003eFig. \u003cspan class=\"InternalRef\"\u003e8\u003c/span\u003eD\u003cstrong\u003e)\u003c/strong\u003e. Finally, through the analysis of the predicted and distributed differences in immune score, stromal score, and estimate score, statistically significant variances in immune scores were observed among the two risk cohorts (p\u0026thinsp;\u0026lt;\u0026thinsp;0.05), with the latter exhibiting a notably higher score \u003cstrong\u003e(\u003c/strong\u003eFig. \u003cspan class=\"InternalRef\"\u003e8\u003c/span\u003eE\u003cstrong\u003e)\u003c/strong\u003e.\u003c/p\u003e\n \u003cp\u003eCurrently, ICIs have emerged as a promising approach in cancer therapy. By scrutinizing the distinctions in ICIs between the high- and low-risk cohorts, 10 differential expression immune checkpoints were identified, defined as different ICIs \u003cstrong\u003e(\u003c/strong\u003eFig. \u003cspan class=\"InternalRef\"\u003e8\u003c/span\u003eF\u003cstrong\u003e)\u003c/strong\u003e. Lastly, correlation analysis revealed the highest positive correlation between risk scores and BCL2 \u003cstrong\u003e(\u003c/strong\u003eFig. \u003cspan class=\"InternalRef\"\u003e8\u003c/span\u003eG\u003cstrong\u003e)\u003c/strong\u003e.\u003c/p\u003e\n\u003c/div\u003e\n\u003ch3\u003e3.9 Functional analysis, lncRNA-miRNA-mRNA regulatory network, and expression validation of prognostic genes\u003c/h3\u003e\n\u003cp\u003eWe predicted through GeneMANIA that the related genes were all co-expression genes, without predicting their functions \u003cstrong\u003e(\u003c/strong\u003eFig. \u003cspan class=\"InternalRef\"\u003e9\u003c/span\u003eA\u003cstrong\u003e)\u003c/strong\u003e. Subsequently, using Friends analysis, we identified the functional similarity among prognostic genes. The highest score was obtained by AKAP12, indicating the strongest functional similarity with other prognostic genes \u003cstrong\u003e(\u003c/strong\u003eFig. \u003cspan class=\"InternalRef\"\u003e9\u003c/span\u003eB\u003cstrong\u003e)\u003c/strong\u003e.\u003c/p\u003e\n\u003cp\u003eTo elucidate the regulatory mechanisms in which prognostic genes are involved in PTC, we intersected the predicted results from the TarBase v8.0 and miRDB databases to obtain 52 crucial miRNAs. Then, using the miRNet database, we predicted the lncRNAs and circRNAs that regulate these crucial miRNA. Finally, we assembled a complete lncRNA-miRNA-mRNA regulatory network by Cytoscape, which consisted of 141 nodes and 259 edges \u003cstrong\u003e(\u003c/strong\u003eFig. \u003cspan class=\"InternalRef\"\u003e9\u003c/span\u003eC\u003cstrong\u003e)\u003c/strong\u003e.\u003c/p\u003e\n\u003cp\u003eDuring the analysis of prognostic gene expression disparities between PTC cohort and control cohort, we observed heightened expression of AK1 and SERGEF in PTC cohort relative to control cohort, while the other 5 genes showed the opposite pattern \u003cstrong\u003e(\u003c/strong\u003eFig. \u003cspan class=\"InternalRef\"\u003e9\u003c/span\u003eD\u003cstrong\u003e)\u003c/strong\u003e. Furthermore, we downloaded the immunohistochemical staining patterns of prognostic genes from the HPA database to validate the differences in protein levels among the PTC cohort and the normal cohort. The findings illuminated discernible variances in the protein levels of prognostic genes between PTC tissues and normal tissues \u003cstrong\u003e(\u003c/strong\u003eFig. \u003cspan class=\"InternalRef\"\u003e9\u003c/span\u003eE\u003cstrong\u003e)\u003c/strong\u003e.\u003c/p\u003e\n\u003cdiv id=\"Sec31\" class=\"Section2\"\u003e\n \u003ch2\u003e3.10 Key genes validation through RT-qPCR\u003c/h2\u003e\n \u003cp\u003eWe checked the concentration of the extracted RNA and found it to be within the standard range (\u003cstrong\u003eSupplementary Table\u0026nbsp;3\u003c/strong\u003e). Next, we performed RT-qPCR on all samples. By comparing the expression levels of key genes between the PTC group and the control group, we found significant differences in the expression of TSHZ3, SERGEF, AKAP12, SGPP2, and ASGR1 \u003cstrong\u003e(\u003c/strong\u003eFig. \u003cspan class=\"InternalRef\"\u003e10\u003c/span\u003eA-\u003cspan class=\"InternalRef\"\u003e10\u003c/span\u003eE\u003cstrong\u003e)\u003c/strong\u003e, with ASGR1 showing the most significant difference (p\u0026thinsp;\u0026lt;\u0026thinsp;0.0001). However, the expression differences of AK1 and PELI2 between the two groups were not significant \u003cstrong\u003e(\u003c/strong\u003eFig. \u003cspan class=\"InternalRef\"\u003e10\u003c/span\u003eF,\u003cspan class=\"InternalRef\"\u003e10\u003c/span\u003eG\u003cstrong\u003e)\u003c/strong\u003e, but this may have been due to the limited sample size.\u003c/p\u003e\n\u003c/div\u003e"},{"header":"Discussion","content":"\u003cp\u003eIn this study, 1007 DEGs were obtained by GEO database analysis, and 6628 module genes were obtained by ssGSEA and WGCNA after the expansion of 9 parthanatos related genes. 82 intersection genes were obtained by the intersection of DEGs and module genes. Then GO and KEGG enrichment analysis and PPI protein interaction analysis were performed. The biological functions of GO enriched mainly included cellular cation homeostasis, divalent inorganic cation homeostasis, cellular metal ion homeostasis, synapse organization, negative regulation of transport, cell cortex, clathrin-coated vesicle membrane, coated vesicle membrane, clathrin-coated vesicle, platelet dense tubular network, protein kinase A binding, SH3 domain binding, clathrin adaptor activity, cargo adaptor activity and nucleobase-containing compound kinase activity.\u003c/p\u003e\u003cp\u003eWe performed univariate Cox regression and LASSO regression based on these 82 genes, and finally used 7 genes to construct the risk model, namely TSHZ3, SERGEF, AKAP12, SGPP2, ASGR1, AK1, and PELI2. Among the 7 genes, SERGEF, AKAP12 and AK1 had been reported in thyroid cancer, especially AK1 had been covered in PTC. In an analysis of the gene expression and CNA profiles of Chernobyl Pediatric Patients, AK1 was found only in children who had been exposed to radiation, which means radiation altered AK1 gene expression\u003csup\u003e[\u003cspan citationid=\"CR36\" class=\"CitationRef\"\u003e36\u003c/span\u003e]\u003c/sup\u003e. In addition, AK1 is associated with prognosis in renal cell carcinoma\u003csup\u003e[\u003cspan citationid=\"CR37\" class=\"CitationRef\"\u003e37\u003c/span\u003e]\u003c/sup\u003e, myeloid leukemia\u003csup\u003e[\u003cspan citationid=\"CR38\" class=\"CitationRef\"\u003e38\u003c/span\u003e]\u003c/sup\u003e, and lung adenocarcinoma\u003csup\u003e[\u003cspan citationid=\"CR39\" class=\"CitationRef\"\u003e39\u003c/span\u003e]\u003c/sup\u003e. Hugo Gomez-Rueda \u003cem\u003eet al.\u003c/em\u003e used the differences in gene expression values and a multivariate search tool coupled with a nearest centroid classifier to identify SERGEF and other gene as the biomarker which could improve the diagnosis of cytologically indeterminate thyroid cancers\u003csup\u003e[\u003cspan citationid=\"CR40\" class=\"CitationRef\"\u003e40\u003c/span\u003e]\u003c/sup\u003e. Similar to our study, Yin Tian \u003cem\u003eet al\u003c/em\u003e. also used bioinformatics methods to study the potential regulatory mechanism of thyroid cancer prognosis, and ultimately constructed lasso risk factors using 9 genes including AKAP12\u003csup\u003e[\u003cspan citationid=\"CR41\" class=\"CitationRef\"\u003e41\u003c/span\u003e]\u003c/sup\u003e. Although the other prognostic genes have not been reported in thyroid cancer, they had been studied in other cancers. TSHZ3 had been constructed a prognostic model in Bladder Cancer\u003csup\u003e[\u003cspan citationid=\"CR42\" class=\"CitationRef\"\u003e42\u003c/span\u003e]\u003c/sup\u003e, ASGR1 was studied in Hepatocellular carcinoma\u003csup\u003e[\u003cspan citationid=\"CR43\" class=\"CitationRef\"\u003e43\u003c/span\u003e]\u003c/sup\u003e, PELI2 was studied in myeloma\u003csup\u003e[\u003cspan citationid=\"CR44\" class=\"CitationRef\"\u003e44\u003c/span\u003e]\u003c/sup\u003e and gastric cancer\u003csup\u003e[\u003cspan citationid=\"CR45\" class=\"CitationRef\"\u003e45\u003c/span\u003e]\u003c/sup\u003e, and SGPP2 was studied in prostate cancer\u003csup\u003e[\u003cspan citationid=\"CR46\" class=\"CitationRef\"\u003e46\u003c/span\u003e]\u003c/sup\u003e. The most extensively studied prognostic genes are still used to construct prognostic related models, which is similar to our research. This indicates that using bioinformatics to construct risk models is a direction for predicting disease prognosis.\u003c/p\u003e\u003cp\u003eAfterwards, we conducted in-depth explore on prognostic genes and risk scores, including correlation analysis with clinical features. The results of clinical feature manifest that the survival rate of the high-risk cohorts is significant inferior to the low-risk cohorts in aged, gender, stage3, T3, N, and M0. The prognostic genes TSHZ3 and ASGR1 have also been found in the different stages of the studies with other cancers. TSHZ3 and other 5 gene were selected as differentially expressed homeobox genes and were used to construct prognostic model in Bladder cancer\u003csup\u003e[\u003cspan citationid=\"CR42\" class=\"CitationRef\"\u003e42\u003c/span\u003e]\u003c/sup\u003e. In the analysis of clinical status, the expression of TSHZ3 was positively correlated with T stage, N stage, and BLCA stage. Yinjiang Zhang \u003cem\u003eet al\u003c/em\u003e. found that the expression of ASGR1 was significantly correlated with T stage in Hepatocellular Carcinoma through public database, and the higher expression level exhibited a longer disease-free survival rate\u003csup\u003e[\u003cspan citationid=\"CR47\" class=\"CitationRef\"\u003e47\u003c/span\u003e]\u003c/sup\u003e. Although TSHZ3 and ASGR1 have been studied in different types of cancer, these findings suggest that TSHZ3 and ASGR1 have prognostic roles in different clinical stages or TNM stages.\u003c/p\u003e\u003cp\u003eStudies has shown that immune system is associated with the occurrence, invasion, and metastasis of thyroid cancer, and this may affect the treatment and prognosis of patients\u003csup\u003e[\u003cspan citationid=\"CR48\" class=\"CitationRef\"\u003e48\u003c/span\u003e]\u003c/sup\u003e. In our study, analysis of immune infiltration closely followed analysis of clinical features. 28 types of immune infiltrating cells show differences between high and low-risk cohorts, and the expression of 28 immune cells revealed a positive correlation, which demonstrated that this signature was significantly correlated with immune infiltration. In the following analysis of the correlation between differential immune cells and prognostic genes, we found that CD56 and TSHZ3 showed the highest negative correlation and AK1 showed the highest positive correlation. TSHZ3 and AK1 were associated with immune cells were also found in some non-thyroid cancer studies. In a study of Bladder Cance, there were differences in the infiltration of most immune cells between high-risk cohorts and low-risk cohorts, and TSHZ3 was related to almost all immune cell types\u003csup\u003e[\u003cspan citationid=\"CR42\" class=\"CitationRef\"\u003e42\u003c/span\u003e]\u003c/sup\u003e. AK1 was associated with CD8+ T cells, and it might be involved in immune-related mechanisms were coverd in dilated cardiomyopathy with fibrosis\u003csup\u003e[\u003cspan citationid=\"CR49\" class=\"CitationRef\"\u003e49\u003c/span\u003e]\u003c/sup\u003e.\u003c/p\u003e\u003cp\u003eIn addition to the above mentioned contents, we also analyzed immune checkpoint analysis, GSEA functional enrichment analysis, drug sensitivity, functional similarity analysis, the construction of GeneMANIA network and ceRNA network for prognostic genes, and the mRNA and protein expression of prognostic genes between PTC group and normal group.\u003c/p\u003e\u003cp\u003eSince the data of this study were from TCGA and GEO, its reliability cannot be guaranteed. In addition, although we verified the mRNA expression levels of prognostic genes, our sample size was insufficient (n = 10) and a large number of additional samples were needed. In order to further study the molecular mechanism of these 7 prognostic genes, we can not only detect their mRNA expression, but also carry out further studies in the laboratory through cell or animal experiments. For the accuracy of the risk model, we had validated with the data from TCGA, but But the amount of data needs to be replenished.\u003c/p\u003e\u003c/div\u003e\u003cp\u003eIn this study, 1007 DEGs were obtained by GEO database analysis, and 6628 module genes were obtained by ssGSEA and WGCNA after the expansion of 9 parthanatos related genes. 82 intersection genes were obtained by the intersection of DEGs and module genes. Then GO and KEGG enrichment analysis and PPI protein interaction analysis were performed. The biological functions of GO enriched mainly included cellular cation homeostasis, divalent inorganic cation homeostasis, cellular metal ion homeostasis, synapse organization, negative regulation of transport, cell cortex, clathrin-coated vesicle membrane, coated vesicle membrane, clathrin-coated vesicle, platelet dense tubular network, protein kinase A binding, SH3 domain binding, clathrin adaptor activity, cargo adaptor activity and nucleobase-containing compound kinase activity.\u003c/p\u003e\u003cp\u003eWe performed univariate Cox regression and LASSO regression based on these 82 genes, and finally used 7 genes to construct the risk model, namely TSHZ3, SERGEF, AKAP12, SGPP2, ASGR1, AK1, and PELI2. Among the 7 genes, SERGEF, AKAP12 and AK1 had been reported in thyroid cancer, especially AK1 had been covered in PTC. In an analysis of the gene expression and CNA profiles of Chernobyl Pediatric Patients, AK1 was found only in children who had been exposed to radiation, which means radiation altered AK1 gene expression\u003csup\u003e[\u003cspan citationid=\"CR36\" class=\"CitationRef\"\u003e36\u003c/span\u003e]\u003c/sup\u003e. In addition, AK1 is associated with prognosis in renal cell carcinoma\u003csup\u003e[\u003cspan citationid=\"CR37\" class=\"CitationRef\"\u003e37\u003c/span\u003e]\u003c/sup\u003e, myeloid leukemia\u003csup\u003e[\u003cspan citationid=\"CR38\" class=\"CitationRef\"\u003e38\u003c/span\u003e]\u003c/sup\u003e, and lung adenocarcinoma\u003csup\u003e[\u003cspan citationid=\"CR39\" class=\"CitationRef\"\u003e39\u003c/span\u003e]\u003c/sup\u003e. Hugo Gomez-Rueda \u003cem\u003eet al.\u003c/em\u003e used the differences in gene expression values and a multivariate search tool coupled with a nearest centroid classifier to identify SERGEF and other gene as the biomarker which could improve the diagnosis of cytologically indeterminate thyroid cancers\u003csup\u003e[\u003cspan citationid=\"CR40\" class=\"CitationRef\"\u003e40\u003c/span\u003e]\u003c/sup\u003e. Similar to our study, Yin Tian \u003cem\u003eet al\u003c/em\u003e. also used bioinformatics methods to study the potential regulatory mechanism of thyroid cancer prognosis, and ultimately constructed lasso risk factors using 9 genes including AKAP12\u003csup\u003e[\u003cspan citationid=\"CR41\" class=\"CitationRef\"\u003e41\u003c/span\u003e]\u003c/sup\u003e. Although the other prognostic genes have not been reported in thyroid cancer, they had been studied in other cancers. TSHZ3 had been constructed a prognostic model in Bladder Cancer\u003csup\u003e[\u003cspan citationid=\"CR42\" class=\"CitationRef\"\u003e42\u003c/span\u003e]\u003c/sup\u003e, ASGR1 was studied in Hepatocellular carcinoma\u003csup\u003e[\u003cspan citationid=\"CR43\" class=\"CitationRef\"\u003e43\u003c/span\u003e]\u003c/sup\u003e, PELI2 was studied in myeloma\u003csup\u003e[\u003cspan citationid=\"CR44\" class=\"CitationRef\"\u003e44\u003c/span\u003e]\u003c/sup\u003e and gastric cancer\u003csup\u003e[\u003cspan citationid=\"CR45\" class=\"CitationRef\"\u003e45\u003c/span\u003e]\u003c/sup\u003e, and SGPP2 was studied in prostate cancer\u003csup\u003e[\u003cspan citationid=\"CR46\" class=\"CitationRef\"\u003e46\u003c/span\u003e]\u003c/sup\u003e. The most extensively studied prognostic genes are still used to construct prognostic related models, which is similar to our research. This indicates that using bioinformatics to construct risk models is a direction for predicting disease prognosis.\u003c/p\u003e"},{"header":"Declarations","content":"\u003cp\u003e\u003cstrong\u003eFunding\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eNot applicable.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eConflict of Interest\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eAll authors declare that they have none of the conflicts of interest.\u0026nbsp;\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eEthics\u0026nbsp;\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eThe protocol was approved by the ethics committee of Third Hospital of Shanxi Medical University(Approval notice number: YXLL-2024-173) in accordance with the Declaration of Helsinki. Informed consents (Consent to Participate and Consent to Publish) were obtained from all participants or, if participants are under 18, from a parent and/or legal guardian.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eAuthor contributions\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eRui Wang and Shuxin Wen\u0026nbsp;conceived and designed the project.\u0026nbsp;Rui Wang\u0026nbsp;and\u0026nbsp;Li Zhang\u0026nbsp;conducted the experiments.\u0026nbsp;Rui Wang\u0026nbsp;analyzed and interpreted the data.\u0026nbsp;Rui Wang,\u0026nbsp;Li Zhang\u0026nbsp;and\u0026nbsp;Shuxin Wen\u0026nbsp;wrote the manuscript.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eAcknowledgment\u0026nbsp;\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eNot applicable.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eClinical trial number\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eNot applicable.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eConsent to Publish declaration\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eNot applicable.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eData Availability Statement\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eThe GSE33630 data are available as open data via the Gene Expression Omnibus database：https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE33630. The TCGA-THCA dataset are available as open data via the Cancer Genome Atlas Program database: https://portal.gdc.cancer.gov/projects/TCGA-THCA.\u003c/p\u003e"},{"header":"References","content":"\u003col\u003e\u003cli\u003e\u003cspan\u003eNone. Guidelines for the diagnosis and treatment of thyroid cancer (2022)[J]. Chin J Practical Surg. 2022;42(12):16.\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eIto Y, Miyauchi A, Fujishima M, Masuoka H, Higashiyama T, Kihara M, Onoda N, Miya A. Prognostic significance of patient age in papillary thyroid carcinoma with no high-risk features. Endocr J. 2022;69(9):1131\u0026ndash;6.\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eLa Vecchia C, Malvezzi M, Bosetti C, Garavello W, Bertuccio P, Levi F, Negri E. Thyroid cancer mortality and incidence: a global overview. Int J Cancer. 2015;136(9):2187\u0026ndash;95.\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eSchindler AM, van Melle G, Evequoz B, Scazziga B. Prognostic factors in papillary carcinoma of the thyroid. Cancer. 1991;68(2):324\u0026ndash;30.\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eShaha AR. Implications of prognostic factors and risk groups in the management of differentiated thyroid cancer. Laryngoscope. 2004;114(3):393\u0026ndash;402.\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eWei X, Wang X, Xiong J, Li C, Liao Y, Zhu Y, Mao J. Risk and Prognostic Factors for BRAFV600E Mutations in Papillary Thyroid Carcinoma. Biomed Res Int. 2022;2022:9959649.\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eXing M. Prognostic utility of BRAF mutation in papillary thyroid cancer. Mol Cell Endocrinol. 2010;321(1):86\u0026ndash;93.\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eChen B, Shi Y, Xu Y, Zhang J. The predictive value of coexisting BRAFV600E and TERT promoter mutations on poor outcomes and high tumour aggressiveness in papillary thyroid carcinoma: A systematic review and meta-analysis. Clin Endocrinol (Oxf). 2021;94(5):731\u0026ndash;42.\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eLiu R, Xing M. TERT promoter mutations in thyroid cancer. Endocr Relat Cancer. 2016;23(3):R143\u0026ndash;55.\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eScheffel RS, Maia AL. The long and still uncertain journey of BRAF as a prognostic tool in patients with papillary thyroid cancer. Arch Endocrinol Metab. 2019;63(2):95\u0026ndash;6.\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eNasirden A, Saito T, Fukumura Y, Hara K, Akaike K, Kurisaki-Arakawa A, Asahina M, Yamashita A, Tomomasa R, Hayashi T, Arakawa A, Yao T. In Japanese patients with papillary thyroid carcinoma, TERT promoter mutation is associated with poor prognosis, in contrast to BRAF V600E mutation. Virchows Arch. 2016;469(6):687\u0026ndash;96.\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eHuang P, Chen G, Jin W, Mao K, Wan H, He Y. Molecular Mechanisms of Parthanatos and Its Role in Diverse Diseases. Int J Mol Sci. 2022;23(13):7292.\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eKoehler RC, Dawson VL, Dawson TM. Targeting Parthanatos in Ischemic Stroke. Front Neurol. 2021;12:662034.\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eXu X, Sun B, Zhao C. Poly (ADP-Ribose) polymerase 1 and parthanatos in neurological diseases: From pathogenesis to therapeutic opportunities. Neurobiol Dis. 2023;187:106314.\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eWang Y, An R, Umanah GK, Park H, Nambiar K, Eacker SM, Kim B, Bao L, Harraz MM, Chang C, Chen R, Wang JE, Kam TI, Jeong JS, Xie Z, Neifert S, Qian J, Andrabi SA, Blackshaw S, Zhu H, Song H, Ming GL, Dawson VL, Dawson TM. A nuclease that mediates cell death induced by DNA damage and poly(ADP-ribose) polymerase-1. Science. 2016;354(6308):aad6872.\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eWang X, Ge P. Parthanatos in the pathogenesis of nervous system diseases. Neuroscience. 2020;449:241\u0026ndash;50.\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eJiang X, Deng W, Tao S, Tang Z, Chen Y, Tian M, Wang T, Tao C, Li Y, Fang Y, Pu C, Gao J, Wang X, Qu W, Gai X, Ding Z, Fu Y, Zheng Y, Cao S, Zhou J, Huang M, Liu W, Xu J, Fan J, Shi Y. A RIPK3-independent role of MLKL in suppressing parthanatos promotes immune evasion in hepatocellular carcinoma. Cell Discov. 2023;9(1):7.\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eZou Y, Xie J, Zheng S, Liu W, Tang Y, Tian W, Deng X, Wu L, Zhang Y, Wong CW, Tan D, Liu Q, Xie X. Leveraging diverse cell-death patterns to predict the prognosis and drug sensitivity of triple-negative breast cancer patients after surgery. Int J Surg. 2022;107:106936.\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eLi C, Zhang J, Wu Q, Kumar A, Pan G, Kelvin DJ. Nifuroxazide Activates the Parthanatos to Overcome TMPRSS2:ERG Fusion-Positive Prostate Cancer. Mol Cancer Ther. 2023;22(3):306\u0026ndash;16.\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eLi C, Xue Y, Wu J, Zhang L, Yang T, Ai M, Han J, Zheng X, Wang R, Boldogh I, Ba X. MTH1 inhibition synergizes with ROS-inducing agents to trigger cervical cancer cells undergoing parthanatos. Biochim Biophys Acta Mol Basis Dis. 2024;1870(5):167190.\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eChen DL, Cai JH, Wang CCN. Identification of Key Prognostic Genes of Triple Negative Breast Cancer by LASSO-Based Machine Learning and Bioinformatics Analysis. Genes (Basel). 2022;13(5):902.\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eWu B, Xi S. Bioinformatics analysis of differentially expressed genes and pathways in the development of cervical cancer. BMC Cancer. 2021;21(1):733.\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eZhang Y, Shan L, Li D, Tang Y, Qian W, Dai J, Du M, Sun X, Zhu Y, Wang Q, Zhou L. Identification of key biomarkers associated with immune cells infiltration for myocardial injury in dermatomyositis by integrated bioinformatics analysis. Arthritis Res Ther. 2023;25(1):69.\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eQin H, Abulaiti A, Maimaiti A, Abulaiti Z, Fan G, Aili Y, Ji W, Wang Z, Wang Y. Integrated machine learning survival framework develops a prognostic model based on inter-crosstalk definition of mitochondrial function and cell death patterns in a large multicenter cohort for lower-grade glioma. J Transl Med. 2023;21(1):588.\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eRitchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W, Smyth GK. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. 2015;43(7):e47.\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eGu Z, Eils R, Schlesner M. Complex heatmaps reveal patterns and correlations in multidimensional genomic data. Bioinformatics. 2016;32(18):2847\u0026ndash;9.\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eH\u0026auml;nzelmann S, Castelo R, Guinney J. GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinformatics. 2013;14:7.\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eLangfelder P, Horvath S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics. 2008;9:559.\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eWu T, Hu E, Xu S, Chen M, Guo P, Dai Z, Feng T, Zhou L, Tang W, Zhan L, Fu X, Liu S, Bo X, Yu G. clusterProfiler 4.0: A universal enrichment tool for interpreting omics data. Innov (Camb). 2021;2(3):100141.\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eShannon P, Markiel A, Ozier O, Baliga NS, Wang JT, Ramage D, Amin N, Schwikowski B, Ideker T. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 2003;13(11):2498\u0026ndash;504.\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eFriedman J, Hastie T, Tibshirani R. Regularization Paths for Generalized Linear Models via Coordinate Descent. J Stat Softw. 2010;33(1):1\u0026ndash;22.\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eLove MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15(12):550.\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eFu D, Zhang B, Wu S, Zhang Y, Xie J, Ning W, Jiang H. Prognosis and Characterization of Immune Microenvironment in Acute Myeloid Leukemia Through Identification of an Autophagy-Related Signature. Front Immunol. 2021;12:695865.\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eMas MT, Chen CY, Hitzeman RA, Riggs AD. Active human-yeast chimeric phosphoglycerate kinases engineered by domain interchange. Science. 1986;233(4765):788\u0026ndash;90.\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eYu G. Gene Ontology Semantic Similarity Analysis Using GOSemSim. Methods Mol Biol. 2020;2117:207\u0026ndash;15.\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eStein L, Rothschild J, Luce J, Cowell JK, Thomas G, Bogdanova TI, Tronko MD, Hawthorn L. Copy number and gene expression alterations in radiation-induced papillary thyroid carcinoma from chernobyl pediatric patients. Thyroid. 2010;20(5):475\u0026ndash;87.\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eLiu J, Liu B, Guo Y, Chen Z, Sun W, Gao W, Wu H, Wang Y. Key miRNAs and target genes played roles in the development of clear cell renal cell carcinoma. Cancer Biomark. 2018;23(2):279\u0026ndash;90.\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eQin T, Zhao H, Shao Y, Hu N, Shi J, Fu L, Zhang Y. High expression of AK1 predicts inferior prognosis in acute myeloid leukemia patients undergoing chemotherapy. Biosci Rep. 2020;40(6):BSR20200097.\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eJan YH, Lai TC, Yang CJ, Huang MS, Hsiao M. A co-expressed gene status of adenylate kinase 1/4 reveals prognostic gene signature associated with prognosis and sensitivity to EGFR targeted therapy in lung adenocarcinoma. Sci Rep. 2019;9(1):12329.\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eGomez-Rueda H, Palacios-Corona R, Guti\u0026eacute;rrez-Hermosillo H, Trevino V. A robust biomarker of differential correlations improves the diagnosis of cytologically indeterminate thyroid cancers. Int J Mol Med. 2016;37(5):1355\u0026ndash;62.\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eTian Y, Xie T, Sun X. Analysis of the regulatory mechanisms of prognostic immune factors in thyroid cancer. Front Oncol. 2022;12:1059591.\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eDong B, Liang J, Li D, Song W, Song J, Zhu M, Zhao S, Ma Y, Yang T. Identification of a Prognostic Signature Associated With the Homeobox Gene Family for Bladder Cancer. Front Mol Biosci. 2021;8:688298.\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eRoa-Colomo A, L\u0026oacute;pez Garrido M\u0026Aacute;, Molina-Vallejo P, Rojas A, Sanchez MG, Aranda-Garc\u0026iacute;a V, Salmeron J, Romero-Gomez M, Muntane J, Padillo J, Alamo JM, Lorente JA, Serrano MJ, Garrido-Navas MC. Hepatocellular carcinoma risk-stratification based on ASGR1 in circulating epithelial cells for cancer interception. Front Mol Biosci. 2022;9:1074277.\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eMasuda T, Haji S, Nakashima Y, Tsuda M, Kimura D, Takamatsu A, Iwahashi N, Umakoshi H, Shiratsuchi M, Kikutake C, Suyama M, Ohkawa Y, Ogawa Y. Identification of a drug-response gene in multiple myeloma through longitudinal single-cell transcriptome sequencing. iScience. 2022;25(8):104781.\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eZhang C, Jing LW, Li ZT, Chang ZW, Liu H, Zhang QM, Zhang QY. Identification of a prognostic 28-gene expression signature for gastric cancer with lymphatic metastasis. Biosci Rep. 2019;39(5):BSR20182179.\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eZhang Y, Kong X, Xin S, Bi L, Sun X. Discovery of Lipid Metabolism-Related Genes for Predicting Tumor Immune Microenvironment Status and Prognosis in Prostate Cancer. J Oncol. 2022;2022:8227806.\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eZhang Y, Wei H, Fan L, Fang M, He X, Lu B, Pang Z. CLEC4s as Potential Therapeutic Targets in Hepatocellular Carcinoma Microenvironment. Front Cell Dev Biol. 2021;9:681372.\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eYin H, Tang Y, Guo Y, Wen S. Immune Microenvironment of Thyroid Cancer. J Cancer. 2020;11(16):4884\u0026ndash;96.\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eLi W, Liu P, Liu H, Zhang F, Fu Y. Integrative analysis of genes reveals endoplasmic reticulum stress-related immune responses involved in dilated cardiomyopathy with fibrosis. Apoptosis. 2023;28(9\u0026ndash;10):1406\u0026ndash;21.\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eSupplementary.\u003c/span\u003e\u003c/li\u003e\u003c/ol\u003e"}],"fulltextSource":"","fullText":"","funders":[],"hasAdminPriorityOnWorkflow":false,"hasManuscriptDocX":true,"hasOptedInToPreprint":true,"hasPassedJournalQc":"","hasAnyPriority":false,"hideJournal":false,"highlight":"","institution":"","isAcceptedByJournal":true,"isAuthorSuppliedPdf":false,"isDeskRejected":"","isHiddenFromSearch":false,"isInQc":false,"isInWorkflow":false,"isPdf":false,"isPdfUpToDate":true,"isWithdrawnOrRetracted":false,"journal":{"display":true,"email":"[email protected]","identity":"discover-oncology","isNatureJournal":false,"hasQc":true,"allowDirectSubmit":false,"externalIdentity":"dion","sideBox":"Learn more about [Discover Oncology](https://www.springer.com/12672)","snPcode":"","submissionUrl":"","title":"Discover Oncology","twitterHandle":"","acdcEnabled":true,"dfaEnabled":true,"editorialSystem":"stoa","reportingPortfolio":"Discover Series","inReviewEnabled":true,"inReviewRevisionsEnabled":true},"keywords":"Papillary thyroid carcinoma, Parthanatos, Prognostic gene, Risk model, Immune checkpoint inhibitors","lastPublishedDoi":"10.21203/rs.3.rs-6673087/v1","lastPublishedDoiUrl":"https://doi.org/10.21203/rs.3.rs-6673087/v1","license":{"name":"CC BY 4.0","url":"https://creativecommons.org/licenses/by/4.0/"},"manuscriptAbstract":"\u003ch2\u003eObject:\u003c/h2\u003e\u003cp\u003eThis study aimed to elucidate the role of parthanatos-related genes (PRGs) in papillary thyroid carcinoma (PTC) and construct a prognostic risk model to guide personalized treatment.\u003c/p\u003e\u003ch2\u003eMethods\u003c/h2\u003e\u003cp\u003eUsing the GSE33630 dataset, differentially expressed PRGs were identified and analyzed via weighted gene co-expression network analysis (WGCNA) to pinpoint key module genes. Regression analysis selected seven prognostic genes for risk model construction. The model\u0026rsquo;s performance was validated, and a nomogram was developed for survival prediction. Further analyses included clinical feature correlations, immune infiltration, drug sensitivity, gene set enrichment analysis (GSEA), and experimental validation via RT-qPCR.\u003c/p\u003e\u003ch2\u003eResults\u003c/h2\u003e\u003cp\u003eSeven prognostic genes (TSHZ3, SERGEF, AKAP12, SGPP2, ASGR1, AK1, PELI2) were identified. The risk model demonstrated robust predictive accuracy, stratifying patients into high- and low-risk groups with significant survival differences. GSEA revealed 29 enriched pathways (e.g., ribosome, focal adhesion), while immune infiltration analysis highlighted CD56\u0026thinsp;+\u0026thinsp;NK cells and AK1 as key immune correlates. Drug sensitivity screening identified 111 differential therapeutics. Functional analysis indicated AKAP12 had the strongest functional similarity among prognostic genes.\u003c/p\u003e\u003ch2\u003eConclusion\u003c/h2\u003e\u003cp\u003e This study comprehensively mapped PRGs in PTC, established a validated risk model, and provided insights into immune-microenvironment interactions and therapeutic targets, advancing precision oncology for PTC.\u003c/p\u003e","manuscriptTitle":"Construction and validation of risk models of prognostic genes associated with parthanatos in papillary thyroid carcinoma based on bioinformatics","msid":"","msnumber":"","nonDraftVersions":[{"code":1,"date":"2025-07-22 16:34:58","doi":"10.21203/rs.3.rs-6673087/v1","editorialEvents":[{"type":"communityComments","content":0},{"type":"decision","content":"Revision requested","date":"2025-09-22T07:56:05+00:00","index":"","fulltext":""},{"type":"editorInvitedReview","content":"","date":"2025-09-07T03:08:44+00:00","index":"hide","fulltext":""},{"type":"editorInvitedReview","content":"","date":"2025-08-12T09:28:27+00:00","index":"hide","fulltext":""},{"type":"reviewerAgreed","content":"262575402701597698731259324887644286300","date":"2025-08-02T01:44:03+00:00","index":"hide","fulltext":""},{"type":"reviewerAgreed","content":"114400738291318079486947575973694348544","date":"2025-07-27T17:33:38+00:00","index":"hide","fulltext":""},{"type":"reviewersInvited","content":"","date":"2025-07-15T14:58:50+00:00","index":"","fulltext":""},{"type":"editorAssigned","content":"","date":"2025-07-15T14:51:35+00:00","index":"","fulltext":""},{"type":"editorInvited","content":"","date":"2025-06-30T13:58:25+00:00","index":"","fulltext":""},{"type":"checksComplete","content":"","date":"2025-06-29T04:01:40+00:00","index":"","fulltext":""},{"type":"submitted","content":"Discover Oncology","date":"2025-06-29T03:57:12+00:00","index":"","fulltext":""}],"status":"published","journal":{"display":true,"email":"[email protected]","identity":"discover-oncology","isNatureJournal":false,"hasQc":true,"allowDirectSubmit":false,"externalIdentity":"dion","sideBox":"Learn more about [Discover Oncology](https://www.springer.com/12672)","snPcode":"","submissionUrl":"","title":"Discover Oncology","twitterHandle":"","acdcEnabled":true,"dfaEnabled":true,"editorialSystem":"stoa","reportingPortfolio":"Discover Series","inReviewEnabled":true,"inReviewRevisionsEnabled":true}}],"origin":"","ownerIdentity":"1ac6eb7e-bd17-4706-9c21-7d06af20430c","owner":[],"postedDate":"July 22nd, 2025","published":true,"recentEditorialEvents":[],"rejectedJournal":[],"revision":"","amendment":"","status":"under-review","subjectAreas":[],"tags":[],"updatedAt":"2025-11-26T13:38:43+00:00","versionOfRecord":[],"versionCreatedAt":"2025-07-22 16:34:58","video":"","vorDoi":"","vorDoiUrl":"","workflowStages":[]},"version":"v1","identity":"rs-6673087","journalConfig":"researchsquare"},"__N_SSP":true},"page":"/article/[identity]/[[...version]]","query":{"redirect":"/article/rs-6673087","identity":"rs-6673087","version":["v1"]},"buildId":"8U1c8b4HqxoKbykW_rLl7","isFallback":false,"isExperimentalCompile":false,"dynamicIds":[84888],"gssp":true,"scriptLoader":[]}