Screening of potential biomarkers for polycystic ovary syndrome and identification of expression and immune characteristics.

(1) Data sources. Granulosa cell expression profiles of PCOS were acquired from Gene Expression Omnibus (GEO) database: GSE137684 , GSE34526 , GSE138518 , GSE168404 , GSE155489 ( Table 1 ) [ 17 ]. RNA-Seq expression profiles ( GSE107746 ) of granulosa cells at five follicular development stages were acquired from GEO database. (2) Data preprocessing. GSE137684 and GSE34526 (microarray) were combined as a training set. GSE138518 , GSE168404 , and GSE155489 (RNA-Seq) were combined as a test set, and the batch effect was removed using the ComBat function of the sva package (Figure S1) [ 18 ]. RNA-Seq expression profiles were converted into TPM format data. We used limma package to identify DEGs, and the threshold was set to p log2(1.5) [ 19 ]. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses were performed on all upregulated DEGs in PCOS using the R software clusterProfiler package to assess their pathway processes [ 20 ], biological processes (BP), cellular components (CC), and molecular functions (MF). Gene set enrichment analysis (GSEA) was performed to further assess the enrichment of relevant functional pathways between different groups [ 21 ]. The recursive feature elimination (RFE) algorithm of the caret package is used for characteristic gene selection and 10x cross-validation [ 22 ]. The key parameters we chose are: functions = "rfFuncs", method = "cv", number = 10. We then used LASSO regression of the glmnet package for characteristic gene screening [ 23 ]. The parameters chosen for the LASSO regression analysis were: nfolds = 5, family = ’binomial’, type. Measure = ’deviance’. The potential biomarkers were acquired by comprehensive analysis of the two algorithms. To further assess the diagnostic performance of candidate biomarkers, the receiver operating characteristic (ROC) curve and its area under the curve (AUC), accuracy, sensitivity, and specificity were computed using the R software timeROC package and confusion matrix [ 24 ]. The gene of Insulin secretion pathway was downloaded from KEGG, and the IR score of each sample was computed by ssGSEA method [ 25 ]. The MCP counter algorithm was used to estimate the relative proportion of 10 kinds of immune cells [ 26 ]. MCP counter is an algorithm that allows stable quantification of the absolute abundance of 2 stromal and 8 immune cell populations in heterogeneous tissues from transcriptomic data [ 27 ]. The infiltration of 28 immune cells was represented by an enrichment score, which was computed using a single sample gene set enrichment analysis (ssGSEA) of the GSVA package [ 21 ]. We acquired the ssGSEA scores of each sample in the HALLMARK Pathway, and computed the correlation between the biomarkers and the scores of each pathway. In addition, we obtained 103 immunomodulation-related genes based on previous studies and evaluated the relationship between PCOS-related genes and these genes using the spearman method [ 28 ]. The difference between the two groups was tested by Wilcoxon rank sum test [ 29 ]. Kruskal-Wallis test was used for the difference between the three groups and more [ 30 ]. Correlation analysis was computed using pearson method [ 31 ]. The P less than 0.05 is defined as a significant difference. Notably, * means p-value less than 0.05; ** means p-value less than 0.01; *** means p-value less than 0.001, and **** means p-value less than 0.0001. ns means there is no significant difference between two groups.

Polycystic ovary syndrome (PCOS) is a widespread disorder of endocrine and metabolic abnormalities that affects up to 15% of women of adolescent and gestational age [ 1 ]. PCOS can lead to abnormalities in the female reproductive system and immune system, and the clinical manifestations are generally irregular menstruation and even amenorrhea, excessive male hormones, polycystic ovary morphology [ 2 ]. Severe PCOS can lead to insulin resistance (IR), female anovulatory infertility, pathological obesity, autoimmune diseases, local or whole body inflammation and many other diseases [ 3 , 4 ]. The specific causes of PCOS are complex and multi-systemic, lack of targeted treatment, different diagnostic criteria and methods, and need to be adjusted according to different clinical symptoms [ 5 ]. In addition, it has been shown that follicular excess and increased size of the ovaries in PCOS are associated with insulin resistance (IR) and that the number of follicles may also be a predictor of IR [ 6 ]. However, the pathogenesis of PCOS is not yet perfect, and there are many hypotheses about different pathogenesis factors, such as chronic low-grade inflammation, genetic factors, lifestyle, post-translational modification mechanism of PCOS, and ovarian autophagy mechanism [ 7 – 9 ]. Ovarian granulosa cells (GC) are the main constituent cells of follicles, which can divide and differentiate in different degrees with the growth of follicles, and also play an important role in hormone secretion [ 10 ]. Abnormalities in GC can lead to ovary-related diseases, such as a deficiency of Rpg5 (an autophagy protein) in GC that can lead to premature ovarian failure [ 11 ]. Inflammation caused by high levels of male hormones can also lead to GC pyroptosis, causing ovarian function decline and tissue fibrosis [ 12 ]. The abnormal communication of oocyte and follicular cells and the damage of GC may be the key to PCOS, such as the damage of transregional projection in communication [ 13 ]. Therefore, understanding the abnormal mechanism of GC is very important for the study of PCOS. The purpose of this study was to search for potential biomarkers of PCOS and explore the association between biomarkers and immune status in PCOS. In this research, Recursive Feature Elimination (RFE) algorithm was used for selection of feature gene. The prediction strength and practicability of this algorithm are relatively good, and it has been used many times in the screening of biomarkers for different diseases, such as chronic obstructive pulmonary disease, keloid disorder, and triple negative breast cancer [ 14 – 16 ].

DEGs in PCOS group and normal control group were identified in microarray cohort and RNA-Seq cohort. A total of 1044 up-regulated genes and 695 and down-regulated genes were acquired in the microarray cohort ( Fig 1A ). A total of 464 up-regulated genes and 103 down-regulated genes were acquired in the RNA-Seq cohort ( Fig 1B ). Cross analysis between the two cohorts revealed 93 common DEGs ( Fig 1C ). A: PCOS data of microarray and DEGs volcano map of control group; B: PCOS data of the RNA-Seq cohort and DEGs volcano map of the control group; C: 93 common differential genes in the two datasets; D: GO and KEGG enrichment analysis of up-regulated genes in all PCOS data. By KEGG enrichment analysis, Tuberculosis, Epstein-Barr virus infection, Phagosome and other pathway processes were found to be significant. Neutrophil activation, neutrophil activation involved in immune response, neutrophil degranulation, neutrophil mediated immunity is a process with significant enrichment in BP. Secretory granule membrane was the most significantly enriched component in CC. Actin binding and GTPase binding were the two most enriched functions in MF. ( Fig 1D ). These results suggest that PCOS is primarily associated with immune-related pathways and biological functions in patients. In addition, GSEA analysis indicated that the enriched gene set was closely related to antigen processing and presentation, B cell receptor signaling pathway, and neutrophil-mediated immunity in the PCOS group ( Fig 2A and 2C ). Metabolism-related pathway activity was significantly higher in the normal control group ( Fig 2B ). A, C: GSEA showed immune signaling pathways and biological processes in the PCOS group; B, D: The 5 pathways most associated with the control group. RFE and LASSO regression were used to screen 93 differential genes. Three genes were acquired by RFE algorithm: JDP2, HMOX1 and BNIP3 ( Fig 3A ). Nine genes were acquired by LASSO regression: TP53I11, MXRA8, JDP2, CDH3, STC1, HMOX1, COLQ, LTBP1 and LGMN ( Fig 3B and 3C ). JDP2 and HMOX1 are the common genes acquired by the two algorithms. A: Screening biomarkers based on RFE algorithm; B, C: Characteristic gene selection by LASSO regression; D: ROC curve of JDP2 and HMOX1 in training set; E: ROC curve of the training set combination JDP2 and HMOX1; F: Training set model performance evaluation; G: ROC curve of JDP2 and HMOX1 in test set; H: ROC curve of the test set combination JDP2 and HMOX1; I: Test set model performance evaluation. ROC analysis indicated that JDP2 and HMOX1 were of high value in distinguishing PCOS from normal samples, and the AUC values of the training set and the test set were ideal ( Fig 3D and 3G ). The ROC curve after the combination of the two genes still had a high AUC value, with the training set AUC = 0.929 and the test set AUC = 0.917( Fig 3E and 3H ). Confusion matrix results showed that the Sensitivity of the training set and the test set were both greater than 0.8, and the Accuracy was greater than 0.9. In addition, we demonstrated using PCA analysis that JDP2 and HMOX1 were able to distinguish between PCOS samples and normal samples ( S1 Fig ). These results indicating that the two biomarkers identified had good classification performance ( Fig 3F and 3I ). Next, we explored the differences in the expression of two genes in PCOS and normal controls. In both the training set and the test set, JDP2 and HMOX1 showed high expression in the PCOS group, P < 0.001 ( Fig 4A and 4B ). At the same time, the specific expression of markers at different developmental stages was explored. Here, their expression in granulosa cells at different follicular development stages was compared by dataset GSE107746 . The results indicated that JDP2 and HMOX1, which should be inhibited during the pre-ovulation phase, were elevated in PCOS tissues, suggesting that high levels of JDP2 and HMOX1 expression in granulosa cells may also play an important role in follicular development in PCOS ( Fig 4C ). A: The expression of JDP2 and HMOX1 in the training set; B: The expression of JDP2 and HMOX1 in the test set; C: Expression of JDP2 and HMOX1 in granulosa cells at different follicular development stages; D: The difference of IR score between PCOS and normal control group; E: JDP2 and HMOX1 were correlated with IR scores. Insulin resistance (IR) plays an important role in the development of PCOS, and the IR score for each sample was computed using the ssGSEA method. The results indicated that PCOS had a higher IR score ( Fig 4D ). The expressions of JDP2 and HMOX1 showed a significant positive correlation with IR score, and the correlation of HMOX1 was stronger, P = 0.00076 ( Fig 4E ). These results demonstrate the presence of higher levels of insulin resistance in PCOS patients and that patients with higher IR scores will have high expression of JDP2 and HMOX1. We compared 28 kinds of immune cell scores between PCOS and normal control samples, and the results indicated that 21 kinds of immune cell scores were significantly higher in PCOS ( Fig 5A , S1 Table ). The immune infiltration results acquired by MCP Counter algorithm are consistent with the above, and the scores of B lineage, Monocytic lineage, and Neutrophils are higher than those of PCOS ( Fig 5B ). Correlation analysis was performed to assess the relationship between biomarkers and infiltrating immune cells. JDP2 showed significant positive correlation with T cells and cytotoxic lymphocytes, with P values less than 0.05 ( Fig 5C ). HMOX1 showed significant positive correlation with monocytes and myeloid dendritic cells. The maximum value of P is 0.036 ( Fig 5D ). In addition, we further explored the relationship of PCOS-related genes with chemokine-related genes, immunoinhibitor-related genes, immunostimulator-related genes, MHC-related genes and receptor-related genes. As shown in S2 Fig , we found that PCOS-associated genes were significantly associated with most of the immunomodulatory genes, such as CCL4, CXCL17, CSF1R, CD48, IL6R, HLA, CCR1 and CCR5. A: The immune scores between PCOS and control group were compared by ssGSEA method. B: MCP Counter method was used to compare the immune scores between PCOS and control group. C: The correlation between JDP2 and immune cells; D: Correlation between HMOX1 and immune cells. We computed the ssGSEA score of each sample in the HALLMARK Pathway, and computed the correlation between JDP2 and HMOX1 and the score of each pathway, respectively. In JDP2, INTERFERON ALPHA RESPONSE, INTERFERON GAMMA RESPONSE, ALLOGRAFT REJECTION and COMPLEMENT are four pathways with significant positive correlation ( Fig 6A ). In HMOX1, REACTIVE OXYGEN SPECIES PATHWAY, IL2 STAT5 SIGNALING, PI3K AKT MTOR SIGNALING, and IL6 JAK STAT3 SIGNALING and COMPLEMENT are the four pathways with significant positive correlation ( Fig 6B ). Comprehensive analysis indicated that the two gene organisms had strong correlation with immune-related pathways. A: The correlation between JDP2 and pathway score; B: The correlation between HMOX1 and pathway score.

93 DEGs were identified in PCOS, and enrichment analysis revealed that the set of immune-related genes was differentially activated in PCOS. Recursive feature elimination (RFE) and LASSO regression were used for further feature selection of these genes, and finally JDP2 and HMOX1 were acquired as candidate biomarkers for PCOS.

JDP2 and HMOX1 could distinguish between PCOS (A) and normal (B) samples in both the training and validation cohorts. (PDF) Click here for additional data file. Relationship of PCOS-related genes to chemokine (A), immunoinhibitory (B), immunostimulatory (C), MHC (D) and receptor-related genes (E). (PDF) Click here for additional data file. (CSV) Click here for additional data file.