1. Introduction
The worldwide incidence of thyroid cancer has been rising rapidly in the past three decades, with the largest contribution being in papillary thyroid cancer (PTC) in all countries analyzed in a global population-based assessment study [1]. While it remains to be investigated whether the over-diagnosis is attributed to increased screening programs, South Korea has the world’s highest rate of thyroid cancer, which increased 15-fold between 1993 and 2011 [2]. Nevertheless, a 10-year disease-specific survival (DSS) by the American Joint Committee on Cancer/Union for International Cancer Control (AJCC/UICC) defined stages I to IV ranging from 81% (stage IV) to 100% (stage II) [3]. Despite being considered as an indolent tumor with low occurrence of local invasion, recurrences and regional or distant metastases, PTC, similar to other types of cancer, show inter- and intra-tumoral heterogeneity with varying degree of genetic diversity, which could have a significant impact on prognosis and on the response to targeted therapy [4,5]. The existence of a small population of tumors having the more aggressive variants of PTC, with distinct clinical, pathological, and molecular features suggest the need for a robust predictive marker for patient stratification disease management [4,5,6].
In thyroid cancer, main driver genetic alterations include BRAF (V600E) and RAS (NRAS involving codons 12 and 61) mutations, RET gene fusions, and PAX8-PPARG gene fusions [7,8]. A PAX8-PPARγ fusion protein (PPFP) is produced when the promoter and most of the PAX8 gene, which encodes an important transcription factor for normal thyroid gland development, fuse with the coding exons of PPARγ, a member of the steroid/thyroid nuclear receptor family [8,9]. PPFP is found in about one-third of follicular thyroid carcinomas (FTC) and can act as an oncoprotein, as evidenced by in vivo and in vitro studies [8,10,11].
Similar to other subtypes of PPAR, PPARγ forms heterodimers with retinoid X receptor alpha (RXRα) and binds to specific DNA sequences, termed as peroxisome proliferator response elements (PPREs), in a ligand-responsive manner to trans-activate target genes [12]. Previous studies have revealed cell-specific maps and 3D structure of the intact PPARγ-RXRα complex, including binding sites within gene promoters and intergenic or intronic regions [12,13,14]. Variants in the PPARG gene were found to reduce the receptor-binding affinity to the PPREs and modulate transcriptional activity of its target genes, such as acyl coenzyme A (acyl-CoA) oxidase, regulating insulin sensitivity and adipose tissue differentiation [15].
While PPARγ is expressed at extremely low levels in the normal thyroid [11], its expression is highly induced in a variety of immune cells, including monocytes and macrophages, in addition to the well-studied adipocytes, governing cellular phenotype and function such as lipid metabolism and secretome through transactivation of PPARγ target genes [11,16,17]. Several key molecular mechanisms regulated by PPARγ in macrophages include differentiation [18], M1-to-M2 polarization [19], lipid metabolism [20], suppression of the production of pro-inflammatory cytokines (e.g., TNFα, IL-1B, and IL-6), and the expression of inflammation-related genes (e.g., iNOS and MMP9) [8,21].
Here, we hypothesized that the PPARγ target genes could serve as a prognosticator of outcome in patients with PTC due to the common variants found in an intronic region located between exon 5 and exon 6 of the PPARG gene encoding the ligand-binding domain (LBD) of nuclear hormone—a key domain for transactivation of PPARγ targets [11,22]—in three pairs of monozygotic twins harboring PTC. A newly developed personalized scoring index computed based on the expression levels of PPARγ target genes applied to merged microarray data (MMD) consisting exclusively of thyroid cancers and normal tissues and TCGA-THCA (thyroid carcinoma) data revealed its robust diagnostic and prognostic ability in predicting tissue-based disease-specific survival (DSS). Leveraging single-cell RNA sequencing (scRNA-seq) data and machine learning algorithms, we further identified specific cell types contributing predominantly to the PPARGi and demonstrated clinical significance and applicability of PPARGi assay in a routine clinical setting, regardless of profiling platform.
2. Materials and Methods
2.1. Subjects
We performed a retrospective study targeting identical twins who were treated with surgery for papillary thyroid cancer at Gangnam Severance Hospital, Yonsei University College of Medicine, between May 2009 and March 2020. Three pairs of identical twins with papillary thyroid carcinoma, denotated as TP1, TP2, and TP3 were identified. All methods were carried out in accordance with relevant guidelines and regulations. This study protocol was approved by the Institutional Review Board of Yonsei University (IRB 3-2020-0281). The IRB of Yonsei University waived the requirement for patient-informed consent as the study is retrospective by design. The clinicopathologic characteristics of the patients and tumors including age, sex, operation date, surgical extent, tumor size, and central and lateral lymph node metastasis are summarized in Table S1.
2.2. Targeted DNA Sequencing and Analysis
Genomic DNA was extracted from FFPE-fixed tumor tissues. Sequencing libraries were prepared by Macrogen (Seoul, Korea) using SureSelect Target Enrichment kit (Agilent Technologies). Two distinct target panels were designed to detect fusion genes and mutations in coding exons (Table S2) using SureSelect Custom DNA Target Enrichment Probes. The libraries were subjected to Illumina platform in paired-end (2 × 150 bp) mode. Analytical platforms used by Macrogen include (1) FASTQC, fastp (quality check and trimming); (2) BWA, PICARD, SAMTOOLS, and BEDTOOLS (alignment); and (3) MuTect2 (GATK) and LUMPY (variant calling). The human assembly GRCh37/hg19 was used for reference genome. The variant call format (VCF v4.2) provided for each sample by Macrogen was used to identify variants (annotated with SnpEff v4.3), in which only the passing variants annotated as PASS were considered as true variants.
2.3. Mapping of Protein Mutation Data
MutationMapper tool [23] provided by cBioPortal (
2.4. MSKCC-ATC and TCGA-THCA Data
Nonsynonymous mutation data and clinical data were obtained directly from the cBioPortal (
2.5. MMD-THCA Data
Merged microarray-acquired Data (MMD) were generated for THCA comprising non-tumor (NT), ATC, PDTC, and PTC tissues, as previously carried out for other major cancer types [30,31,32,33]. The minimum information about a microarray experiment (MIAME) compliant datasets were carefully selected using gene expression omnibus (GEO), a public functional genomics data repository (
2.6. ScRNA-Seq Data
Expression matrices with molecule counts per gene per cell index of PTC tissues from 14 patients were obtained directly from GEO under the accession code GSE158291 [36]. For QC and preprocessing of scRNA-seq data, cells having unique feature counts over 10,000 or less than 200 counts and genes having less than 10 molecules across the cells were filtered out. The filtered data were used to perform normalization, feature selection (i.e., identification of highly variable genes), linear transformation (i.e., scaling the data), dimensional reduction (i.e., principal component analysis), cell clustering, and non-linear dimensional reduction (i.e., UMAP) using the R Seurat package [37] (v4.0.1). Data from nodular goiter were excluded from further analyses. Cell type identity was assigned manually to each cluster based on differentially expressed features using the FindAllMarkers function (Table S3). Expression levels of the identified cell cluster-specific markers were assessed in a scRNA-seq study named “ICA: Ileum Lamina Propria Immunocyte (Sinai)” through the Single Cell Portal (
2.7. DE, GSEA, and GO Analysis
Differential expression (DE) analyses were performed using the linear modeling features of the R limma package [27] (v3.46.0). Statistical cutoffs of |log2FC| > 1 and adjusted p-value < 0.05 were applied to determine genes to be differentially expressed in both TCGA-THCA and MMD-THCA data. Gene set enrichment analysis (GSEA) was performed using the R fgsea package [38] (v1.16.0) to assess the enrichment of “GOBP_THYROID_HORMONE_GENERATION” gene set, which was downloaded from the Molecular Signatures Database (MSigDB), in PPARGihigh tumors. Gene ontology (GO) analysis and the enrichment analysis of disease-gene associations were performed using the enrichGO function and the enrichDGN function in the R clusterProfiler package [39] (v3.18.1), respectively.
2.8. CIBERSORT
The proportion of immune cell types in TCGA-THCA and MMD-THCA data was estimated using CIBERSORT (
2.9. ROC and Survival Analysis
The receiver operator characteristic (ROC) area under the curve (AUC) and the most optimal threshold (i.e., the threshold with the highest sum of sensitivity and specificity) were computed to evaluate the diagnostic accuracy of the PPARGi in deciding the presence of the disease using the R pROC package [40] (v1.17.0.1). Univariate Cox proportional hazards (PH) regression and Kaplan–Meier (KM) survival analyses were performed using the R survival package (v3.2.11). Samples with missing either disease-specific survival (DSS) data or gene expression data were excluded from the analysis. The regression coefficients, the Wald statistic p-values, hazard ratios (HR), confidence intervals (CIs) of the HR, and log-rank statistics of PPARγ target genes in predicting DSS in TCGA-THCA are listed in Table S4. As previously described [30,32,41], the optimal cut-off used in the survival analysis was determined using Cutoff Finder [42] (
2.10. PPARGi Derivation
The PPARgene database (
2.11. Machine Learning Algorithms for Disease Selection
Orange (v3.29.3) was used for t-SNE visualization, evaluation of ML models, and generation of confusion matrices and ROC curves. The schematic workflow is shown in Figure S1. The processed TCGA-THCA, MMD-THCA (PTC), and MMD-THCA (ATC) data were sent to the ‘test and score’ widget, in which multiple ML models were tested with the following defined parameters: (1) k-nearest neighbors (kNN, number of neighbors = 5, metric = Euclidean, weight = uniform); (2) support vector machine (SVM, cost = 1.00, kernel = RBF, numerical tolerance = 0.001, iteration limit = 100); (3) neural network (neurons in hidden layers = 100, activation = ReLu, solver = Adam with regularization = 0.0001, maximum number of iterations = 200); (4) logistic regression (regularization type = lasso (L1), strength = C1); and (5) random forest (number of trees = 10, subsets smaller than 5 were not split). These five models were evaluated using 10-fold cross-validation. The evaluation results are summarized in Tables 1–3 using the following defined parameters: (1) AUC = area under the ROC; (2) CA = classification accuracy, defined as the proportion of correctly classified examples; (3) F1 = weighted harmonic mean of precision and recall; (4) precision = proportion of TP among instances classified as positive; and (5) recall = proportion of TP among all positive instances in the given data. For cross-platform analyses, the ML models trained on TCGA-THCA data, which were TDM-transformed using the R TDM package (v0.3), as recently shown [33], were applied to MMD-THCA data consisting exclusively of PTC.
3. Results
3.1. Targeted NGS of Thyroid Cancer-Related Genes in Monozygotic Twins with PTC
Formalin-fixed paraffin-embedded (FFPE) tissues derived from three pairs of identical twins with PTC, denoted as TP1, TP2, and TP3, were subjected to targeted next generation sequencing (NGS) for 103 thyroid cancer-related genes (see Section 2.2 of Materials and Methods; Table S2). A total number of variants observed in twins varied from 98 (TP1) to 155 (TP2), in which the pairwise concordance varied from 34.5% (TP3) to 69.0% (TP2) (Figure S2A). The SnpEff-annotated effects of the identified variants were mostly intronic variants, while other variants including 3′ and 5′ UTR variants, disruptive in-frame deletion, missense variants, synonymous variants, upstream and downstream gene variants, 5′ UTR premature start codon gain variants, and splice region variants were detected at lower frequency (Figure S2B).
Notably, all three pairs of twins with PTC had four common variants, of which three (75%) occurred in the PPARG gene (Figure 1A) and the other variant occurred in the TERT gene. Using the MutationMapper at cBioPortal (see Section 2.3 of Materials and Methods), we found that two of the variants of PPARG (chr3: 12468710 and chr3: 12470239) occurred within the intronic region located between exon 5 and exon 6 of the PPARG gene encoding the “Hormone_recep” domain (PF00104, domain source: Pfam), which is defined as ligand-binding domain (LBD) of nuclear hormone receptor (Figure S2C). These genomic positions mapped onto the previously reported 3D structure of the intact PPAR gamma–retinoid X receptor (RXR) alpha (PPARγ-RXRα) nuclear receptor complex (PDB identifier: 3E00) are shown in Figure S2D. Thyroid carcinomas harboring the CDS mutations in the LBD of PPARG were further found in COSMIC GRCh37 v94 database (Figure S3) [24].
Of the SnpEff-annotated variants, those predicted with moderate impact in at least one of the twin pairs were found in ARID1A, ALK, MSH2, FN1, KMT2A, TUBA3C, BRAF, STRN, KMT2D, and NF1 (Figure S4A). The first twin pair (TP1) had the greatest number of missense variants in these genes, including an in-frame deletion in MSH2. They were discordant for missense variants in TUBA3C (p. Tyr262Cys) and BRAF (p. Val640Glu). For TP2, concordant variants were found in BRAF (p. Val640Glu) and STRN (p. Val620Leu), while missense variants in KMT2D (p. Gln827His) were present in only one of the twins. TP3 harbored missense variants in FN1 (p. Asn172Asp) and NF1 (p. Glu836Ala), which was present only in one of the twins who had recurrent cancer. To explore the frequency of nonsynonymous mutations in these identified genes in thyroid cancers, we leveraged cBioPortal-derived MSKCC and TCGA data (see Section 2.4 of Materials and Methods) for thyroid carcinoma (TCGA-THCA, n = 500) and poorly differentiated and anaplastic thyroid cancers (MSKCC-ATC, n = 117). Except for BRAF mutation occurring in 60% and 37% of PTC and ATC, respectively, nonsynonymous variants of ARID1A, ALK, MSH2, FN1, KMT2A, TUBA3C, STRN, KMT2D, and NF1 were found in about 0.2–1.4% and 2.6–6% of patients in TCGA-THCA and MSKCC-ATC, respectively (Figure S4B,C). These variants include putative driver mutations (in-frame variant, missense variant, and truncating variant) as well as amplification and structural variants (Figure S5). The PathwayMapper tool provided by cBioPortal [29] further identified RTK-RAS pathway as the most altered signaling pathway with alterations in ALK, BRAF, and NF1 (score = 3.00; Figure S6).
3.2. Clinical Significance of Pparγ Target Genes in Thyroid Cancer
To examine the diagnostic and prognostic significance of target genes of PPARγ, in which all pairs of twins with PTC harbored the intronic variants, we next shortlisted the experimentally verified PPARγ target genes using PPARgene database (see Section 2.10 of Materials and Methods). Previous works have reported upregulation of the identified targets by PPARγ in human-derived tissues or cells (Table S5). Univariate Cox regression survival analyses of TCGA-THCA-derived RNA-seq revealed that of the 39 targets, 10 PPARγ targets (25.6%) were favorable prognostic factor for disease-specific survival (DSS; Table S4). To construct a personalized scoring system for patient classification, we developed an index termed as PPARGi, which is computed based on the regression coefficient and expression level of the 10 PPARγ target genes (Figure 1B,C). PPARGi varied greatly across tumor samples (n = 505) and was significantly higher than that of normal tissues (n = 59; Wilcoxon p-value = 2.2 × 10−16; Figure 1D). Using the area under the receiver operating characteristic (ROC) curve (AUC), we found that PPARGi distinguishes cancers from normal tissues with the AUC of 0.876 (Figure 1E), demonstrating the potential diagnostic use of PPARGi in deciding the presence of the disease.
PPARGihigh (n = 458; 93.5%) and PPARGilow (n = 32; 6.5%) tumors stratified by the optimal cut-off index determined by Cutoff Finder (see Section 2.9 of Materials and Methods) had significantly different DSS outcomes (Figure 1F; log-rank p-value = 5.8 × 10−16). Through differential expression (DE) analysis between the two stratified groups, we found a total of 1775 DE genes with fold change (log2-base) > 1 and adjusted p-value < 0.05, as shown in Figure 1G. Intriguingly, gene set enrichment analysis (GSEA) revealed the enrichment of genes related to thyroid hormone generation (GO:0006590) in PPARGilow tumors (normalized enrichment score = 1.72 and adjusted p-value = 0.008), which were found to highly express DIO1, TPO, IYD, DIO2, TG, DUOXA2, FOXE1, PAX8, SLC5A5, and DUOX2 (Figure 1H). Additionally, gene ontology (GO) analysis of downregulated genes revealed “T cell activation” and ”regulation of T cell activation” as top-enriched terms in PPARGilow tumors (Figure 1I). These findings are further corroborated by the CIBERSORT analysis, which identified regulatory T cells (Treg) as one of the immune cell types, in which their composition is positively corelated with PPARGi (Figure 1J). In addition to Treg, relative abundance of M0 and M2 macrophages, resting dendritic cells, and activated mast cells increased with PPARGi, while that of CD4 naïve T cells, monocytes, activated natural killer (NK) cells, eosinophils, plasma cells, and T follicular helper (Tfh) cells decreased (Pearson’s correlation p-value < 0.05).
3.3. Validation of PPARGi in MMD-THCA
To validate diverse expression levels of PPARγ target genes and their potential diagnostic value in thyroid cancer, we next generated a unified, cancer type-specific, merged microarray dataset (MMD-THCA), as previously carried out for other major cancer types [30,31,32,33]. Briefly, four independent transcriptomic datasets (GSE29265, GSE33630, GSE65144, and GSE76039) derived from affymetrix human genome U133 Plus 2.0 (HG-U133_Plus_2) array were normalized, integrated, and corrected for batch effects (see Section 2 Materials and Methods). These datasets comprised samples derived from anaplastic thyroid cancer (ATC; n = 52), poorly-differentiated thyroid cancer (PDTC; n = 17), papillary thyroid cancer (PTC; n = 69), and normal tissues (n = 78). As depicted in UMAP representation (Figure 2A), the newly developed MMD exhibited an overlay of samples colored by the source of data and a clear separation between samples derived from different tissue type (i.e., ATC, PDTC, PTC, and normal), demonstrating successful removal of batch effects arising from different studies. Technical validity of the ComBat-transformed MMD was further examined through DE analysis of PTC and normal tissues, in which the enrichment analysis of disease-gene associations (see Section 2.7 of Materials and Methods) identified “carcinoma, papillary”, “follicular adenoma”, “thyroid gland follicular adenoma”, and “follicular thyroid carcinoma” as top enriched terms in PTC (Figure 2B).
MMD-THCA exhibited heterogeneous expression of PPARGi-comprising genes, as previously seen in TCGA-THCA (Figure 2C). Consistent with our findings in TCGA-THCA, PPARGi of thyroid cancers was significantly higher than that of normal tissues (Wilcoxon p-value ≤ 0.0001), of which PTC showed the highest mean of PPARGi among thyroid cancers (Figure 2D). Further, PPARGi achieved the AUC of 0.760, 0.654, and 0.899 in distinguishing cancers from normal tissues in ATC, PDTC, and PTC, respectively (Figure 2E). To validate the enrichment of genes related to thyroid hormone generation (GO:0006590) in PPARGilow tumors observed in TCGA-THCA, we next stratified MMD-PTC into two groups, such that the proportion of PPARGilow tumors in MMD (n = 5; 7.25%) would be comparable to TCGA stratification. DE analysis identified a total of 549 DE genes (Figure 2F), of which the query gene ontology term was enriched in PPARGilow tumors (normalized enrichment score = 1.64 and adjusted p-value = 0.012), which were found to highly express DIO1, TPO, IYD, PAX8, DIO2, FOXE1, DUOXA2, TG, and SLC5A5 (Figure 2G). While we did not find the correlation between PPARGi and the CIBERSORT-estimated proportion of Tregs previously observed in TCGA-THCA, resting dendritic cells and M2 macrophages repeatedly showed statistical significance (Pearson’s correlation p-value < 0.05), highlighting their robust association with PPARGi-comprising genes (Figure 2H).
3.4. Single-Cell Analysis of Ppargi-Comprising Genes in PTC
To identify specific cell types expressing 10 PPARGi-comprising genes in PTC, we next processed and analyzed a scRNA-seq dataset (GSE158291) derived from PTC tissues from 14 patients. Seurat (v4.0.1) identified a total of 15 cell clusters in PTC comprising 4045 QC-passed cells, including pericytes, plasma cells, group 1 innate lymphoid cells (ILC), T cells, B cells, endothelial cells, fibroblasts, epithelial cells, and macrophages (Figure 3A). Each cell cluster was defined manually using differentially expressed features (Figure 3B, Table S3), in which the expression levels of immune cell type-specific features were validated using an independent scRNA-seq dataset obtained from the Single Cell Portal (see Section 2.6 of Materials and Methods). Notably, PPARGi-comprising genes were expressed predominantly in epithelial cells and macrophages (Figure 3C,D). While GPD1, CYP27A1, and REN were rarely expressed in PTC, the remaining genes including DBI, APOE, SAT1, CDKN1A, KLF4, and PLIN2 showed highly heterogenous levels of expression contributing to varied range of PPARGi (Figure 3E). GO analysis of PPARGi-comprising genes identified lipid localization and transport, cholesterol/sterol/steroid catabolic process, and phospholipid transport as top enriched terms (Figure 3F), consistent with the known and putative roles of PPARγ in governing lipid storage, glucose, and insulin process [12].
3.5. Machine Learning for Disease Selection and Risk Stratification
To evaluate the diagnostic performance of expression profiles of PPARGi-comprising genes, we applied different machine learning (ML) algorithms to TCGA-THCA and MMD-THCA data annotated with the origin of tissue (see Section 2.11 of Materials and Methods; Figure S1). As PPARGi was derived from PTC expression data, cancers (TT) separated clearly from normal (NT) tissues in both TCGA-THCA (n = 564) and MMD-THCA (PTC, n = 147) and to a lesser extent in MMD-THCA (ATC, n = 130), as shown in Figure 4A. Of the tested ML models, LASSO multinomial logistic regression (LASSO) achieved the AUC of 0.990 in classifying TT from NT tissues in TCGA-THCA, demonstrating near-perfect diagnostic performance (Table 1 and Figure 4B). Similarly, thyroid tissues were best classified into TT and NT with SVM (AUC = 0.937) and LASSO (AUC = 0.945) in MMD-THCA PTC and ATC, respectively (Table 1 and Figure 4B).
We next aimed to demonstrate the clinical utility of PPARGi genes across different profiling platforms by developing the ML model from RNA-seq-acquired TCGA-PTC data and applying it to the microarray-acquired MMD-PTC data. Notably, a k-nearest neighbors (kNN) model developed from TDM-transformed TCGA-PTC data achieved the best AUC value of 0.942 in classifying TT from NT tissues in MMD-PTC, while the rest of the tested ML models achieved the AUC ranging from 0.890 (SVM) to 0.925 (neural network), showing the clinical applicability of PPARGi genes regardless of profiling platform (Table 2 and Figure 4C,D).
Lastly, we sought to evaluate the predictive potential of ML models in stratifying patients into PPARGihigh and PPARGilow groups, which were found to have significantly different prognosis in terms of DSS (Table 3 and Figure 4D). The tested ML models achieved the AUC ranging from 0.965 (kNN) to 0.998 (SVM) in distinguishing the two risk groups, highlighting the clinical significance of PPARGi genes predictive of survival outcomes in thyroid cancer. These data altogether show the promising results of ML-optimized models built from expression profiles of PPARGi genes for clinical applications, which can easily be translated to other sequencing platforms.
4. Discussion
Despite the well-established role of the PAX8-PPARγ fusion oncogene, the expression landscape of experimentally validated PPARγ target genes remains unclear in thyroid carcinoma. In this study, we found that three pairs of monozygotic twins with PTC harbored four common variants, of which two variants occurred in the LBD of nuclear hormone in the PPARG gene, which have not been reported in healthy Korean population (Table S6) [44]. We further showed that low expression levels of many of the PPARγ target genes were associated with poor survival outcomes, reinforcing the importance of PPARγ-RXRα pathway in governing immune microenvironment in PTC. Through scRNA-seq data analysis, we found that PPARGi-comprising genes were expressed most strongly in the epithelial cells and macrophages among the 15 cell clusters found in patient-derived PTC tissues. Lastly, ML models developed from RNA-seq-derived TCGA-THCA achieved a near-perfect performance in selecting the disease from benign tissues and in stratifying patients into the two risk groups in microarray-derived MMD-THCA, highlighting potential applicability of expression profiles of PPARGi-comprising genes in the clinical setting for patient management.
The LBD, situated in the C-terminus, is connected to the DNA-binding domain (DBD) via a flexible hinge region, which interacts physically with the DNA in PPARγ [11]. It is a key domain for transactivation and transrepression of PPARγ target genes that play important roles in adipogenesis, insulin sensitization, lipid metabolism, and inflammation, making PPARγ an effective target for the management of metabolic diseases, such as type 2 diabetes, obesity, and atherosclerosis [22,45]. The PPARG gene has 15 transcripts (splice variants), of which four transcripts lack exon 5 and 6, which encode LBD (Figure S3A). Skipping of PPARG exon 5, for example, induced by ligand-mediated PPARγ activation, has been proposed as an alternative splicing event regulating PPARγ activity and PPARγ-related diseases [12].
In cancer, the tumorigenic role of PPARγ remains highly controversial. Studies have observed the inhibitory effects of PPARγ-RXRα signaling pathway on tumor growth, angiogenesis, differentiation, and that of production of inflammatory cytokines and tumor invasiveness, suggesting the anti-tumorigenic role of PPARγ in several cancer types including colon, lung, pancreas, prostate, and breast cancer [11,45]. By contrast, the pro-tumorigenic role of PPARγ has been reported in a variety of cancers, such as bladder tumor, renal pelvic tumors, hemangioma, lipoma, skin fibrosarcoma, mammary adenocarcinoma, and hepatic tumors, through distinct molecular mechanisms activated by PPARγ ligands regulating cancer cell proliferation, angiogenesis, and metastasis [46].
The cancer-related loss-of-function PPARγ mutations have been found predominantly throughout the LBD of the PPARG gene with varying degree of impaired ability in inducing transactivation of target genes [11]. Here, we observed that PPARγ target genes were expressed at significantly higher levels compared with benign tissues, suggesting that the PPARγ activation might increase the risk of developing thyroid cancer. Further, among PTC tissues, the low expression levels of PPARGi-comprising genes were closely associated with poor survival outcomes, indicating that the loss-of-function PPARγ variants might deregulate tumor microenvironment, including immune cell-infiltration, through macrophages in PPARGilow tumors. In line with the findings supporting the role of PPARγ-RXRα pathway in modulating immune cell-infiltration, genomic alterations of PPARγ-RXRα complex (i.e., PPARγHigh/RXRαS427F/Y) induced evasion of immunosurveillance and partial resistance to immunotherapies in muscle-invasive bladder cancer [47].
Through analysis of microarray- and RNA-seq-generated, patient-derived, tissue-level gene expression and scRNA-seq data, we assessed the expression profiles of diagnostic and prognostic PPARγ targets comprehensively in a cell type-specific manner and identified specific cell populations contributing the most to the computation of PPARGi in PTC. The last decade has seen the emergence of PPARγ as a key regulator of inflammatory and immune responses particularly in monocytes and macrophages, with anti-inflammatory effects in several disease models and clinical studies [48,49,50]. These studies have suggested that, in macrophages, PPARγ represses pro-inflammatory genes such as TNFα, IL-1B, IL-6, IL-12 MCP-1, and MMP-9 through unique ligand-dependent transcriptional mechanisms [11,21,48,51,52,53]. We thus carefully speculate that the poor prognosis observed in PPARGilow PTC tumors might be attributed in part to the suppressed repression of the production of these pro-inflammatory genes in which the expression levels of many of the genes are potential predictors for metastasis or shorted survival time in thyroid carcinoma [54,55].
By curating multiple independent GEO datasets comprising a total of 216 patient-derived PTC and PTC-free tissues, we generated a merged PTC-specific MMD annotated with clinical features across different subtypes of thyroid carcinomas, including anaplastic thyroid cancer, poorly-differentiated thyroid cancer, and papillary thyroid cancer. These unified data processed using a uniform R pipeline in this paper source would allow parallel cross-platform analyses with TCGA, as previously shown in other major types of cancer [30,31,32,33]. In this study, we further demonstrated the cross-platform compatibility of the newly developed PPARGi, providing promising results for its clinical application and our informatics pipeline, which can be readily translated to other sequencing platforms.
5. Conclusions
It remains to be investigated whether the identified intronic variants found in monozygotic twins in the LBD of the PPARG gene would induce loss-of-function effects on PPARγ or impaired PPARγ-RXRα signaling pathway. Further, our sample size was small and might only represent a small subset of PTC cohort in Korea, although the expression levels and prognostic performance of PPARGi-comprising genes were extensively validated using public transcriptomic databases. We thus aim to assess the predictive power of PPARGi in a larger validation patient cohort or in a prospectively conducted study. Altogether, the functional PPARGi personalized scoring system may represent a powerful and effective genomic tool to improve patient management in PTC.
Conceptualization, S.B.L., S.Y.K., and S.-M.K.; methodology, J.K., S.-X.M., S.-J.S., S.Y.K., and S.B.L.; software, J.K., S.-X.M., and S.B.L.; validation, J.K., S.-X.M., and S.B.L.; formal analysis, J.K., S.-X.M., and S.B.L.; investigation, S.B.L., S.Y.K., and S.-M.K.; resources, J.K., S.-X.M., S.Y.K., S.-M.K., S.-J.S., Y.S.L., H.C., H.-S.C., C.S.P., and S.B.L.; data curation, J.K., S.-X.M., and S.B.L.; writing—original draft preparation, J.K., S.-X.M., and S.B.L.; writing—review and editing, J.K., S.-X.M., S.Y.K., S.-M.K., and S.B.L.; visualization, J.K., S.-X.M., and S.B.L.; supervision, S.B.L., S.Y.K., and S.-M.K.; project administration, S.B.L., S.Y.K., and S.-M.K.; funding acquisition, S.B.L., S.Y.K., and S.-M.K. All authors have read and agreed to the published version of the manuscript.
This research was funded by the National Research Foundation of Korea (NRF) grant funded by the Korea government (MSIT) (2021R1F1A1064122 and 2017R1E1A1A03070345) and a faculty research fund of Ajou University School of Medicine and Yonsei University College of Medicine (6-2019-0092).
The study was conducted according to the guidelines of the Declaration of Helsinki and approved by the Institutional Review Board of Yonsei University (IRB 3-2020-0281) on 4 July 2021.
Patient consent was waived as this study was retrospective by design.
Datasets used to generate MMD-THCA are available at GEO under the accession code GSE29265, GSE33630, GSE65144, and GSE76039. The scRNA-seq dataset analyzed in this study is available at GEO under the accession code GSE158291.
This research was funded by the National Research Foundation of Korea (NRF) grant funded by the Korea government (MSIT) (2021R1F1A1064122 and 2017R1E1A1A03070345) and a faculty research fund of Ajou University School of Medicine and Yonsei University College of Medicine (6-2019-0092).
The authors declare no conflict of interest.
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Enlarge this image.Figure 1. Diagnostic and prognostic performance of PPARGi in TCGA-THCA. (A) Venn diagrams showing the number of the SnpEff-predicted variants in three pairs of twins (TP1–3). Genomic positions of mutations, including chromosome (chr) number, and gene symbols are shown in the table (right). (B) Expression heatmap (left) and forest plot (right) of PPARGi-comprising genes. The data are sorted in increasing PPARGi. The horizontal axis of the forest plot represents hazard ratio (HR) with 95% confidence intervals (CI) estimated using a Cox proportional hazards model. The regression coefficients (Coef) and the Wald statistic p-values (p) are stated. (C) PPARGi computation for patient stratification. (D) PPARGi of normal (n = 59) and tumor tissues (n = 505). Mann–Whitney–Wilcoxon test p-values (p) and the number of samples (n) are stated. (E) The area under the ROC curve (AUC) of the PPARGi classifier. The AUC value and the optimal cut-off are stated. (F) The hazard ratio (HR), log-rank p-value, and the number of patients successfully stratified (n) determined from univariate Cox regression analysis are shown on the survival Kaplan–Meier (KM) curve. Black and red KM curves represent predicted PPARGilow and PPARGihigh group, respectively. (G) Volcano plot depicting differentially expressed (red and blue) and non-significant (gray) genes in the PPARGi-stratified groups. The number of genes (n) are stated. (H) GSEA plot showing the enrichment of thyroid hormone generation gene set (GO:0006590) in PPARGilow tumors. The cumulative enrichment score (ES) is plotted as the green curve, which is the running sum of the weighted ES as the analysis walks down the limma-generated ranked list. The vertical black lines on the horizontal axis of the plot indicate the position of query genes in the ranked list of genes. The bottom plot shows the value of the fold change (log2-base) as the computation goes down the limma-generated ranked list. Normalized ES (NES) and adjusted p-values (padj) are stated. (I) Dot plot showing top gene sets (downregulated) in PPARGilow tumors. (J) Dot plot showing Pearson’s correlations between PPARGi and CIBERSORT-estimated proportion of immune cell populations.
Enlarge this image.Figure 2. PPARGi in MMD-THCA. (A) UMAP representation of the uncorrected (left) and corrected (right) MMD. Data are colored by study (top) and tissue type (bottom). The number of samples (n) are stated. (B) Disease-gene associations depicting the linkages of genes and the enriched diseases as a network. The color and size of the node represent the value of fold change (FC) and gene count, respectively. (C) Expression heatmap of PPARGi-comprising genes. The data are sorted in increasing PPARGi. (D) PPARGi of normal and tumor tissues (ATC, PDTC, and PTC). The asterisks represent the statistical significance assessed by Mann–Whitney–Wilcoxon test (**** p ≤ 0.0001, ** p ≤ 0.01, ns: p > 0.05). Kruskal–Wallis p-value (p) is stated. (E) The area under the ROC curve (AUC) of the PPARGi classifier. The AUC values are stated for ATC (left), PDTC (middle), and PTC (right). (F) Volcano plot depicting differentially expressed (red and blue) and non-significant (gray) genes in the PPARGi-stratified groups in PTC. The number of genes (n) are stated. (G) GSEA plot showing the enrichment of thyroid hormone generation gene set (GO:0006590) in PPARGilow tumors. The cumulative enrichment score is plotted as the green curve, which is the running sum of the weighted ES as the analysis walks down the limma-generated ranked list. The vertical black lines on the horizontal axis of the plot indicate the position of query genes in the ranked list of genes. The bottom plot shows the value of the fold change (log2-base) as the computation goes down the limma-generated ranked list. Normalized ES (NES) and adjusted p-values (padj) are stated. (H) Dot plot showing Pearson’s correlations between PPARGi and CIBERSORT-estimated proportion of immune cell populations.
Enlarge this image.Figure 3. Analysis of PPARGi-comprising genes at the single-cell level in PTC. (A) UMAP of 4045 QC-passed cells depicting 15 cell clusters. (B) Heatmap of top 5 differentially expressed features across the identified cell clusters. The number (n) and proportion (%) of cells in each cluster are stated. (C) Distribution of PPARGi across different cell clusters. Kruskal–Wallis p-value (p) is stated. (D) Dot plot depicting average expression of PPARGi-comprising genes across different cell clusters. (E) Heatmap showing expression levels of PPARGi-comprising genes and PPARGi in macrophages (left) and epithelial cells (right). The data are sorted in increasing PPARGi. (F) Dot plot showing top enriched GO terms from PPARGi-comprising genes.
Enlarge this image.Figure 4. ML algorithms applied to expression profiles of PPARGi-comprising genes. (A) The t-distributed stochastic neighbor embedding (t-SNE) visualization of cancers (TT) and normal (TT) tissues in TCGA-THCA (left), MMD-THCA (PTC, middle), and MMD-THCA (ATC, right). (B) Confusion matrices of LASSO and SVM models in TCGA-THCA (top right), MMD-PTC (bottom left), and MMD-ATC (bottom right). (C) ROC curves and (D) confusion matrix showing classifying performance for normal (left) and tumor (right) tissues in MMD-PTC using kNN model developed from TDM-transformed TCGA-PTC. (E) t-SNE visualization of PPARGihigh (colored in blue) and PPARGiow (colored in red) tumors in TCGA-THCA. (F) Confusion matrix of SVM model in classifying PPARGihigh from PPARGilow tumors. (G) ROC curve showing classifying performance for PPARGiow tumors.
ML models for disease selection in thyroid cancer.
| Dataset | Model | AUC | CA | F1 | Precision | Recall |
|---|---|---|---|---|---|---|
| TCGA-THCA | kNN | 0.980 | 0.961 | 0.962 | 0.963 | 0.961 |
| SVM | 0.989 | 0.968 | 0.968 | 0.967 | 0.968 | |
| Random forest | 0.975 | 0.956 | 0.954 | 0.954 | 0.956 | |
| Neural network | 0.986 | 0.975 | 0.976 | 0.976 | 0.975 | |
| Logistic regression | 0.990 | 0.972 | 0.972 | 0.972 | 0.972 | |
| MMD-THCA (PTC) | kNN | 0.925 | 0.884 | 0.884 | 0.885 | 0.884 |
| SVM | 0.937 | 0.912 | 0.911 | 0.912 | 0.912 | |
| Random forest | 0.914 | 0.871 | 0.871 | 0.871 | 0.871 | |
| Neural network | 0.928 | 0.912 | 0.911 | 0.913 | 0.912 | |
| Logistic regression | 0.929 | 0.912 | 0.912 | 0.912 | 0.912 | |
| MMD-THCA (ATC) | kNN | 0.944 | 0.838 | 0.835 | 0.841 | 0.838 |
| SVM | 0.929 | 0.862 | 0.861 | 0.861 | 0.862 | |
| Random forest | 0.897 | 0.808 | 0.807 | 0.807 | 0.808 | |
| Neural network | 0.938 | 0.862 | 0.862 | 0.862 | 0.862 | |
| Logistic regression | 0.945 | 0.892 | 0.892 | 0.892 | 0.892 |
Cross-platform evaluation results of the ML models.
| Dataset | Model | AUC | CA | F1 | Precision | Recall |
|---|---|---|---|---|---|---|
| Train-test data | kNN | 0.942 | 0.864 | 0.862 | 0.873 | 0.864 |
| SVM | 0.890 | 0.469 | 0.300 | 0.220 | 0.469 | |
| Random forest | 0.828 | 0.599 | 0.551 | 0.721 | 0.599 | |
| Neural network | 0.925 | 0.837 | 0.836 | 0.860 | 0.837 | |
| Logistic regression | 0.922 | 0.531 | 0.368 | 0.282 | 0.531 |
ML models for risk stratification in thyroid cancer.
| Dataset | Model | AUC | CA | F1 | Precision | Recall |
|---|---|---|---|---|---|---|
| TCGA-THCA | kNN | 0.965 | 0.980 | 0.979 | 0.979 | 0.980 |
| SVM | 0.998 | 0.990 | 0.990 | 0.990 | 0.990 | |
| Random forest | 0.988 | 0.966 | 0.962 | 0.964 | 0.966 | |
| Neural network | 0.997 | 0.982 | 0.982 | 0.982 | 0.982 | |
| Logistic regression | 0.974 | 0.976 | 0.975 | 0.975 | 0.976 |
Supplementary Materials
The following are available online at
References
1. Miranda-Filho, A.; Lortet-Tieulent, J.; Bray, F.; Cao, B.; Franceschi, S.; Vaccarella, S.; Dal Maso, L. Thyroid cancer incidence trends by histology in 25 countries: A population-based study. Lancet Diabetes Endocrinol.; 2021; 9, pp. 225-234. [DOI: https://dx.doi.org/10.1016/S2213-8587(21)00027-9]
2. Park, S.; Oh, C.-M.; Cho, H.; Lee, J.Y.; Jung, K.-W.; Jun, J.K.; Won, Y.-J.; Kong, H.-J.; Choi, K.S.; Lee, Y.J. et al. Association between screening and the thyroid cancer “epidemic” in South Korea: Evidence from a nationwide study. BMJ; 2016; 355, i5745. [DOI: https://dx.doi.org/10.1136/bmj.i5745] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/27903497]
3. Nixon, I.J.; Kuk, D.; Wreesmann, V.; Morris, L.; Palmer, F.L.; Ganly, I.; Patel, S.G.; Singh, B.; Tuttle, R.M.; Shaha, A.R. et al. Defining a valid age cutoff in staging of well-differentiated thyroid cancer. Ann. Surg. Oncol.; 2016; 23, pp. 410-415. [DOI: https://dx.doi.org/10.1245/s10434-015-4762-2] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/26215199]
4. Fugazzola, L.; Muzza, M.; Pogliaghi, G.; Vitale, M. Intratumoral genetic heterogeneity in papillary thyroid cancer: Occurrence and clinical significance. Cancers; 2020; 12, 383. [DOI: https://dx.doi.org/10.3390/cancers12020383]
5. Muzza, M.; Pogliaghi, G.; Persani, L.; Fugazzola, L.; Colombo, C. Combined mutational and clonality analyses support the existence of intra-tumor heterogeneity in papillary thyroid cancer. J. Clin. Med.; 2021; 10, 2645. [DOI: https://dx.doi.org/10.3390/jcm10122645] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/34208446]
6. Toraih, E.A.; Hussein, M.H.; Zerfaoui, M.; Attia, A.S.; Marzouk Ellythy, A.; Mostafa, A.; Ruiz, E.M.L.; Shama, M.A.; Russell, J.O.; Randolph, G.W. et al. Site-specific metastasis and survival in papillary thyroid cancer: The importance of brain and multi-organ disease. Cancers; 2021; 13, 1625. [DOI: https://dx.doi.org/10.3390/cancers13071625]
7. Khatami, F.; Tavangar, S.M. A review of driver genetic alterations in thyroid cancers. Iran. J. Pathol.; 2018; 13, pp. 125-135. [DOI: https://dx.doi.org/10.30699/ijp.13.2.125] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/30697281]
8. Raman, P.; Koenig, R.J. Pax-8–PPAR-γ fusion protein in thyroid carcinoma. Nat. Rev. Endocrinol.; 2014; 10, pp. 616-623. [DOI: https://dx.doi.org/10.1038/nrendo.2014.115] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/25069464]
9. Vu-Phan, D.; Grachtchouk, V.; Yu, J.; Colby, L.A.; Wicha, M.S.; Koenig, R.J. The thyroid cancer PAX8-PPARG fusion protein activates Wnt/TCF-responsive cells that have a transformed phenotype. Endocr. Relat. Cancer; 2013; 20, pp. 725-739. [DOI: https://dx.doi.org/10.1530/ERC-13-0058] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/24025583]
10. Eberhardt, N.L.; Grebe, S.K.; McIver, B.; Reddi, H.V. The role of the PAX8/PPARgamma fusion oncogene in the pathogenesis of follicular thyroid cancer. Mol. Cell. Endocrinol.; 2010; 321, pp. 50-56. [DOI: https://dx.doi.org/10.1016/j.mce.2009.10.013] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/19883731]
11. Hernandez-Quiles, M.; Broekema, M.F.; Kalkhoven, E. PPARgamma in metabolism, immunity, and cancer: Unified and diverse mechanisms of action. Front. Endocrinol.; 2021; 12, 624112. [DOI: https://dx.doi.org/10.3389/fendo.2021.624112] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/33716977]
12. Aprile, M.; Cataldi, S.; Ambrosio, M.R.; D’Esposito, V.; Lim, K.; Dietrich, A.; Bluher, M.; Savage, D.B.; Formisano, P.; Ciccodicola, A. et al. PPARgammaDelta5, a naturally occurring dominant-negative splice isoform, impairs PPARgamma function and adipocyte differentiation. Cell Rep.; 2018; 25, pp. 1577-1592. [DOI: https://dx.doi.org/10.1016/j.celrep.2018.10.035] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/30404011]
13. Lefterova, M.I.; Haakonsson, A.K.; Lazar, M.A.; Mandrup, S. PPARgamma and the global map of adipogenesis and beyond. Trends Endocrinol. Metab.; 2014; 25, pp. 293-302. [DOI: https://dx.doi.org/10.1016/j.tem.2014.04.001]
14. Nielsen, R.; Pedersen, T.A.; Hagenbeek, D.; Moulos, P.; Siersbaek, R.; Megens, E.; Denissov, S.; Borgesen, M.; Francoijs, K.J.; Mandrup, S. et al. Genome-wide profiling of PPARgamma:RXR and RNA polymerase II occupancy reveals temporal activation of distinct metabolic pathways and changes in RXR dimer composition during adipogenesis. Genes Dev.; 2008; 22, pp. 2953-2967. [DOI: https://dx.doi.org/10.1101/gad.501108]
15. Sarhangi, N.; Sharifi, F.; Hashemian, L.; Hassani Doabsari, M.; Heshmatzad, K.; Rahbaran, M.; Jamaldini, S.H.; Aghaei Meybodi, H.R.; Hasanzad, M. PPARG (Pro12Ala) genetic variant and risk of T2DM: A systematic review and meta-analysis. Sci. Rep.; 2020; 10, 12764. [DOI: https://dx.doi.org/10.1038/s41598-020-69363-7]
16. Le Menn, G.; Neels, J.G. Regulation of immune cell function by PPARs and the connection with metabolic and neurodegenerative diseases. Int. J. Mol. Sci.; 2018; 19, 1575. [DOI: https://dx.doi.org/10.3390/ijms19061575] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/29799467]
17. Szeles, L.; Torocsik, D.; Nagy, L. PPARgamma in immunity and inflammation: Cell types and diseases. Biochim. Biophys. Acta; 2007; 1771, pp. 1014-1030. [DOI: https://dx.doi.org/10.1016/j.bbalip.2007.02.005] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/17418635]
18. Moore, K.J.; Rosen, E.D.; Fitzgerald, M.L.; Randow, F.; Andersson, L.P.; Altshuler, D.; Milstone, D.S.; Mortensen, R.M.; Spiegelman, B.M.; Freeman, M.W. The role of PPAR-gamma in macrophage differentiation and cholesterol uptake. Nat. Med.; 2001; 7, pp. 41-47. [DOI: https://dx.doi.org/10.1038/83328] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/11135614]
19. Penas, F.; Mirkin, G.A.; Vera, M.; Cevey, A.; Gonzalez, C.D.; Gomez, M.I.; Sales, M.E.; Goren, N.B. Treatment in vitro with PPARalpha and PPARgamma ligands drives M1-to-M2 polarization of macrophages from T. cruzi-infected mice. Biochim. Biophys. Acta; 2015; 1852, pp. 893-904. [DOI: https://dx.doi.org/10.1016/j.bbadis.2014.12.019] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/25557389]
20. Chawla, A.; Barak, Y.; Nagy, L.; Liao, D.; Tontonoz, P.; Evans, R.M. PPAR-gamma dependent and independent effects on macrophage-gene expression in lipid metabolism and inflammation. Nat. Med.; 2001; 7, pp. 48-52. [DOI: https://dx.doi.org/10.1038/83336] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/11135615]
21. Welch, J.S.; Ricote, M.; Akiyama, T.E.; Gonzalez, F.J.; Glass, C.K. PPARgamma and PPARdelta negatively regulate specific subsets of lipopolysaccharide and IFN-gamma target genes in macrophages. Proc. Natl. Acad. Sci. USA; 2003; 100, pp. 6712-6717. [DOI: https://dx.doi.org/10.1073/pnas.1031789100] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/12740443]
22. Costa, V.; Gallo, M.A.; Letizia, F.; Aprile, M.; Casamassimi, A.; Ciccodicola, A. PPARG: Gene expression regulation and next-generation sequencing for unsolved issues. PPAR Res.; 2010; 2010, 409168. [DOI: https://dx.doi.org/10.1155/2010/409168] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/20871817]
23. Vohra, S.; Biggin, P.C. Mutationmapper: A tool to aid the mapping of protein mutation data. PLoS ONE; 2013; 8, e71711. [DOI: https://dx.doi.org/10.1371/journal.pone.0071711] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/23951226]
24. Tate, J.G.; Bamford, S.; Jubb, H.C.; Sondka, Z.; Beare, D.M.; Bindal, N.; Boutselakis, H.; Cole, C.G.; Creatore, C.; Dawson, E. et al. COSMIC: The catalogue of somatic mutations in cancer. Nucleic Acids Res.; 2019; 47, pp. D941-D947. [DOI: https://dx.doi.org/10.1093/nar/gky1015] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/30371878]
25. Colaprico, A.; Silva, T.C.; Olsen, C.; Garofano, L.; Cava, C.; Garolini, D.; Sabedot, T.S.; Malta, T.M.; Pagnotta, S.M.; Castiglioni, I. et al. TCGAbiolinks: An R/Bioconductor package for integrative analysis of TCGA data. Nucleic Acids Res.; 2016; 44, e71. [DOI: https://dx.doi.org/10.1093/nar/gkv1507] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/26704973]
26. Robinson, M.D.; McCarthy, D.J.; Smyth, G.K. edgeR: A Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics; 2010; 26, pp. 139-140. [DOI: https://dx.doi.org/10.1093/bioinformatics/btp616]
27. Ritchie, M.E.; Phipson, B.; Wu, D.; Hu, Y.; Law, C.W.; Shi, W.; Smyth, G.K. Limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res.; 2015; 43, e47. [DOI: https://dx.doi.org/10.1093/nar/gkv007] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/25605792]
28. Gauthier, N.P.; Reznik, E.; Gao, J.; Sumer, S.O.; Schultz, N.; Sander, C.; Miller, M.L. MutationAligner: A resource of recurrent mutation hotspots in protein domains in cancer. Nucleic Acids Res.; 2016; 44, pp. D986-D991. [DOI: https://dx.doi.org/10.1093/nar/gkv1132] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/26590264]
29. Bahceci, I.; Dogrusoz, U.; La, K.C.; Babur, O.; Gao, J.; Schultz, N. PathwayMapper: A collaborative visual web editor for cancer pathways and genomic data. Bioinformatics; 2017; 33, pp. 2238-2240. [DOI: https://dx.doi.org/10.1093/bioinformatics/btx149] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/28334343]
30. Lim, S.B.; Chua, M.L.K.; Yeong, J.P.S.; Tan, S.J.; Lim, W.-T.; Lim, C.T. Pan-cancer analysis connects tumor matrisome to immune response. NPJ Precis. Oncol.; 2019; 3, 15. [DOI: https://dx.doi.org/10.1038/s41698-019-0087-0] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/31123708]
31. Lim, S.B.; Tan, S.J.; Lim, W.T.; Lim, C.T. A merged lung cancer transcriptome dataset for clinical predictive modeling. Sci. Data; 2018; 5, 180136. [DOI: https://dx.doi.org/10.1038/sdata.2018.136] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/30040079]
32. Lim, S.B.; Tan, S.J.; Lim, W.T.; Lim, C.T. An extracellular matrix-related prognostic and predictive indicator for early-stage non-small cell lung cancer. Nat. Commun.; 2017; 8, 1734. [DOI: https://dx.doi.org/10.1038/s41467-017-01430-6] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/29170406]
33. Lim, S.B.; Tan, S.J.; Lim, W.T.; Lim, C.T. Compendiums of cancer transcriptomes for machine learning applications. Sci. Data; 2019; 6, 194. [DOI: https://dx.doi.org/10.1038/s41597-019-0207-2]
34. Gautier, L.; Cope, L.; Bolstad, B.M.; Irizarry, R.A. Affy—Analysis of Affymetrix GeneChip data at the probe level. Bioinformatics; 2004; 20, pp. 307-315. [DOI: https://dx.doi.org/10.1093/bioinformatics/btg405]
35. Leek, J.T.; Johnson, W.E.; Parker, H.S.; Jaffe, A.E.; Storey, J.D. The sva package for removing batch effects and other unwanted variation in high-throughput experiments. Bioinformatics; 2012; 28, pp. 882-883. [DOI: https://dx.doi.org/10.1093/bioinformatics/bts034] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/22257669]
36. Peng, M.; Wei, G.; Zhang, Y.; Li, H.; Lai, Y.; Guo, Y.; Chen, Y.; Liu, L.; Xiao, H.; Guan, H. et al. Single-cell transcriptomic landscape reveals the differences in cell differentiation and immune microenvironment of papillary thyroid carcinoma between genders. Cell Biosci.; 2021; 11, 39. [DOI: https://dx.doi.org/10.1186/s13578-021-00549-w] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/33588924]
37. Hao, Y.; Hao, S.; Andersen-Nissen, E.; Mauck, W.M., 3rd; Zheng, S.; Butler, A.; Lee, M.J.; Wilk, A.J.; Darby, C.; Zager, M. et al. Integrated analysis of multimodal single-cell data. Cell; 2021; 184, pp. 3573-3587. [DOI: https://dx.doi.org/10.1016/j.cell.2021.04.048]
38. Korotkevich, G.; Sukhov, V.; Budin, N.; Shpak, B.; Artyomov, M.N.; Sergushichev, A. Fast gene set enrichment analysis. bioRxiv; 2021; 060012. [DOI: https://dx.doi.org/10.1101/060012]
39. Yu, G.; Wang, L.-G.; Han, Y.; He, Q.-Y. ClusterProfiler: An R package for comparing biological themes among gene clusters. OMICS; 2012; 16, pp. 284-287. [DOI: https://dx.doi.org/10.1089/omi.2011.0118] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/22455463]
40. Robin, X.; Turck, N.; Hainard, A.; Tiberti, N.; Lisacek, F.; Sanchez, J.C.; Muller, M. pROC: An open-source package for R and S+ to analyze and compare ROC curves. BMC Bioinform.; 2011; 12, 77. [DOI: https://dx.doi.org/10.1186/1471-2105-12-77] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/21414208]
41. Lim, S.B.; Yeo, T.; Lee, W.D.; Bhagat, A.A.S.; Tan, S.J.; Tan, D.S.W.; Lim, W.T.; Lim, C.T. Addressing cellular heterogeneity in tumor and circulation for refined prognostication. Proc. Natl. Acad. Sci. USA; 2019; 116, pp. 17957-17962. [DOI: https://dx.doi.org/10.1073/pnas.1907904116] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/31416912]
42. Budczies, J.; Klauschen, F.; Sinn, B.V.; Gyorffy, B.; Schmitt, W.D.; Darb-Esfahani, S.; Denkert, C. Cutoff Finder: A comprehensive and straightforward Web application enabling rapid biomarker cutoff optimization. PLoS ONE; 2012; 7, e51862. [DOI: https://dx.doi.org/10.1371/journal.pone.0051862] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/23251644]
43. Fang, L.; Zhang, M.; Li, Y.; Liu, Y.; Cui, Q.; Wang, N. PPARgene: A database of experimentally verified and computationally predicted PPAR target genes. PPAR Res.; 2016; 2016, 6042162. [DOI: https://dx.doi.org/10.1155/2016/6042162] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/27148361]
44. Lee, S.; Seo, J.; Park, J.; Nam, J.Y.; Choi, A.; Ignatius, J.S.; Bjornson, R.D.; Chae, J.H.; Jang, I.J.; Lee, S. et al. Korean Variant Archive (KOVA): A reference database of genetic variations in the Korean population. Sci. Rep.; 2017; 7, 4287. [DOI: https://dx.doi.org/10.1038/s41598-017-04642-4] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/28655895]
45. Jang, D.M.; Jang, J.Y.; Kim, H.J.; Han, B.W. Differential effects of cancer-associated mutations enriched in helix H3 of PPARgamma. Cancers; 2020; 12, 3580. [DOI: https://dx.doi.org/10.3390/cancers12123580] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/33266062]
46. Yousefnia, S.; Momenzadeh, S.; Seyed Forootan, F.; Ghaedi, K.; Nasr Esfahani, M.H. The influence of peroxisome proliferator-activated receptor gamma (PPARgamma) ligands on cancer cell tumorigenicity. Gene; 2018; 649, pp. 14-22. [DOI: https://dx.doi.org/10.1016/j.gene.2018.01.018] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/29369787]
47. Korpal, M.; Puyang, X.; Jeremy Wu, Z.; Seiler, R.; Furman, C.; Oo, H.Z.; Seiler, M.; Irwin, S.; Subramanian, V.; Julie Joshi, J. et al. Evasion of immunosurveillance by genomic alterations of PPARgamma/RXRalpha in bladder cancer. Nat. Commun.; 2017; 8, 103. [DOI: https://dx.doi.org/10.1038/s41467-017-00147-w] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/28740126]
48. Ricote, M.; Glass, C.K. PPARs and molecular mechanisms of transrepression. Biochim. Biophys. Acta; 2007; 1771, pp. 926-935. [DOI: https://dx.doi.org/10.1016/j.bbalip.2007.02.013] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/17433773]
49. Heming, M.; Gran, S.; Jauch, S.L.; Fischer-Riepe, L.; Russo, A.; Klotz, L.; Hermann, S.; Schafers, M.; Roth, J.; Barczyk-Kahlert, K. Peroxisome proliferator-activated receptor-gamma modulates the response of macrophages to lipopolysaccharide and glucocorticoids. Front. Immunol.; 2018; 9, 893. [DOI: https://dx.doi.org/10.3389/fimmu.2018.00893]
50. Setoguchi, K.; Misaki, Y.; Terauchi, Y.; Yamauchi, T.; Kawahata, K.; Kadowaki, T.; Yamamoto, K. Peroxisome proliferator-activated receptor-gamma haploinsufficiency enhances B cell proliferative responses and exacerbates experimentally induced arthritis. J. Clin. Investig.; 2001; 108, pp. 1667-1675. [DOI: https://dx.doi.org/10.1172/JCI13202]
51. Lee, C.H.; Chawla, A.; Urbiztondo, N.; Liao, D.; Boisvert, W.A.; Evans, R.M.; Curtiss, L.K. Transcriptional repression of atherogenic inflammation: Modulation by PPARdelta. Science; 2003; 302, pp. 453-457. [DOI: https://dx.doi.org/10.1126/science.1087344] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/12970571]
52. Jiang, C.; Ting, A.T.; Seed, B. PPAR-gamma agonists inhibit production of monocyte inflammatory cytokines. Nature; 1998; 391, pp. 82-86. [DOI: https://dx.doi.org/10.1038/34184] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/9422509]
53. Kuenzli, S.; Saurat, J.H. Effect of topical PPARbeta/delta and PPARgamma agonists on plaque psoriasis. A pilot study. Dermatology; 2003; 206, pp. 252-256. [DOI: https://dx.doi.org/10.1159/000068897]
54. Kobawala, T.P.; Trivedi, T.I.; Gajjar, K.K.; Patel, D.H.; Patel, G.H.; Ghosh, N.R. Significance of interleukin-6 in papillary thyroid carcinoma. J. Thyroid Res.; 2016; 2016, 6178921. [DOI: https://dx.doi.org/10.1155/2016/6178921] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/27034885]
55. Liu, X.; Su, C.; Xu, J.; Zhou, D.; Yan, H.; Li, W.; Chen, G.; Zhang, N.; Xu, D.; Hu, H. Immunohistochemical analysis of matrix metalloproteinase-9 predicts papillary thyroid carcinoma prognosis. Oncol. Lett.; 2019; 17, pp. 2308-2316. [DOI: https://dx.doi.org/10.3892/ol.2018.9850] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/30675296]
© 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/). Notwithstanding the ProQuest Terms and Conditions, you may use this content in accordance with the terms of the License.
Details
; Soo Young Kim 2 ; Shi-Xun, Ma 3 ; Seok-Mo, Kim 4
; Su-Jin, Shin 5
; Lee, Yong Sang 4 ; Chang, Hojin 4 ; Hang-Seok Chang 4 ; Park, Cheong Soo 6 ; Su Bin Lim 1
1 Department of Biochemistry and Molecular Biology, Ajou University School of Medicine, Suwon 16499, Korea;
2 Department of Surgery, Ajou University School of Medicine, Suwon 16499, Korea;
3 Department of Neurology, Institute for Cell Engineering, Johns Hopkins University School of Medicine, Baltimore, MD 21205, USA;
4 Thyroid Cancer Center, Department of Surgery, Institute of Refractory Thyroid Cancer, Gangnam Severance Hospital, Yonsei University College of Medicine, Seoul 06273, Korea;
5 Department of Pathology, Gangnam Severance Hospital, Yonsei University College of Medicine, Seoul 06273, Korea;
6 CHA Ilsan Medical Center, Department of Surgery, Goyang-si 10414, Korea;