Skip to main content
Download PDF
  • Research article
  • Open access
  • Published:

Construction of a ceRNA network of hub genes affecting immune infiltration in ovarian cancer identified by WGCNA

BMC Cancer volume 21, Article number: 970 (2021) Cite this article

Abstract

Background

Ovarian cancer is the leading cause of death among gynecological malignancies. Immunotherapy has demonstrated potential effects in ovarian cancer. However, few studies on immune-related prognostic signatures in ovarian cancer have been reported. This study aimed to identify hub genes associated with immune infiltrates to provide insight into the immune regulatory mechanisms in ovarian cancer.

Methods

Raw data and clinical information were downloaded from The Cancer Genome Atlas (TCGA) and University of California, Santa Cruz (UCSC) Xena websites. Single-sample gene set enrichment analysis (ssGSEA) and weighted gene co-expression network analysis (WGCNA) were used to identify hub genes. Kaplan-Meier analysis and differential expression analysis were applied to explore the real hub genes.

Results

Through ssGSEA and WGCNA, 7 hub genes (LY9, CD5, CXCL9, IL2RG, SLAMF1, SLAMF6, and SLAMF7) were identified. Finally, LY9 and SLAMF1 were recognized as the real hub genes in immune infiltrates of ovarian cancer. LY9 and SLAMF1 are classified as SLAM family receptors involved in the activation of hematopoietic cells and the pathogenesis of multiple malignancies. Furthermore, 12 lncRNAs and 43 miRNAs significantly related to the 2 hub genes were applied to construct a lncRNA-miRNA-mRNA ceRNA network. The lncRNA-miRNA-mRNA ceRNA network shows upstream regulatory sites of the 2 hub genes.

Conclusions

These findings improve our understanding of the regulatory mechanism of and reveal potential immune checkpoints for immunotherapy for ovarian cancer.

Peer Review reports

Background

Ovarian cancer is the leading cause of death among gynecological neoplasms [1]. Cytoreductive surgery followed by platinum-based chemotherapy is the standard treatment for advanced ovarian cancer. Unfortunately, most patients experience recurrence and eventually die. Angiogenesis inhibitors such as bevacizumab and PARP inhibitors have been recently added to the available treatment schemes. In recent years, there has been increasing interest related to the role played by immunotherapy in ovarian cancer progression control.

Unlike chemotherapy or other targeted therapies, immunotherapies have the advantage of triggering the immune response, which clinically exerts specific, systemic, and durable antitumor effects. However, less than 15% of patients with advanced/metastatic ovarian cancer respond to immune checkpoint inhibitors [2, 3]. In other words, this toxic and costly treatment is potentially ineffective in the majority of ovarian cancer patients. Therefore, biomarkers affecting the immune response are needed to guide treatment decisions. Expression of programmed death ligand 1 (PD-L1) on the surface of T cells [4], the microsatellite instability (MSI) status [5] or DNA mismatch repair deficiency (dMMR) [6] and total tumor mutational burden (TMB) [7, 8] are the main predictive markers for the response to immunotherapy. However, because of the complexity of tumor-immune interactions, efforts to capture this complexity via a single analyte, such as PD-L1 expression or tumor mutational load, as a surrogate of potential tumor antigenicity, yield limited and incomplete information about the complex and dynamic nature of the tumor-immune microenvironment [9]. In summary, many efforts have been made to explore tumor immune signatures; however, few efforts to regulate tumor immune signatures have been reported.

With the development of high-throughput sequencing, bioinformatics has played an increasingly important role in the field of cancer research. Many studies consider only differences in the expression of genes between different samples and ignore the underlying connection of each gene. Weighted gene co-expression network analysis (WGCNA) [10] is a systematic biological method used to describe correlation patterns among genes in samples and can identify clusters (modules) of highly correlated genes for the investigation of clinical traits. In the present study, we performed WGCNA on RNA-seq data derived from The Cancer Genome Atlas (TCGA), reconstructed a gene (green) module related to immune infiltrates and identified hub genes in this green module.

To further reveal the regulatory mechanism of immune infiltrate-related hub genes, we constructed a lncRNA-miRNA-mRNA competing endogenous RNA (ceRNA) network. The ceRNA hypothesis states that a pool of long noncoding RNAs (lncRNAs), circular RNAs (circRNAs) and messenger RNAs (mRNAs) compete and bind to microRNAs (miRNAs), regulating their activity [11, 12]. Among the ceRNAs, miRNAs regulate the expression of their target genes by binding to the miRNA response elements on the target mRNAs, and lncRNAs act as molecular sponges to repress the negative regulation of target mRNAs by miRNAs. An R package called multiMiR was used to predict miRNAs for hub genes, and Star7Base [13] was used to predict interactions between lncRNAs and miRNAs. Finally, we identified 12 lncRNAs, 43 miRNAs and 2 mRNAs to construct a lncRNA-miRNA-mRNA ceRNA network.

Overall, in our study, we aimed to screen hub genes associated with immune infiltrates of ovarian cancer using the WGCNA method, construct a ceRNA network, and provide insight into the regulatory mechanisms of immune infiltrates in ovarian cancer.

Methods

TCGA data download

Level 3 HTSeq-Counts data, HTSeq-FPKM (per million fragments mapped) data and level 1 clinical information (shown in Table 1), such as age, histological type, survival and outcome of patients with ovarian cancer, were downloaded from The Cancer Genome Atlas (TCGA; http://cancergenome.nih.gov/). We chose “cystic, mucinous and serous neoplasms” in “Disease Type” as our research object. There were 379 samples from 376 ovarian cancer patients included in the TCGA program. Details on these samples are listed in Table 1. Because the TCGA database does not contain information on normal ovarian tissue, we also downloaded combined TPM (transcripts per kilobase million) data, including non-diseased ovarian tissue downloaded from the GTEx (Genotype-Tissue Expression) project, and ovarian cancer tissue from TCGA, and data from the UCSC Xena website (https://xena.ucsc.edu/), to compare the expression values of genes.

Table 1 Clinical information about TCGA-OV
Full size table

Calculation of the Immunophenoscore

Data were analyzed in the R programming environment, version 4.0.2. The names of genes were converted from the Ensembl ID to the gene symbol through the AnnotationDbi and org. Hs.eg.db packages. The immunophenoscore (normalized enrichment score, NES) of each TCGA OV sample was calculated through the ssGSEA (single-sample gene set enrichment analysis) method using the GSVA package based on the expression of the representative genes in gene sets, which were downloaded from Cell Reports (https://doi.org/10.1016/j.celrep.2016.12.019) [14]. ssGSEA ranks the genes by their expression in each sample and computes the enrichment score by integrating the differences between the empirical cumulative distribution functions of the gene ranks [15]. The results were visualized with the ggplot2 R package.

Survival analysis

Univariate Cox regression analysis was performed to evaluate the association between the overall survival time and research objects, including the genes and immunophenoscores obtained above, using the survival and survminer R packages. P < 0.05 was considered to indicate a statistically significant difference. In addition, the log2(HR), 95% CI and statistical significance of the infiltrated immune cell types were calculated and illustrated using a forest plot through the forest plot package.

The patients were dichotomized based on the best cutoff of the immunophenoscore using the surv_cutpoint function of the survminer R package, which divides the immunophenoscore into high and low groups according to the best separation point and then generates a Kaplan-Meier curve. Log-rank tests were used to compare overall survival between different groups. The P values were adjusted for multiple testing based on the false discovery rate (FDR) according to the Benjamini-Hochberg method.

Co-expression network construction

WGCNA is a systematic biological method used to build gene co-expression networks to mine network modules closely associated with clinical traits [10]. In the present study, we used the immunophenoscores of the infiltrated immune cell types in every sample as the target clinical traits. As read counts follow a negative binomial distribution, the RNAseq data of TCGA OV were normalized with the voom methodology of the limma R package [16]. This method estimates the mean variance of the log counts and generates a precision weight for each observation. The top 25% (6383) of genes with the highest median absolute deviation (MAD) used as a robust measure of variability were selected for WGCNA. Meanwhile, we removed 227 genes that were included in the gene sets we used.

Next, the average linkage method was performed for all pairwise genes to construct a co-expression similarity matrix. The co-expression similarity matrix was then transformed into the adjacency matrix by choosing the power of β = 4 as the soft-thresholding parameter to ensure an unsigned scale-free network. Then, we created a topological matrix using the topological overlap measure (TOM) [17]. To classify genes with similar expression patterns into gene modules, the dynamic hybrid cut method according to TOM-based dissimilarity was performed with the following major parameters: minModuleSize of 30 and deepSplit of 3. Finally, a cut-line (0.25) was selected for the module dendrogram, and some modules were merged according to the dissimilarity of estimated module eigengenes (MEs), which were defined as the first principal components of a given module and represent gene expression patterns in a module [18].

Identification of clinically significant modules and hub genes

The interesting module was identified by calculating the relevance between clinical traits and MEs, which were the first principal components of a given module. Here, we chose the immunophenoscores of the infiltrated immune cell types as the clinical traits. The module that highly correlated with the target clinical trait was selected for further analysis.

Hub genes that were defined as highly interconnected with nodes in a module have been shown to be functionally significant. Three approaches were used to identify hub genes in this study. First, potential hub genes were defined by module connectivity (Pearson’s correlation of module membership > 0.8) and clinical characteristic relations (Pearson’s correlation of gene significance > 0.2). Module membership (MM), which quantifies how close a gene is to a given module, was referred to as the correlation between the ME and the gene expression profile. Gene significance (GS) was defined as the log10 transformation of the p value of each gene in the linear regression between gene expression and the clinical traits. Second, intramodule connectivity represents the relationship between genes within a specific module. The top 30 modules of interest with a certain intramodular connectivity value were identified as candidate hub genes. Genes with a connectivity degree ≥5 (the connectivity weight threshold was set to 0.25) in the co-expression network were defined as hub genes, and the modules of interest were constructed using Cytoscape 3.8.0 [19]. The genes with high MM, GS, intramodular connectivity and connectivity degree were considered hub genes in the module of interest.

Pathway enrichment analysis

To further explore the biological significance of the hub genes, Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis was conducted on the hub genes based on the clusterProfiler package. An enriched pathway with a p-value ≤0.05 was considered to be statistically significant.

Validation of hub genes

To further explore the prognostic value of hub genes in ovarian cancer, a survival analysis of hub genes was conducted using univariate Cox regression and Kaplan-Meier plotter, as mentioned above. To determine the expression values of hub genes in normal ovarian tissues and ovarian cancer tissues, we downloaded data from the UCSC Xena website, in which TPM values of TCGA and GTEx data were extracted, log2(x + 0.001) transformed, and combined. The hub genes with different expression levels and influence on survival were considered the “real” hub genes.

CeRNA regulatory network

The lncRNA-miRNA-mRNA ceRNA network was constructed based on the ceRNA hypothesis that lncRNAs regulate the activity of mRNAs by sequestering and binding miRNAs, thereby acting as miRNA sponges. One R package called multiMiR was used to predict interactions between the mRNAs of hub genes and miRNAs. The multiMiR database [20] contains human and mouse data from 14 external databases that are categorized into three components: the three validated miRNA–target databases (miRecords [21], miRTarBase [22] and TarBase [23]), the eight predicted miRNA–target databases and the three disease−/drug-related miRNA databases. Only experimentally validated miRNAs were collected to construct the ceRNA network. StarBase [13] (http://starbase.sysu.edu.cn/) was used to predict interactions between lncRNAs and miRNAs. We chose lncRNAs whose ensemble IDs were included in the GENCODE database, those with a pancancerNum > 10 (pan-cancer types in which miRNA targets show anticorrelation relationships) and those with a clipExpNum > 0 (the number of CLIP-seq experiments) to construct the ceRNA network in Cytoscape 3.8.0 [19].

Results

Calculation of the Immunophenoscore and survival analysis

The immunophenoscores of 379 TCGA-OV samples were calculated through the GSVA package and used as clinical traits for the survival analysis and WGCNA. Because there were 374 primary tumor samples and 5 recurrent tumor samples, we extracted the immunophenoscores of the 374 patients with primary tumor for survival analysis and used all of them for WGCNA. The immunophenoscores of 28 infiltrated immune cell types from 374 primary ovarian cancer samples are shown in the boxplot in Fig. 1A.

Fig. 1
figure 1

A Boxplot of 28 immune cell immunophenoscores in 379 TCGA-OV samples. B Forest plot of the hazard ratios of each infiltrated immune cell type for OS. Three significant infiltrated immune cell types are underlined with a red dashed line. C Kaplan-Meier curves for OS times. Survival analysis according to 3 infiltrated immune cell types. Red lines represent high immunophenoscores of the real hub genes, and blue lines represent low immunophenoscores

Full size image

To determine the predictive value of immune cell infiltration, we performed univariate Cox regression analysis. We found that three infiltrated immune cell types, namely, activated B cells (HR: 0.1500; P = 0.047), CD56 bright natural killer cells (HR: 0.0082; P = 0.036) and memory B cells (HR: 34.00; P = 0.024), were significantly correlated with OS (underlined with the red dotted line in the forest plot in Fig. 1B). Activated B cells and CD56 bright natural killer cells could be seen as protective factors; however, the HR (hazard ratio) of memory B cells was greater than 1, which means that this cell type is a poor prognostic factor. Figure 1C shows the Kaplan-Meier survival curves based on the best cutoff of the immunophenoscore in these three infiltrated immune cell types.

Weighted co-expression network construction

The 379 samples with clinical information (immunophenoscores of 28 infiltrated immune cell types) were clustered using the average linkage method and Pearson’s correlation method. The expression values of 6156 genes were used to construct the co-expression gene networks, and a sample dendrogram and trait heatmap were constructed (Fig. 2A). In this study, a power of β = 4 was selected to ensure a scale-free network (Fig. 2B, C) (scale-free R2 = 0.97, slope = − 1.97).

Fig. 2
figure 2

A Sample tree and trait heatmap. The clinical trait information includes the immunophenoscores of 28 infiltrated immune cell types in each sample. The larger the immunophenoscore is, the darker the color. B Analysis of the scale-free fit index for various soft-thresholding powers (β) and the mean connectivity for the soft-thresholding powers. C Histogram of connectivity distribution when β = 4 and determining the scale-free topology when β = 4

Full size image

After merging some modules through a cut-line (0.25) (Fig. 3A), a total of 14 modules were identified by the dynamic tree cut method. The clustering dendrograms of genes are shown in Fig. 3B. A heat map for the module-trait relationship is shown in Fig. 3C.

Fig. 3
figure 3

A Clustering of module eigengenes. A cut-line (0.25) was selected for the module dendrogram, and some modules were merged according to the dissimilarity of estimated module eigengenes. B A cluster dendrogram that presents 14 gene co-expression modules was built based on the dissimilarity of the topological overlap. The gray module indicates no co-expression between the genes. C Correlated heatmap of the adjacency of modules in the WGCNA. In the correlated heatmap plot, light blue represents low adjacency, while red represents high adjacency

Full size image

Notably, the green module had the highest association with immune infiltration in the heat map of the module-trait relationship, and this module was selected as the clinically significant module for further analysis. However, only activated B cells among the three significant immune infiltration cell types had a strong correlation (0.75) with the green module. Therefore, 45 genes with high connectivity with activated B cells in the module, referred to as genes_MM, were selected based on the cutoff criteria (|MM| > 0.8 and |GS| > 0.2) (Fig. 4A). Meanwhile, we chose the top 30 genes with the highest intramodular connectivity, referred to as genes_inmodule. In addition, 23 genes were identified through Cytoscape (Fig. 4C) to have a connectivity degree ≥5 and a connectivity weight > 0.25 in the co-expression network; these genes were referred to as genes_TOM. Finally, 18 hub genes were identified as having high MM, GS, intramodular connectivity and connectivity degrees (Fig. 4B).

Fig. 4
figure 4

A Scatter plot of module eigengenes in the green module related to histologic grade. B Venn diagram. Selection of hub genes via WGCNA with 3 methods. C Visualization of the weighted gene correlation network in the green module. Red represents the real hub genes, green represents other hub genes identified through WGCNA, and light blue represents nodes in genes_TOM

Full size image

Pathway enrichment analysis

According to Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis, our results demonstrated that these hub genes are mainly involved in the following pathways: cytokine−cytokine receptor interaction, viral protein interaction with cytokine and cytokine receptor, T cell receptor signaling pathway and primary immunodeficiency (Fig. 5A). These results indicate that the hub genes in the clinically significant module are mainly involved in the regulation of the immune system.

Fig. 5
figure 5

A Bubble diagram of the KEGG pathways analysis. B Forest plot of the hazard ratios of each hub gene for OS. Seven significant real hub genes are underlined with red dashed lines. C Kaplan-Meier curves for OS times. Survival analysis according to 2 real hub genes. Red lines represent high expression of the real hub genes, and blue lines represent low expression. D The expression of 2 real hub genes between TCGA and GTEx samples is shown in a boxplot. Red plots represent gene expression in TCGA samples, and blue plots represent gene expression in GTEx samples

Full size image

Validation of hub genes

There were 7 genes among the 18 hub genes with statistical significance when conducting the survival analysis, as shown in the forest plot (Fig. 5B) (CD5, CXCL9, IL2RG, SLAMF6 and SLAMF7). All of them had P values < 0.05; however, their HRs were close to 1, which meant that these genes had little effect on overall survival (Fig. 5B). LY9 (HR = 0.71) and SLAMF1 (HR = 0.76) were selected as the real hub genes. Moreover, based on TCGA and GTEx data, the expression levels of these 2 real hub genes were significantly higher in tumor tissues than in normal tissues (Fig. 5C).

A lncRNA-miRNA-mRNA ceRNA network is constructed

During the following step, we found that 43 experimentally validated miRNAs had relationships with the 2 hub genes. The network also included 12 lncRNAs that were obtained from StarBase [13] and connected with these miRNAs. Finally, a lncRNA-miRNA-mRNA ceRNA network was constructed in Cytoscape 3.8.0 and is shown in Fig. 6.

Fig. 6
figure 6

A lncRNA-miRNA-mRNA ceRNA network of immune infiltration in ovarian cancer was constructed with 12 lncRNAs, 43 miRNAs and 2 mRNAs

Full size image

Discussion

Ovarian cancer is one of the most lethal malignant gynecologic tumors worldwide [1] and is usually diagnosed at an advanced stage. Nearly all patients suffering from ovarian cancer inevitably undergo multistep development, such as occurrence, progression, regression, recurrence and chemotherapy resistance. Platinum-based chemotherapy remains the mainstay of treatment for ovarian cancer. In patients with recurrence, the platinum-free interval (PFI) is the most important factor for PFS and OS. The longer the platinum-free interval is, the higher the response rate (RR) and the longer the duration of response to treatment after first-line therapy [24]. There is an unmet need to prolong the PFI by maintenance therapy. PARP inhibitors are emerging as a promising maintenance treatment for high-grade serous ovarian cancers (HGSOCs) with germline or somatic BRCA1/2 mutations. The median progression-free survival time was approximately 36 months longer in the olaparib group than in the placebo group in the SOLO1 clinical trial [25]. Unfortunately, only approximately 20–30% of HGSOCs have BRCA1/2 mutations [26, 27]. The FDA approved bevacizumab, a humanized anti-vascular endothelial growth factor monoclonal antibody, in combination with chemotherapy for frontline and maintenance therapy for women with newly diagnosed ovarian cancer. However, the costly treatment and the use of bevacizumab during and up to 10 months after carboplatin and paclitaxel chemotherapy prolong the median progression-free survival by only 4 months in patients with advanced epithelial ovarian cancer [28]. Thus, new treatment strategies for ovarian cancer are urgently needed. Currently, immunotherapies such as anti-PD-1/PD-L1 antibodies are attracting attention worldwide. However, the response rate to anti-PD-1 antibodies is only 20–30% in various cancer types, and the response rate in ovarian cancer is lower than that in other malignancies (11.5–15%) [29]. Moreover, immune checkpoint inhibitors have serious adverse effects [30], such as diarrhea, rash, alopecia, pneumonitis, and hepatic or gastrointestinal adverse events. We urgently need new immune checkpoint inhibitors or to understand regulatory mechanisms targeted to ovarian cancer to improve the response rate and reduce adverse effects.

In this study, we calculated the immunophenoscores of 28 infiltrated immune cell types in cystic, mucinous and serous neoplasm ovarian cancer samples from the TCGA. Through univariate Cox regression analysis, we found that activated B cells, CD56 bright natural killer cells and memory B cells significantly influenced the survival prognosis of patients with ovarian cancer. Through WGCNA, we found that immune infiltration obviously had the highest association with the green module; however, only activated B cells among the three significantly infiltrating immune cell types had a strong correlation with the green module (0.75). Therefore, we selected activated B cells for further analysis. B and T lymphocyte attenuators (BTLAs) were reported to be identified mostly on B lymphocytes rather than on T lymphocytes and natural killer cells, and the combination of chemotherapy and the anti-BTLA antibody reduced the peritoneal tumor volume and extended survival in tumor-bearing mice [31]. Anne Montfort, Oliver Pearce and Eleni Maniati reported that B cells mainly infiltrated lymphoid structures in the stroma of high-grade serous ovarian cancer metastases and that there was a strong B-cell memory response directed at a restricted repertoire of antigens and the production of tumor-specific IgGs by plasma cells [32].

To further explore the genes related to the regulation of immune infiltration, we conducted WGCNA using the immunophenoscores of 28 infiltrated immune cell types as the target clinical traits. The green module was found to have the highest association with immune infiltration. Finally, 7 genes were screened from the module. However, the HRs of CD5, CXCL9, IL2RG, SLAMF6 and SLAMF7 were very close to 1. We selected LY9 and SLAMF1 as the real hub genes. LY9 (SLAMF3) and SLAMF1 are members of the signaling lymphocyte activation molecule family (SLAMF), which is a collection of nine surface receptors expressed mainly on hematopoietic cells. SLAMF receptors are expressed on B cells in the healthy and disease states and play a pivotal role in the control of malignant cell survival, interaction with cells in their tumor microenvironment, and retention in the supporting niches [33]. As shown in Fig. 5D, based on TCGA and GTEx data, the mRNA expression levels of LY9 and SLAMF1 were significantly higher in tumor tissues than in normal ovarian tissues. Moreover, as shown in Fig. 5C, the expression of LY9 and SLAMF1 was significantly associated with the overall survival of ovarian cancer patients. Thus, the expression of LY9 and SLAMF1 was higher in ovarian cancer tissues than in normal ovarian tissues, and the high expression of these 2 genes was relevant to the prognosis of ovarian cancer. Because they are mainly expressed on hematopoietic cells, they are a major research topic in the field of leukemia. The SLAMF1 receptor was reported to be an important modulator of the BCR (B cell receptor) signaling axis and may improve immune control in chronic lymphocytic leukemia by interfering with NK cells [34]. SLAMF1-deficient cells are resistant to drugs that activate autophagy, and these results indicate that SLAMF1 expression potentially affects drug responses in chronic lymphocytic leukemia [35]. LY9 overexpression in cancerous cells could represent a potential therapeutic strategy to improve the drug sensitivity of resistant cells [36]. These results indicate that LY9 and SLAMF1 might be potential therapeutic targets of ovarian cancer.

Accumulating evidence has indicated that ceRNAs are involved in the occurrence and development of cancer. Here, we identified 12 lncRNAs, 43 miRNAs and 2 mRNAs to construct a lncRNA-miRNA-mRNA ceRNA network. MiRNAs are the center of the ceRNA network, and all of these miRNAs were experimentally validated. Hsa-miR-15b-5p was reported to have a significantly higher expression level in ovarian cancer tissue than in normal tissue [37]. LINC00511 is highly expressed in breast cancer and correlated with a poor prognosis [38]. Junjie Zhang also predicted that LINC00511 bound to hsa-miR-195-5p in a ceRNA network of hepatocellular carcinoma [39]. LINC00511 plays an oncogenic function by interacting with EZH2 and repressing P21 expression in ovarian cancer cells [40]. The lncRNA NEAT1 regulates the proliferation and apoptosis of ovarian cancer cells, providing a potential therapeutic approach for ovarian cancer [41]. LINC00662 is highly expressed in ovarian cancer tissues, and the increased expression of LINC00662 is associated with short overall survival [42]. LINC01133 was reported to repress ovarian cancer cell proliferation, invasion, migration, and tumorigenic ability [43]. MEG3 regulates the PTEN gene in ovarian cancer cells to prohibit cell proliferation, promote apoptosis and block cell cycle progression [44]. Although our study lacked in vivo and in vitro validation, the results of our study provide possible prognostic markers for immunotherapy for ovarian cancer.

Conclusions

Immune infiltrating activated B cells are associated with the survival prognosis of ovarian cancer patients. Moreover, LY9 and SLAMF1 were identified as the real hub genes associated with the overall survival of ovarian cancer patients by affecting the infiltration of activated B cells. LY9 and SLAMF1 might be potential therapeutic targets of ovarian cancer. The lncRNA-miRNA-mRNA ceRNA network demonstrates the molecular regulatory mechanism of the 2 hub genes. These findings improve our understanding of the regulatory mechanism of and provide potential therapeutic targets for immunotherapy for ovarian cancer.

Availability of data and materials

The datasets used and/or analyzed during the current study are available from the corresponding author on reasonable request.

Abbreviations

TCGA:

The Cancer Genome Atlas

UCSC:

University of California, Santa Cruz

ssGSEA:

Single-sample gene set enrichment analysis

WGCNA:

weighted gene co-expression network analysis

PARP:

Poly (ADP-ribose) polymerase

PD-L1:

programmed death ligand 1

ceRNA:

competing endogenous RNA

lncRNAs:

long noncoding RNAs

circRNAs:

circular RNAs

mRNAs:

messenger RNAs

miRNAs:

microRNAs

NES:

normalized enrichment score

FDR:

false discovery rate

MEs:

module eigengenes

GS:

Gene significance

KEGG:

Kyoto Encyclopedia of Genes and Genomes

GTEx:

Genotype-Tissue Expression

BTLAs:

B and T lymphocyte attenuators

SLAMF:

signaling lymphocyte activation molecule family

BCR:

B cell receptor

MAD:

median absolute deviation

OS:

Overall Survival

PFS:

Progression Free Survival

References

  1. Miller KD, Nogueira L, Mariotto AB, Rowland JH, Yabroff KR, Alfano CM, et al. Cancer treatment and survivorship statistics, 2019. CA Cancer J Clin. 2019;69(5):363–85. https://doi.org/10.3322/caac.21565.

    Article  PubMed  Google Scholar 

  2. Hamanishi J, Mandai M, Ikeda T, Minami M, Kawaguchi A, Murayama T, et al. Safety and antitumor activity of anti-PD-1 antibody, Nivolumab, in patients with platinum-resistant ovarian Cancer. J Clin Oncol. 2015;33(34):4015–22. https://doi.org/10.1200/JCO.2015.62.3397.

    Article  CAS  PubMed  Google Scholar 

  3. Matulonis UA, Shapira-Frommer R, Santin AD, Lisyanskaya AS, Pignata S, Vergote I, et al. Antitumor activity and safety of pembrolizumab in patients with advanced recurrent ovarian cancer: results from the phase II KEYNOTE-100 study. Ann Oncol. 2019;30(7):1080–7. https://doi.org/10.1093/annonc/mdz135.

    Article  CAS  PubMed  Google Scholar 

  4. Ma W, Gilligan BM, Yuan J, Li T. Current status and perspectives in translational biomarker research for PD-1/PD-L1 immune checkpoint blockade therapy. J Hematol Oncol. 2016;9(1):47. https://doi.org/10.1186/s13045-016-0277-y.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  5. Marcus L, Lemery SJ, Keegan P, Pazdur R. FDA approval summary: Pembrolizumab for the treatment of microsatellite instability-high solid tumors. Clin Cancer Res. 2019;25(13):3753–8. https://doi.org/10.1158/1078-0432.CCR-18-4070.

    Article  CAS  PubMed  Google Scholar 

  6. Le DT, Durham JN, Smith KN, Wang H, Bartlett BR, Aulakh LK, et al. Mismatch repair deficiency predicts response of solid tumors to PD-1 blockade. Science. 2017;357(6349):409–13. https://doi.org/10.1126/science.aan6733.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  7. Yarchoan M, Hopkins A, Jaffee EM. Tumor mutational burden and response rate to PD-1 inhibition. N Engl J Med. 2017;377(25):2500–1. https://doi.org/10.1056/NEJMc1713444.

    Article  PubMed  PubMed Central  Google Scholar 

  8. Chalmers ZR, Connelly CF, Fabrizio D, Gay L, Ali SM, Ennis R, et al. Analysis of 100,000 human cancer genomes reveals the landscape of tumor mutational burden. Genome Med. 2017;9(1):34. https://doi.org/10.1186/s13073-017-0424-2.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  9. Danaher P, Warren S, Lu R, Samayoa J, Sullivan A, Pekker I, et al. Pan-cancer adaptive immune resistance as defined by the tumor inflammation signature (TIS): results from the Cancer genome atlas (TCGA). J Immunother Cancer. 2018;6(1):63. https://doi.org/10.1186/s40425-018-0367-1.

    Article  PubMed  PubMed Central  Google Scholar 

  10. Langfelder P, Horvath S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics. 2008;9(1):559. https://doi.org/10.1186/1471-2105-9-559.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  11. Qi X, Zhang DH, Wu N, Xiao JH, Wang X, Ma W. ceRNA in cancer: possible functions and clinical implications. J Med Genet. 2015;52(10):710–8. https://doi.org/10.1136/jmedgenet-2015-103334.

    Article  CAS  PubMed  Google Scholar 

  12. Ala U, Karreth FA, Bosia C, Pagnani A, Taulli R, Leopold V, et al. Integrated transcriptional and competitive endogenous RNA networks are cross-regulated in permissive molecular environments. Proc Natl Acad Sci U S A. 2013;110(18):7154–9. https://doi.org/10.1073/pnas.1222509110.

    Article  PubMed  PubMed Central  Google Scholar 

  13. Li JH, Liu S, Zhou H, Qu LH, Yang JH. starBase v2.0: decoding miRNA-ceRNA, miRNA-ncRNA and protein-RNA interaction networks from large-scale CLIP-Seq data. Nucleic Acids Res. 2014;42(Database issue):D92–7. https://doi.org/10.1093/nar/gkt1248.

    Article  CAS  PubMed  Google Scholar 

  14. Charoentong P, Finotello F, Angelova M, Mayer C, Efremova M, Rieder D, et al. Pan-cancer Immunogenomic analyses reveal genotype-Immunophenotype relationships and predictors of response to checkpoint blockade. Cell Rep. 2017;18(1):248–62. https://doi.org/10.1016/j.celrep.2016.12.019.

    Article  CAS  PubMed  Google Scholar 

  15. Finotello F, Trajanoski Z. Quantifying tumor-infiltrating immune cells from transcriptomics data. Cancer Immunol Immunother. 2018;67(7):1031–40. https://doi.org/10.1007/s00262-018-2150-z.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  16. Law CW, Chen Y, Shi W, Smyth GK. voom: Precision weights unlock linear model analysis tools for RNA-seq read counts. Genome Biol. 2014;15(2):R29.

    Article  PubMed  PubMed Central  Google Scholar 

  17. Li A, Horvath S. Network neighborhood analysis with the multi-node topological overlap measure. Bioinformatics. 2007;23(2):222–31. https://doi.org/10.1093/bioinformatics/btl581.

    Article  CAS  PubMed  Google Scholar 

  18. Langfelder P, Horvath S. Eigengene networks for studying the relationships between co-expression modules. BMC Syst Biol. 2007;1(1):54. https://doi.org/10.1186/1752-0509-1-54.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  19. Shannon P, Markiel A, Ozier O, Baliga NS, Wang JT, Ramage D, et al. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 2003;13(11):2498–504. https://doi.org/10.1101/gr.1239303.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  20. Ru Y, Kechris KJ, Tabakoff B, Hoffman P, Radcliffe RA, Bowler R, et al. The multiMiR R package and database: integration of microRNA-target interactions along with their disease and drug associations. Nucleic Acids Res. 2014;42(17):e133. https://doi.org/10.1093/nar/gku631.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  21. Xiao F, Zuo Z, Cai G, Kang S, Gao X, Li T. miRecords: an integrated resource for microRNA-target interactions. Nucleic Acids Res. 2009;37(Database issue):D105–10. https://doi.org/10.1093/nar/gkn851.

    Article  CAS  PubMed  Google Scholar 

  22. Hsu SD, Lin FM, Wu WY, Liang C, Huang WC, Chan WL, et al. miRTarBase: a database curates experimentally validated microRNA-target interactions. Nucleic Acids Res. 2011;39(Database issue):D163–9. https://doi.org/10.1093/nar/gkq1107.

    Article  CAS  PubMed  Google Scholar 

  23. Vergoulis T, Vlachos IS, Alexiou P, Georgakilas G, Maragkakis M, Reczko M, et al. TarBase 6.0: capturing the exponential growth of miRNA targets with experimental support. Nucleic Acids Res. 2012;40(Database issue):D222–9. https://doi.org/10.1093/nar/gkr1161.

    Article  CAS  PubMed  Google Scholar 

  24. Chuang YT, Chang CL. Extending platinum-free interval in partially platinum-sensitive recurrent ovarian cancer by a non-platinum regimen: its possible clinical significance. Taiwan J Obstet Gynecol. 2012;51(3):336–41. https://doi.org/10.1016/j.tjog.2012.07.003.

    Article  PubMed  Google Scholar 

  25. Moore K, Colombo N, Scambia G, Kim BG, Oaknin A, Friedlander M, et al. Maintenance Olaparib in patients with newly diagnosed advanced ovarian Cancer. N Engl J Med. 2018;379(26):2495–505. https://doi.org/10.1056/NEJMoa1810858.

    Article  CAS  PubMed  Google Scholar 

  26. Alsop K, Fereday S, Meldrum C, DeFazio A, Emmanuel C, George J, et al. BRCA mutation frequency and patterns of treatment response in BRCA mutation-positive women with ovarian cancer: a report from the Australian ovarian Cancer study group. J Clin Oncol. 2012;30(21):2654–63. https://doi.org/10.1200/JCO.2011.39.8545.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  27. Gupta S, Nag S, Aggarwal S, Rauthan A, Warrier N. Maintenance therapy for recurrent epithelial ovarian cancer: current therapies and future perspectives - a review. J Ovarian Res. 2019;12(1):103. https://doi.org/10.1186/s13048-019-0579-0.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  28. Burger RA, Brady MF, Bookman MA, Fleming GF, Monk BJ, Huang H, et al. Incorporation of bevacizumab in the primary treatment of ovarian cancer. N Engl J Med. 2011;365(26):2473–83. https://doi.org/10.1056/NEJMoa1104390.

    Article  CAS  PubMed  Google Scholar 

  29. Callahan MK, Postow MA, Wolchok JD. Targeting T cell co-receptors for Cancer therapy. Immunity. 2016;44(5):1069–78. https://doi.org/10.1016/j.immuni.2016.04.023.

    Article  CAS  PubMed  Google Scholar 

  30. Topalian SL, Hodi FS, Brahmer JR, Gettinger SN, Smith DC, McDermott DF, et al. Safety, activity, and immune correlates of anti-PD-1 antibody in cancer. N Engl J Med. 2012;366(26):2443–54. https://doi.org/10.1056/NEJMoa1200690.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  31. Chen YL, Lin HW, Chien CL, Lai YL, Sun WZ, Chen CA, et al. BTLA blockade enhances Cancer therapy by inhibiting IL-6/IL-10-induced CD19(high) B lymphocytes. J Immunother Cancer. 2019;7(1):313. https://doi.org/10.1186/s40425-019-0744-4.

    Article  PubMed  PubMed Central  Google Scholar 

  32. Montfort A, Pearce O, Maniati E, Vincent BG, Bixby L, Bohm S, et al. A strong B-cell response is part of the immune landscape in human high-grade serous ovarian metastases. Clin Cancer Res. 2017;23(1):250–62. https://doi.org/10.1158/1078-0432.CCR-16-0081.

    Article  CAS  PubMed  Google Scholar 

  33. Shachar I, Barak A, Lewinsky H, Sever L, Radomir L. SLAMF receptors on normal and malignant B cells. Clin Immunol. 2019;204:23–30. https://doi.org/10.1016/j.clim.2018.10.020.

    Article  CAS  PubMed  Google Scholar 

  34. von Wenserski L, Schultheiß C, Bolz S, Schliffke S, Simnica D, Willscher E, et al. SLAMF receptors negatively regulate B cell receptor signaling in chronic lymphocytic leukemia via recruitment of prohibitin-2. Leukemia. 2021;35(4):1073-86.

  35. Bologna C, Buonincontri R, Serra S, Vaisitti T, Audrito V, Brusa D, et al. SLAMF1 regulation of chemotaxis and autophagy determines CLL patient response. J Clin Invest. 2016;126(1):181–94. https://doi.org/10.1172/JCI83013.

    Article  PubMed  Google Scholar 

  36. Fouquet G, Debuysscher V, Ouled-Haddou H, Eugenio MS, Demey B, Singh AR, et al. Hepatocyte SLAMF3 reduced specifically the multidrugs resistance protein MRP-1 and increases HCC cells sensitization to anti-cancer drugs. Oncotarget. 2016;7(22):32493–503. https://doi.org/10.18632/oncotarget.8679.

    Article  PubMed  PubMed Central  Google Scholar 

  37. Wang L, Zhu MJ, Ren AM, Wu HF, Han WM, Tan RY, et al. A ten-microRNA signature identified from a genome-wide microRNA expression profiling in human epithelial ovarian cancer. PLoS One. 2014;9(5):e96472. https://doi.org/10.1371/journal.pone.0096472.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  38. Lu G, Li Y, Ma Y, Lu J, Chen Y, Jiang Q, et al. Long noncoding RNA LINC00511 contributes to breast cancer tumourigenesis and stemness by inducing the miR-185-3p/E2F1/Nanog axis. J Exp Clin Cancer Res. 2018;37(1):289. https://doi.org/10.1186/s13046-018-0945-6.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  39. Zhang J, Lou W. A key mRNA-miRNA-lncRNA competing endogenous RNA triple sub-network linked to diagnosis and prognosis of hepatocellular carcinoma. Front Oncol. 2020;10:340. https://doi.org/10.3389/fonc.2020.00340.

    Article  PubMed  PubMed Central  Google Scholar 

  40. Wang J, Tian Y, Zheng H, Ding Y, Wang X. An integrated analysis reveals the oncogenic function of lncRNA LINC00511 in human ovarian cancer. Cancer Med. 2019;8(6):3026–35. https://doi.org/10.1002/cam4.2171.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  41. Ding N, Wu H, Tao T, Peng E. NEAT1 regulates cell proliferation and apoptosis of ovarian cancer by miR-34a-5p/BCL2. Onco Targets Ther. 2017;10:4905–15. https://doi.org/10.2147/OTT.S142446.

    Article  PubMed  PubMed Central  Google Scholar 

  42. Tao LM, Gong YF, Yang HM, Pei JH, Zhao XJ, Liu SS. LINC00662 promotes glycolysis and cell survival by regulating miR- 375/HIF-1alpha axis in ovarian cancer. J Biol Regul Homeost Agents. 2020;34(3):467–77. https://doi.org/10.23812/19-300-A-18.

    Article  CAS  PubMed  Google Scholar 

  43. Liu M, Shen C, Wang C. Long noncoding RNA LINC01133 confers tumor-suppressive functions in ovarian Cancer by regulating leucine-rich repeat kinase 2 as an miR-205 sponge. Am J Pathol. 2019;189(11):2323–39. https://doi.org/10.1016/j.ajpath.2019.07.020.

    Article  CAS  PubMed  Google Scholar 

  44. Wang J, Xu W, He Y, Xia Q, Liu S. LncRNA MEG3 impacts proliferation, invasion, and migration of ovarian cancer cells through regulating PTEN. Inflamm Res. 2018;67(11–12):927–36. https://doi.org/10.1007/s00011-018-1186-z.

    Article  CAS  PubMed  Google Scholar 

Download references

Acknowledgments

We thank Dr. Jianming Zeng (University of Macau) and all the members of his bioinformatics team and biotrainees, especially Xiaojie Sun, for generously sharing their experience and codes. We also thank Guozi Learn Bioinformatics (WeChat Official Accounts).

Funding

This study was supported by grants from the National Natural Science Foundation of China (grant number 81802590 to Jiangdong Xiang and grant number 82002736 to Chengjuan Jin). The funder had no role in the study design, data collection and analysis, decision to publish, or preparation of the manuscript. None of the authors received a salary from the funder.

Author information

Author notes

  1. Rongjia Su and Chengjuan Jin contributed equally to this work.

Authors and Affiliations

  1. Department of Obstetrics and Gynecology, Shanghai General Hospital, Shanghai Jiaotong University School of Medicine, 650 Xin Songjiang Road, Shanghai, 201620, China

    Rongjia Su, Chengjuan Jin, Lina Zhou, Yannan Cao, Menghua Kuang, Linxia Li & Jiangdong Xiang

Authors
  1. Rongjia Su
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Chengjuan Jin
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Lina Zhou
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Yannan Cao
    View author publications

    You can also search for this author in PubMed Google Scholar

  5. Menghua Kuang
    View author publications

    You can also search for this author in PubMed Google Scholar

  6. Linxia Li
    View author publications

    You can also search for this author in PubMed Google Scholar

  7. Jiangdong Xiang
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

Conceptualization and study design: RS and JX; Writing the manuscript draft: RS, CJ and MK; Review and editing: JX, LZ and RS; Data acquisition, analysis and Data interpretation: CJ, YC, LZ and LL; Supervision: JX. All authors read and approved the final manuscript.

Corresponding author

Correspondence to Jiangdong Xiang.

Ethics declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Su, R., Jin, C., Zhou, L. et al. Construction of a ceRNA network of hub genes affecting immune infiltration in ovarian cancer identified by WGCNA. BMC Cancer 21, 970 (2021). https://doi.org/10.1186/s12885-021-08711-w

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12885-021-08711-w

Keywords

BMC Cancer

ISSN: 1471-2407

Contact us