An individualized protein-based prognostic model to stratify pediatric patients with papillary thyroid carcinoma

Abstract

Pediatric papillary thyroid carcinomas (PPTCs) exhibit high inter-tumor heterogeneity and currently lack widely adopted recurrence risk stratification criteria. Hence, we propose a machine learning-based objective method to individually predict their recurrence risk. We retrospectively collect and evaluate the clinical factors and proteomes of 83 pediatric benign (PB), 85 pediatric malignant (PM) and 66 adult malignant (AM) nodules, and quantify 10,426 proteins by mass spectrometry. We find 243 and 121 significantly dysregulated proteins from PM vs. PB and PM vs. AM, respectively. Function and pathway analyses show the enhanced activation of the inflammatory and immune system in PM patients compared with the others. Nineteen proteins are selected to predict recurrence using a machine learning model with an accuracy of 88.24%. Our study generates a protein-based personalized prognostic prediction model that can stratify PPTC patients into high- or low-recurrence risk groups, providing a reference for clinical decision-making and individualized treatment.

Fig. 1. Study overview.

a Study design of analyzed cohort and experiment workflow. Created through Biorender.com. b Enrollment and exclusion criteria for pediatric papillary thyroid carcinoma (PTC), pediatric benign nodule and adult PTC patients.

Fig. 2. Analysis of the clinical recurrence risk factors for pediatric papillary thyroid carcinoma (PPTC).

a Kaplan–Meier survival curves of two groups (red: patients below the median; blue: otherwise) for three significant factors: age, total lymph node metastasis (TLNN) and lateral lymph node metastasis number (LLNN). P values are derived from Log-Rank Test. b, c Forest plots for two multivariate CoxPH models using (b) continuous non-negative integer age and (c) categorical age, respectively. Data are presented as hazard ratio value with 95% confidence interval. P values are tested by Cox proportional hazard model.

Fig. 3. Functional analyses of the dysregulated proteins.

a, b Differentially expressed proteins (DEPs) are shown in the volcano plots: a pediatric malignant nodules (PM) vs. pediatric benign nodules (PB) and b PM vs. adult malignant nodules (AM). The cutoff is defined by requiring the fold change (FC) to be greater than 1.5, with adjusted P < 0.05 (BH-adjusted two-sided Welch’s t test). The names of the up/downregulated proteins with the top five largest FC are reported in the plots. c The scatter diagram shows the FC distribution of the dysregulated proteins in two pairwise comparisons: PM vs. PB and PM vs. AM. The overlapping significantly co-dysregulated proteins are colored in red. The proteins significantly dysregulated in PM/PB are colored in orange, and those dysregulated in PM/AM are colored in blue. Here, the DEP lists were derived from FC threshold 1.5. d The heatmap shows 37 proteins: they are co-upregulated/co-downregulated proteins and the five most up-/down- regulated ones from the volcano plots (a, b). Proteins were clustered using hierarchical clustering. e Pathway enrichment of the 243 DEPs from the volcano plot with the PM/PB comparisons (FC > 1.5). The red and blue bars represent the active and inhibited pathways, respectively. f Results of the gene ontology enrichment of the biological processes using the DEPs in PM/PB. g Pathway enrichment of 121 DEPs from the volcano plot with the PM/AM comparisons (FC > 1.5). The red and blue bars represent the active and inhibited pathways, respectively. P values are derived from one-sided Fisher’s Exact Test for pathway and gene ontology enrichment.

Fig. 4. In silico immune infiltration analysis and expression levels of immune checkpoints.

a Relative proportions of seven types of immune cells in pediatric benign (PB, N = 83) and pediatric malignant (PM, N = 85) samples imputed by CIBERSORTx. Boxes are first and third quartiles, the center line is median, whiskers are ±1.5 interquartile range, and dots are individual data points. Abundance outliers and missing values are not included in the boxplot, throughout Fig. 4. b Representative multiplex immunohistochemistry staining in PB (N = 10) and PM samples (N = 10). The scale bar represents 50 μm. c The protein expression abundances of poliovirus receptor (PVR) and interleukin 10 receptor B (IL10RB) in PB (N = 83), PM-NR (non-recurrence, N = 73) and PM-R (recurrence, N = 12) groups. The significance throughout Fig. 4 is determined by two-sided Welch’s t test.

Fig. 5. Pediatric papillary thyroid carcinoma (PPTC) prognostic prediction.

a The C-indexes of our five models were calculated on training, threefold cross-validation, and test sets. b Density curves of the training  continuous risk ranking (Crank) scores of two groups (recurrence (1) or no recurrence (0)). The Fisher decision boundary was used to differentiate the low- from the high-risk groups. c The Kaplan–Meier survival curves of the low- and high-risk groups, calculated on the training and test sets, show significant differences. d Permutation importance of the 19 proteins from the ProtRsf model. e Network showing the 19 features of the ProtRsf model with the connected proteins enriched using the Ingenuity Pathway Analysis software. f The relative protein abundances of galectin-3 (LGALS3), the hub protein of the network (e), in the four groups. Boxes are first and third quartiles, the center line is median, whiskers are ±1.5 interquartile range, and dots are individual data points. Abundance outliers and missing values are not shown in the boxplot. Biologically independent samples shown in boxplot: PB, N = 77; Low-risk PM, N = 67; High-risk PM, N = 18; AM, N = 60. The mild outliers were removed, and a two-sided unpaired Wilcoxon rank-sum test was used, without continuity correction, to calculate the P values.

Fig. 6. Risk stratification.

a Predicted risk stratification for pediatric papillary thyroid carcinoma (PPTC) patients. The sample indexes of false positives (N = 8) and the false negatives (N = 2) are labeled. b, c The predicted survival curves of the two false negatives (b) and eight false positives (c) with their continuous risk ranking (Crank) scores, sample indexes and recurrence or latest follow-up times (shown by the vertical lines).