PPARG: A New Prognostic Biomarker and Correlates with Immune Infiltrates in Breast Invasive Carcinoma Based on Bioinformatics Analysis

International Clinical and Medical Case Reports Journal Research Article (ISSN: 2832-5788) Int Clinc Med Case Rep Jour (ICMCRJ) 2023 | Volume 2 | Issue 4 PPARG: A New Prognostic Biomarker and Correlates with Immune Infiltrates in Breast Invasive Carcinoma Based on Bioinformatics Analysis Yaohui Jiang†, Xiaoyu Zhou†, Yalin Ju, Yuxin Ma,Yujiao Dong, Shaobo Wu* Health Science Center, Xi’an Jiaotong University, Xi’an 710061, Shaanxi Province, PR China †These authors contributed equally to this work and should be considered co-first authors Citation: Yaohui Jiang, Xiaoyu Zhou, Yalin Jv, Yuxin Ma,Yujiao Dong, Shaobo Wu. PPARG: A New Prognostic Biomarker and Correlates with Immune Infiltrates in Breast Invasive Carcinoma Based on Bioinformatics Analysis. Int Clinc Med Case Rep Jour. 2023;2(4):1-15. Received Date: 30 January, 2023; Accepted Date: 03 February, 2023; Published Date: 05 February, 2023 *Corresponding author: Shaobo Wu. Health Science Center, Xi’an Jiaotong University, Xi’an 710061, Shaanxi Province, PR China Copyright: © Shaobo Wu, Open Access 2023. This article, published in Int Clinc Med Case Rep Jour (ICMCRJ) (Attribution 4.0 International), as described by http:// creativecommons.org/licenses/by/4.0/. 1. ABSTRACT PPARG (Peroxisome Proliferator Activated Receptor Gamma) is a Protein Coding gene. This gene encodes a member of the peroxisome proliferator-activated receptor (PPAR) subfamily of nuclear receptors. The protein encoded by this gene is PPAR-gamma and is a regulator of adipocyte differentiation. However, the relationship between PPARG expression, tumor immunity, and prognosis in Breast invasive carcinoma (BRCA) remain unclear. PPARG expression and its influence on tumor prognosis were analyzed by the ONCOMINE, Tumor Immune Estimation Resource (TIMER) and Kaplan-Meier plotter. The relationship between BRCA expression and tumor immunity were analyzed by TIMER and Gene Expression Profiling Interactive Analysis (GEPIA). PPARG expression in BRCA tissues correlated with prognosis. PPARG expression was significantly lower in several human cancers, including HCC, than in corresponding normal tissues. Moreover, we screened hub genes by constructing a protein‐protein interaction (PPI) network. These findings demonstrate that PPARG is a potential prognostic biomarker that determines cancer progression and correlated with tumor immune cells infiltration in BRCA. Keywords: PPARG; Biomarker; Immune infiltrates; Breast cancer 2. INTRODUCTION Breast cancer is the most frequently diagnosed cancer in women and ranks second among causes for cancer related death in women.[1] Breast cancer is the most frequent malignancy in women worldwide and is curable in ~70–80% of patients with early-stage, non-metastatic disease. Advanced breast cancer with distant organ metastases is considered incurable with currently available therapies. On the molecular level, breast cancer is a heterogeneous disease; molecular features include activation of human epidermal growth factor receptor 2 (HER2, encoded by ERBB2), activation of hormone receptors (oestrogen receptor and progesterone receptor) and BRCA mutations. Treatment strategies differ according to molecular subtype. Breast cancer is the second International Clinical and Medical Case Reports Journal Research Article (ISSN: 2832-5788) Int Clinc Med Case Rep Jour (ICMCRJ) 2023 | Volume 2 | Issue 4 highest incidence rate of cancer.[2] Therefore, a deeper understanding of the molecular mechanisms involved in BRCA pathogenesis is critical for developing new treatments to improve survival rates. The Genechip, a high‐throughput platform and effective technique for gene expression is widely used in multiple disease research areas, and can detect the expression of thousands of genes and the association between disease and genes.[3-5] Although several bioinformatical studies on BRCA have been reported in recent years, we could sieve through different target genes through analyzing distinct databases, which could assist us in further exploring and better studying the underlying mechanisms. In this study, we downloaded two gene expression profiles (GSE15852, GSE42568) from the Gene Expression Omnibus (GEO) database. A total of 207 breast tissue samples were included in this study, consisting of 60 normal breast tissue samples and 147 BRCA samples. Then, the GEO2R online tool, volcano map software, and Venn diagram online tool were applied for DEGs in the four data sets above. Furthermore, Gene Ontology (GO) analysis and Kyoto Encyclopedia of Gene and Genome (KEGG) pathway analysis were performed to analyze these DEGs’ molecular function (MF), cellular component (CC), biological process (BP), and KEGG pathways using the Database for Annotation, Visualization, and Integrated Discovery (DAVID) and KOBAS, respectively. We constructed a protein‐protein interaction (PPI) network with Cytoscape MCODE (Molecular Complex Detection) for identifying 11 hub genes with high connectivity degrees. Moreover, we verified the differential expression of 11 screened hub genes and prognosis for survival of patients with BRCA in the large‐ scale gene data sets regarding BRCA from The Cancer Genome Atlas (TCGA) database. As a result, only PPARG qualified for the prognosis target gene. Finally, we carried out gene set enrichment analysis (GSEA) for these three genes and potential molecular mechanisms in BRCA. In conclusion, this study provides some additional effective biomarkers for BRCA patients and contributes to the understanding of the molecular mechanisms of BRCA progression. Then, we comprehensively analyzed the expression of PPARG, its correlation with prognosis in different types of tumors including BRCA, and the status of different tumorinfiltrating immune cells based on expression of specific markers using the ONCOMINE, Kaplan-Meier plotter, Gene Expression Profiling Interactive Analysis (GEPIA), and Tumor Immune Estimation Resource (TIMER) databases. Our results shed light on the important role of PPARG in BRCA prognosis and provided an underlying mechanism that PPARG expression might modulate tumor immunity by regulating the infiltration of immune cells in BRCA. 3. MATERIALS AND METHODS 3.1 Microarray data information The GEO (http://www.ncbi.nlm.nih.gov/geo/) database is a free public database of microarrays and is used for gene expression data sets and platform records. The gene expression profiles of GSE15852 and GSE42568 were chosen from the GEO database. GSE15852 was based on the GPL96 [HG-U133A] Affymetrix Human Genome U133A Array, containing 43 breast cancer tissues and 43 normal breast tissues. GSE42568 was based on the GPL570 [HG-U133_Plus_2] Affymetrix Human Genome U133 Plus 2.0 Array, containing 104 breast cancer tissues and 17 normal breast tissues. The downloaded data was processed using the R package, then calibrated, standardized, and log2 transformation was performed on all data of gene expression. International Clinical and Medical Case Reports Journal Research Article (ISSN: 2832-5788) Int Clinc Med Case Rep Jour (ICMCRJ) 2023 | Volume 2 | Issue 4 3.2 Screening for DEGs DEGs between BRCA samples and normal endometrium samples were identified using the GEO2R web tool (http://www.ncbi.nlm.nih.gov/geo/geo2r/).[6] GEO2R is an interactive online tool, which compares two groups of samples in most GEO series to obtain genes with different expressions under the same experimental conditions. The DEGs with |log FC|≥ 1 and adjust P< 0.05 were considered as the cut‐off criteria. Then, we applied for Venn software online (http://bioinfogp.cnb.csic.es/tools/venny/) to obtain the common DEGs in all four in-dependent cohorts. 3.3 GO term and KEGG pathway enrichment analysis of DEGs GO analysis, as a common approach, was conducted to identify the unique biological properties of DEGs, including biological processes, cellular components, and molecular functions.[7] KEGG pathway enrichment analysis was performed to explore key pathways to initiation and progression of BRCA.[8] The GO annotation enrichment of DEGs were performed using DAVID Bioinformatics Resources 6.8 (https://david.ncifcrf.gov/),[9] which consisted of systematic and integrative functional annotation tools to integrate functional genomic annotations. The KEGG pathway analysis of DEGs was performed using the KOBAS 3.0 online analysis tool (http://kobas.cbi.pku.edu.cn/).[10] We regarded P<0.05 and counts more than 2 as a statistically significant difference and significant enrichment, respectively. 3.4 PPI network and module analysis The PPI network is essential to interpret the key genes and important gene modules in cancer development. The network analysis software Cytoscape (www.cytoscape.org; version 3.2.1)[11] was used to visualize networks and screen hub genes from these DEGs according to degrees. Subsequently, the MCODE (Molecular Complex Detection) app in Cytoscape software was used to screen modules within the PPI network with the default parameters “Degree Cutoff = 7,” “Node Score Cutoff = 0.2,” “K‐Core = 2,” and “Max.Depth = 100.” Finally, the interactions of all hit genes were constructed using String (https://string-db.org, version 11.0).[12] 3.5 Validation of the hub genes in TCGA database The expression of upregulated and downregulated hub genes was validated in the TCGA_BRCA data set after hub genes were determined from these four GEO datasets. Furthermore, we analyzed the association between expression of hub genes and biological processes/pathways and phenotypes using GSEA (GSEA v2.0, available online: http://www.broad.mit.edu/gsea/). 3.6 Statistical analysis The results of RNA seq v2 of unpaired samples were analyzed using an independent sample t test. Hub genes levels from G1 to G4 was studied using a two‐sided Student's test for two‐group comparisons and by one‐way analysis of variance, followed by a Bonferroni post hoctest, for multiple group comparisons. The Kaplan‐Meier (KM) curve was conducted to assess the association between the expression level of hub genes and survival time of patients with BRCA. All statistical analyses were performed with the GraphPad Prism 5 (GraphPad Software, Inc, La Jolla, CA). Data were presented as the mean ±standard deviation. Statistical significance was considered significant when P< 0.05. Gene expression data from the Oncomine database were analyzed using the P-values, fold changes, and ranks. Survival curves were produced by the Kaplan-Meier plots and GEPIA database. The correlation of gene expression was evaluated in the TIMER and GEPIA databases using Spearman’s correlation analysis. P-values <0.05 were considered as statistically significant. International Clinical and Medical Case Reports Journal Research Article (ISSN: 2832-5788) Int Clinc Med Case Rep Jour (ICMCRJ) 2023 | Volume 2 | Issue 4 3.7 CCL14 gene expression analysis The mRNA levels of PPARG in several cancers including BRCA were identified from the Oncomine database (https://www.oncomine.org/resource/login.html).[13] The threshold was determined as follows: fold change of 1.5, P-value of 0.001, and gene ranking of all. 3.8 Kaplan-Meier survival curve analysis Kaplan-Meier survival curve analysis was performed to assess the correlation between the expression of the 54,000 genes on the survival rates in 21 different cancers using more than 10,000 cancer samples, including 371 liver, 1,440 gastric, 3,452 lung, 2,190 ovarian, and 6,234 breast cancer samples. Kaplan-Meier plots (http://kmplot.com/analysis/) were used to analyze the relationship between CCL14 gene expression and survival rates in liver, gastric, breast, pancreatic, ovarian, and lung cancers based on the hazard ratios (HR) and log-rank P-values.[14] 3.9 TIMER analysis TIMER database was used to systematically analyze the tumor-infiltrating immune cells (TIICs) in 32 cancer types using more than 10,000 samples from The Cancer Genome Atlas (TCGA) (https://cistrome.shinyapps. io/timer/) database.[15] TIMER determines the abundance of tumor-infiltrating immune cells (TIICs) based on the statistical analysis of gene expression profiles.[16] We analyzed the association between the level of CCL14 gene expression and the abundance of infiltrating immune cells, including CD4+ T cells, CD8+ T cells, B cells, neutrophils, dendritic cells and macrophages based on expression of specific marker genes in different cancers including HCC. The marker genes used for analysis of tumor-infiltrating immune cells including T cells, B cells, TAMs, monocytes, M1 macrophages, M2 macrophages, natural killer (NK) cells, neutrophils, dendritic cells (DCs), T-helper (Th) cells, T-helper 17 (Th17) cells, follicular helper T (Tfh) cells, exhausted T cells, and Tregs were based on data from previous studies.[17] CCL14 gene was on the x-axis and related marker genes are on the y-axis. 3.10 GEPIA analysis The Gene Expression Profiling Interactive Analysis (GEPIA) database (http://gepia.cancer-pku.cn/index. html) was used to analyze the RNA sequencing expression data from 113 normal and 1109 tumor tissue samples from the TCGA and GTEx projects.[18] We also used GEPIA to generate survival curves and determine OS and DFS rates and their correlation to specific gene expression in 33 different types of cancer to further confirm the significantly correlated genes in the TIMER analysis. 4. RESULTS 4.1 Microarray data information and identification of DEGs in BRCA Data from each GEO data set was respectively analyzed using GEO2R online tools to screen DEGs (|log FC|≥1 and adjust P<0.05). In this study, 104 BRCA tissues and 17 normal breast tissues from four GEO data sets were involved. A total of 1,731 upregulated and 1,688 downregulated genes were filtered from GSE42568; 70 upregulated and 150 downregulated genes from GSE15852. The DEGs of each gene expression profile data were visualized using a volcano plot (Figure 1A). Subsequently, 52 overlapping upregulated genes and 146 overlapping downregulated genes were screened in the BRCA group compared with the normal control group using Venn software online (Figure 1B). International Clinical and Medical Case Reports Journal Research Article (ISSN: 2832-5788) Int Clinc Med Case Rep Jour (ICMCRJ) 2023 | Volume 2 | Issue 4 4.2 GO functional enrichment analysis Aiming at explore the biological functions of the integrated DEGs in BRCA, GO term enrichment analysis of the 198 DEGs (52 upregulated genes and 146 downregulated genes) was performed using the DAVID online analysis tool. In biological process terms, the DEGs were significantly enriched in multicellular organismal homeostasis (Figure 1C). For molecular functions, the DEGs were mainly enriched in sulfur compound binding (Figure 1D). In cellular component terms, the DEGs were mainly involved in the collagen−containing extracellular matrix (Figure 1E). The significant GO enrichment terms are shown in bubble maps using the ggplot2.R package with P< 0.05. 4.3 KEGG pathway analysis As shown in Figure 1F, KEGG pathway analysis of the 198 aberrant DEGs was conducted in BRCA using the KOBAS 3.0 online analysis tool, of which the top was Fatty acid degradation. The analysis indicated that the 122 dysregulated genes were mainly related to human cancer. 4.4 PPI network construction, module analysis, and hub gene selection The PPI network was constructed for further investigation of the interaction among the integrated DEGs via Cytoscape software, based on the STRING database. The PPI network consisted of 155 nodes and 461 edges. We regarded the top 11 DEGs with a high degree of connectivity as the hub genes of BRCA (Figure 1G). The identified hub genes were all downregulated (PCK1, FABP4, LPL, CEBPA, CIDEC, ADIPOQ, CFD, LEP PLIN1, CIDEA and PPARG) genes. The String database was used to construct an interactive network of all hit genes. Interestingly, all the hub genes formed a network of interactions (Figure 1H). 4.5 Hub gene validation To further evaluate the 11 candidate hub genes that were dysregulated in BRCA, publicly available transcriptome data were downloaded from the TCGA_BRCA database. We analyzed RNA‐Seq profiles of 113 normal and 1109 cancer samples from the TCGA_BRCA data set. All hit genes were remarkably differently expressed in BRCA compared with control group, their expression tendency in TCGA data set was consistent with the messenger RNA (mRNA) level from GEO data set (Figure 2A‐K). Then we download pathological section of BRCA from human protein atlas (https://www.proteinatlas.org/search/HAMP) and we found that the PPARG expression products level was obviously decreased in breast cancer tissues compared with normal lung tissues (Figure 1I). 4.6 Prognostic significance of hubgenes in patients with BRCA To further elucidate whether these hub genes were potentially prognostic markers for BRCA, we analyzed the overall survival (OS) and disease‐free survival (DFS) for each hub gene by an online KM survival analysis tool (KMplot, http://kmplot.com/analysis/) and the patients were split by auto select best cutoff. Among those 11 hub genes, only the downregulated gene PPARG was associated with BRCA prognosis. The low expression of PPARG mRNA level was significantly associated with worse OS and DSS in patients with BRCA (Figures 3A and 3B). Therefore, PPARG could be prognostic markers for BRCA. 4.7 Clinical stage analyses To clarify the clinical significance of PPARG expression in BRCA, we analyzed their expression in different stages from TCGA. As shown in Figure 3F, increased PPARG expression were significantly associated with advanced tumor stage. We observed PPARG could distinguish patients with poor prognosis for DSS, PFI International Clinical and Medical Case Reports Journal Research Article (ISSN: 2832-5788) Int Clinc Med Case Rep Jour (ICMCRJ) 2023 | Volume 2 | Issue 4 (Figure 3C and 3E). These results demonstrated that might be associated with BRCA progression and might have prognostic significance for BRCA patients. 4.8 Hub gene GSEA analysis To further identify the possible mechanism of these hub genes in BRCA, GSEA was conducted to obtain the biological pathway from a database to a gene set. The patients from the TCGA-BRCA data set were divided into high and low-expression groups based on the median value of PPARG. As shown in Figure 2L, GSEA analysis revealed that the low expression of these hub genes was enriched in the “cytokine receptor interaction” and “cell surface interactions at the vascular” genesets from TCGA_BRCA data sets. 4.9 The levels of PPARG expression correlate with the infiltration levels of immune cells in BRCA The survival times of patients in several cancers is determined by the quantity and activity status of tumor infiltrating lymphocytes. Therefore, we analyzed the association of PPARG expression with prognosis and immune infiltration in BRCA. The analysis of immune infiltration by genomic methods is highly influenced by tumor purity in clinical samples. Moreover, GEPIA and TIMER databases contain most of the homologous TCGA. Therefore, we selected BRCA in the TIMER database that show significant correlation between PPARG expression and tumor purity, and in the GEPIA database that PPARG expression is relevant to tumor prognosis. We observed that low PPARG expression correlated with poorer prognosis and high infiltration of most immune cell types in BRCA (Figure 4). The level of PPARG expression negatively correlated with the infiltration levels of B cells (r=-0.039, P=2.28e-01), CD4+ T cells (r=0.25, P=3.18e-15), CD8+ T cells (r=0.279, P=5.96e-19), macrophages (r=-0.226, P=2.10e-17), neutrophils (r=0.176, P=4.75e-08), and DCs (r=-0.186, P=8.11e-09) in BRCA tissues. The data suggest that CCL14 expression modulates infiltration of immune cells into tumor tissues. 4.10 Correlation analysis between mRNA levels of PPARG and markers of different subsets of immune cells Next, we investigated the correlation between PPARG expression and the status of tumor-infiltrating immune cells based on the levels of immune marker gene expression in BRCA tissues using the TIMER and GEPIA databases. The immune cells analyzed in BRCA tissues included CD8+ T cells, T cell (general), Monocyte, tumor-associated macrophages (TAMs), monocytes, M1 and M2 macrophages, neutrophils, DCs, and natural killer (NK) cells. Moreover, different subsets of T cells, namely, T helper 1 (Th1), Th2, follicular helper T (Tfh), Th17, regulatory T (Tregs), and exhausted T cells were also analyzed. Since tumor purity of clinical samples influences the analysis of immune infiltration, the correlation analysis was adjusted for purity (Table 1). Analysis of the TIMER and GEPIA databases showed that PPARG expression in BRCA tissues significantly correlated with the expression of marker genes from tumor-infiltrating M1 and M2 Macrophage, Monocyte, Th1, Tfh, Neutrophils, B cell, CD8+ T cell (Figure 5 and Table 2). International Clinical and Medical Case Reports Journal Research Article (ISSN: 2832-5788) Int Clinc Med Case Rep Jour (ICMCRJ) 2023 | Volume 2 | Issue 4 Figure 1: Identification of overlapping DEGs, functional analysis and interaction network construction in BRCA. A, Volcano plots of gene expression profile data in GSE15852 and GSE42568, respectively. B, Venn plots of upregulated overlapping DEGs and Venn plots of downregulated overlapping DEGs. C, Biological process of the DEGs. D, Molecular function of the DEGs. E, Cell component of the DEGs. F, KEGG pathway of the DEGs. G, The PPI network of the DEGs. H, Top 11 hub genes PPI network constructed by STRING online database. I, Immunohistochemical staining of PPARG was performed in breast cancer and normal tissue. DEG, differentially expressed gene; BRCA, Breast invasive carcinoma; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Gene and Genome; PPI, protein‐protein interaction. International Clinical and Medical Case Reports Journal Research Article (ISSN: 2832-5788) Int Clinc Med Case Rep Jour (ICMCRJ) 2023 | Volume 2 | Issue 4 Figure 2: Expression validation of nine hub targets in BRCA compared with adjacent tissues from TCGA data sets. A, ADIPOQ. B, CEBPA. C, CFD. D, CIDEA. E, CIDEC. F, FABP4. G, LEP. H, LPL. I, PCK1. J, PLIN1. K, PPARG. L, Enrichment plots from GSEA. “cytokine receptor interaction” and “cell surface interactions at the vascular” were significantly enriched. BRCA, Breast invasive carcinoma; TCGA, The Cancer Genome Atlas; GSEA, Gene Set Enrichment Analysis. Figure 3: Prognostic gene signature in the BRCA patients. A, Association between PPARG expression and overall survival (OS) were performed with the online Kaplan‐Meier survival analysis (KMplot, http://kmplot.com/analysis/) (P<0.05). Low expression of PPARG was significantly associated with shorter OS time of BRCA patients. B, Online Kaplan-Meier survival analysis revealed a shorter DSS time with low PPARG expression level .C, the mRNA levels of the potential prognostic markers (PPARG) is extracted from the TCGA_BRCA data set and mRNA levels between alive and dead patients for PFI, OS(D) and DFS(E). F, The clinical significance of PPARG mRNA with T stage in TCGA. DSS, Disease Specific Survival; PFI, platinumfree interval; BRCA, Breast invasive carcinoma; mRNA, messenger RNA; OS, overall survival; TCGA, The Cancer Genome Atlas. International Clinical and Medical Case Reports Journal Research Article (ISSN: 2832-5788) Int Clinc Med Case Rep Jour (ICMCRJ) 2023 | Volume 2 | Issue 4 Figure 4: Correlation analysis of PPARG expression and infiltration levels of immune cells in BRCA tissues using the TIMER database. PPARG expression in BRCA tissues negatively correlates with tumor purity and infiltration levels of B cells, and positively correlates with CD8+ T cells, CD4+ T cells, macrophages, neutrophils, and dendritic cells.