Research Paper Volume 12, Issue 21 pp 21874—21889
Integrated analysis of transcriptomic and metabolomic data demonstrates the significant role of pyruvate carboxylase in the progression of ovarian cancer
- 1 Department of Gynecology, Affiliated Hangzhou First People's Hospital, Zhejiang University School of Medicine, Hangzhou 310006, Zhejiang Province, China
- 2 Department of Gynecology, Affiliated Hangzhou First People's Hospital, Nanjing Medical University, Hangzhou 310006, Zhejiang Province, China
Received: April 3, 2020 Accepted: August 14, 2020 Published: November 7, 2020https://doi.org/10.18632/aging.104004
How to Cite
Copyright: © 2020 Shang et al. This is an open access article distributed under the terms of the Creative Commons Attribution License (CC BY 3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
The aim of this study was to explore prognosis-related biomarkers and underlying mechanisms during ovarian carcinoma progression and development. mRNA expression profiles and GSE49997 dataset were downloaded. Survival analyses were performed for genes with high expression levels. Expression level of candidate genes was explored in four ovarian cancer cells lines. Pyruvate carboxylase (PC) was found to be one of significantly differentially expressed gene (DEG). The role of PC knockdown was analyzed in SKOV cells using cell proliferation, flow cytometric, and Transwell migration and invasion assays. DEGs and metabolites in PC-shRNA (shPC)-treated samples vs. control groups were identified. PC was a prognosis-related gene and related to metabolic pathway. Knockdown of PC regulated cell proliferation, cell cycle progression, and migration and invasion of SKOV-3 cells. Transcriptome sequencing analyses showed STAT1 and TP53 gained higher degrees in PPI network. A total of 44 metabolites were identified. These DEGs and metabolites in PC samples were related with neuroactive ligands receptor interaction, glycine, serine and threonine metabolism, and ABC transporter pathways. PC may affect the tumor biology of ovarian cancer through the dysregulation of glycine, serine, and threonine metabolism, and ABC transporter pathways, as well as STAT1 and TP53 expression.
Ovarian cancer is 9th common malignancies in women worldwide and 8th most common cause of cancer death. According to the GLOBOCAN 2018 database, there are approximately 295,414 new cases and 184,799 deaths of ovarian cancer . There are multiple risk factors in the development of ovarian cancer, such as increased ovulation over their lifetime, hormone therapy after menopause, fertility medication, and obesity [2–4]. Epithelial ovarian carcinoma (EOC) is the most common type of ovarian cancer accounting for 90% of all ovarian cancer cases  and a major cause of death from gynecologic cancers . Ovarian cancer is usually asymptomatic in the early stages and typically diagnosed in the advanced stage, which makes treatment more challenging. Furthermore, EOC is prone to recur and migrate, which often leads to poor prognosis.
Currently, surgery and cytotoxic chemotherapy can effectively improve clinical outcomes and are indicated as primary therapies for patients with ovarian cancer; however, these treatment modalities often fail to cure the disease at terminal stages . The standard treatment for advanced ovarian cancer is aggressive cytoreductive surgery combined with platinum and taxane-based chemotherapy. Despite several advances in the treatment of ovarian cancer, the survival rate is discouragingly low in recurrent chemo-resistant ovarian cancers . Therefore, it is important to discover new biomarkers for improving the diagnosis and prognosis of patients with ovarian cancer. Metabolomics has been widely used to discover key molecular changes underlying disease pathophysiology [9, 10]. Using integrated analyses of metabolomic and transcriptomic data, researchers can greatly increase the understanding of metabolic networks and biological systems. Recently, a series of studies uncovered potential biomarkers and biological processes in several cancer types using an integrated analysis of metabolomic and transcriptomic data [11–13]. Susan et al. demonstrated that the metabolism of branched-chain amino acids (BCAAs) was altered in type 2 diabetes (T2D), caused by reduced catabolism, and proposed this alteration as a novel biomarker for T2D . An additional study showed that dysregulation of the lipolytic pathway, involving lipases, contributes to the development of pancreatic cancer, and several saturated free fatty acids (FFAs) were closely related to the proliferation of pancreatic cancer cells . However, integrated analyses of metabolomics and transcriptomics data in biomarker and pathway discovery for ovarian cancers are rare. In this study, we aimed to identify potential and reliable biomarkers with prognostic value using an integrated analysis of metabolomic and transcriptomics data in ovarian cancer. According to preliminary bioinformatic analyses, pyruvate carboxylase (PC) was found to be a prognosis-related gene significantly enriched in the metabolic pathway. The potential role of PC in ovarian cancer development was investigated by assaying cell proliferation, invasion, and metastasis in ovarian tumor cells. Related mechanisms were explored using an integrated analysis of transcriptomic and metabolomics data. Our study demonstrates a promising approach to investigate the metabolic mechanism of ovarian cancer with the aim of discovering more reliable biomarkers with prognostic value.
Identifying candidate genes significantly associated with survival of EOC
To identify the candidate genes with prognostic value in ovarian cancer, we downloaded the gene expression profiles from the TCGA and GEO databases, which was followed by screening of DEGs and text mining.
Overall, we screened 1153 and 1022 genes related with prognosis from TCGA and GSE49997 dataset, respectively. After Venn diagram analysis, 66 overlapping genes were identified as candidate genes related to EOC progression (Figure 1).
Figure 1. Prognosis-related genes screened from TCGA and GSE49997 datasets. The gene expression profiles from TCGA and GEO databases were subjected to survival analyses. A total of 1153 and 1022 genes related with prognosis were obtained from TCGA and GSE49997 datasets, respectively. After Venn diagram analysis, 66 overlapping genes were identified as candidate genes related to EOC progression.
After text mining, a total of 16 genes were found to be associated with EOC (Table 1). KEGG pathway analyses showed that the metabolic pathway (hsa01100) was the most significant pathway involved with GGPS1, ME1, DSE, NTPCR (also known as C1orf57, MGC13186), PPOX, and PC (Table 2). Among the metabolic pathway-related genes, ME1 and DSE have been reported to be associated with ovarian cancer. The KM survival curve for the remaining four genes is shown in Figure 2. High expression levels of GGPS1 and NTPCR correlated with a longer patient survival time, while high expression levels of PC and PPOX were closely related to a shorter patient-survival time.
Figure 2. Survival curve for four candidate genes closely related with the prognosis of epithelial ovarian cancer. Results show that GGPS1, NTPCR (also known as C1orf57 or MGC13186), PC, and PPOX are significantly associated with the prognosis of epithelial ovarian cancer. The survival curves based on TCGA and GSE49997 datasets are listed, separately.
Table 1. Text mining analysis for 16 genes in PubMed dataset associated with ovarian cancer.
|ME1||3||PMID: 8888953|PMID: 2305835 |PMID: 3276635|
|DSE||2||PMID: 11092984|PMID: 4028011|
Table 2. Functional enrichment analysis for candidate genes screened from TCGA-ovarian cancer and GSE49997 datasets.
|hsa04310:Wnt signaling pathway||3||PRICKLE2;TCF7L1;PPP3CA||0.0088142|
|hsa05200:Pathways in cancer||3||lpar6;TCF7L1;CSF3R||0.0420139|
|hsa05412:Arrhythmogenic right ventricular cardiomyopathy (ARVC)||2||TCF7L1;CACNA1C||0.022791|
|hsa05410:Hypertrophic cardiomyopathy (HCM)||2||ACE;CACNA1C||0.024746|
|hsa04921:Oxytocin signaling pathway||2||CACNA1C;PPP3CA||0.0439704|
|hsa04022:cGMP-PKG signaling pathway||2||CACNA1C;PPP3CA||0.0456096|
|hsa04020:Calcium signaling pathway||2||CACNA1C;PPP3CA||0.0481754|
|hsa04060:Cytokine-cytokine receptor interaction||2||CSF3R;CXCL14||0.00571639|
|hsa04010:MAPK signaling pathway||2||CACNA1C;PPP3CA||0.00610035|
|hsa04151:PI3K-Akt signaling pathway||2||lpar6;CSF3R||0.00722717|
The dysregulation of metabolic pathway-related genes in ovarian cancer cells
To evaluate the role of the metabolic pathway-related genes in ovarian cancer progression, we initially conducted real-time qPCR analyses to examine the expression levels of genes in different ovarian tumor cells. The results showed that PC was significantly upregulated in four cancer cell lines (SKOV3, CAOV-3, OV-1063, and OVCAR-3) compared to that in control IOSE80 cells (Figure 3A). Other genes were dysregulated in the different cell types. For example, GGPS1 was upregulated in the OV-1063 and OVCAR-3 cells and downregulated in CAOV-3 cells; PPOX was only upregulated in OV-1063 and OVCAR-3 cells. Considering these findings, the siRNA knockdown of PC was designed to explore its roles in ovarian cancer development. The knockdown of PC was confirmed using qRT-PCR analysis, which revealed that the expression of PC was significantly reduced in SKOV3 cells-transfected with shPC-1 and shPC-3, and shPC-2 presented off-target effects (Figure 3B). The degree of knockdown for shPC-1 and shPC-3 were 66.7% and 58%, respectively. Thus, shPC-1 transfected SKOV3 cells were used for subsequent experiments.
Figure 3. Knockdown of PC can significantly inhibit cell proliferation, cell cycle progression, and cell migration and invasion. (A) The relative expression of GGPS1, NTPCR, PPOX, PC, PRICKLE2, TCF7L1, and PPP3CA in four cancer cell lines (SKOV3, CAOV-3, OV-1063, and OVCAR-3). Gene expression levels of candidate genes (GGPS1, NTPCR, PPOX, PC, PRICKLE2, TCF7L1, and PPP3CA) were examined in different ovarian cancer cells lines using real-time qPCR analysis. (B) The relative expression of PC in SKOV-3 cell after PC knockdown using shRNAs by real-time qPCR analysis. (C) Decreased expression of PC can significantly inhibit SKOV-3 cell proliferation. (D) Cell cycle analysis of SKOV-3 cells following PC knockdown. The effect of PC on cell cycle progression was examined using flow cytometry. PC knockdown significantly inhibited cell cycle transition from G1 to S phase. (E) Cell migration and invasion of ovarian cancer cells after PC knockdown. The effects of PC on migration and invasion of SKOV3 cells were evaluated using the Transwell system. PC knockdown significantly inhibited the invasive and metastatic abilities of ovarian tumor cells.
PC significantly inhibits proliferation, cell cycle progression, and cell migration and invasion in SKOV3 cells
The CCK-8 assay was conducted to evaluate the changes in proliferation of SKOV3 cells after PC knockdown. Depletion of PC significantly decreased proliferation in the shPC group compared to that of the shNC group (Figure 3C). The effect of PC on cell cycle progression was also examined using flow cytometry. The percentage of cells in the G1 phase was increased from 49.56% in the control group to 68.91% in shPC cells. In contrast, the percentage of cells in the S phase was decreased from 32.41% in the control group to 17.31% in the shPC group (Figure 3D). These results indicate that PC-knockdown significantly inhibited cell cycle transition from the G1 to S stage.
Invasion and metastasis are major events in the progression of cancers. Therefore, we further investigated the effect of PC on ovarian cell migration and invasion using the Transwell system. Results revealed that decreased expression of PC inhibited the invasive and metastatic abilities of ovarian tumor cells (Figure 3E).
Transcriptome sequencing reveals genes and pathways regulated by PC
After filtering out genes with low expression levels, a total of 22,838 gene expression matrices were obtained. The correlation between samples based on gene expression was measured using the Pearson’s correlation coefficient. A heat map of sample correlation is shown in Figure 4A. PCA analysis was performed for validation (Figure 4B) and the results showed that the samples in the different groups could be distinguished. Using thresholds of |log2FC| > 1.585 and FDR < 0.05, 1404 DEGs were identified in the shPC group compared to the shNC group, including 586 up-regulated genes and 818 down-regulated genes. Bidirectional clustering analysis illustrates that the expression profiles of DEGs are significantly different between the shPC and shNC groups (Figure 4C). We further analyzed the molecular functions of the genes and found that DEGs in shPC samples were mainly associated with steroid hormone biosynthesis and drug metabolism (Figure 4D). GSEA showed that the down-regulation of genes in the shPC group was more associated with pathogenesis of ovarian cancer compared to that in the shNC group (Figure 4E).
Figure 4. Differentially expressed genes and pathways enrichment analysis after PC knockdown based on transcriptome sequencing. (A) Sample correlation based on differential gene expression. The correlation between samples was analyzed using Pearson’s correlation coefficient based on gene expression values. There were significant positive correlations between samples. (B) Principal component analysis results. The different colored dots represent the sample group under the condition. (C) Heatmap of differentially expressed genes between shPC and shNC samples. Two-dimension clustering analysis results were visualized using heatmaps for differentially expressed genes from PC-knockdown samples compared to the normal control group. The gene expression profiles were significantly different between groups. Red represents high expression levels while blue represents low expression levels. (D) KEGG pathway enrichment analysis for differentially expressed genes. The top ten pathway terms ranked by p-value were visualized using dot plot. The vertical axis represents KEGG pathways and the horizontal axis shows differentially expressed genes. A category with a smaller p-value represents a more significant difference. (E) Gene set enrichment analysis results. The red line refers to the highest enrichment score.
A PPI network was constructed, which consisted of 1495 edges connecting 431 nodes (Figure 5A). In the PPI network, STAT1 (degree = 54) and TP53 (degree = 46) had higher degree values than other genes, which indicated the significant role of STAT1 and TP53 in ovarian cancer progression. In the miRNA-TF-gene network, several miRNAs and TFs were identified as hub factors associated with ovarian cancer. Hsa-miR-124-3p, EZH2, and EP300 were the significant nodes in the miRNA-TF-gene network (Figure 5B).
Figure 5. PPI analysis and miRNA regulatory networks. (A) PPI networks were constructed to visualize the relationships of differentially expressed genes screened from PC-knockdown samples compared to control samples. Red dots represent up-regulated genes and green dots refer to down-regulated genes. The points with a blue border refer to ovarian cancer-related genes. (B) miRNA regulatory network. A complex regulatory network was constructed to visualize connections of miRNAs, transcription factors, and genes related to ovarian cancer. The dots in red and green represent up- and down-regulated genes, respectively. Triangular nodes refer to miRNAs and v-shaped nodes represent transcription factors.
Metabolomics sequencing reveals the dysregulated metabolic profile induced by PC knockdown
After normalization, the data were subjected to PCA and PLS-DA discriminant analyses under ESI+ and ESI- modes, separately. PCA showed that the samples in the different groups could be clustered together, and there were small offsets for the QC samples, which suggested data stability (Figure 6A). PLS-DA score plots showed obvious group separation between samples, indicating that the data were stable and reliable (Figure 6B). The cumulative R2Y and Q2Y values were both close to 1.0 (Figure 6C), which indicated that the PLS-DA models were stable, credible, and supported the differences between groups.
Figure 6. Data normalization and partial least squares discriminant analysis for metabolomics data. Two-dimensional principal component analysis (A) and three-dimensional partial least squares discriminant analysis (B) were conducted to ensure the detection stability of metabolomics data analysis. (C) Orthogonal partial least squares discriminant analysis by the X-Score model.
As for the shPC vs. shNC groups, a total of 44 metabolites met the cutoff value of FDR < 0.01, |log2FC| > 1, and VIP > 1, including 22 up-regulated and 22 down-regulated metabolites. Bidirectional clustering heat maps and volcano plots were used to visualize the results (Figure 7A and 7B). The differential metabolites in the shPC samples were mainly enriched in the metabolic pathways, ABC transporters, and nicotinate and niacinamide metabolism (Figure 7C).
Figure 7. Metabolomics data analysis for the identification of potential metabolic pathways regulated by PC. (A) Heatmap of differential metabolites in shPC vs. shNC groups in positive mode and negative mode. A bidirectional clustering analysis heat map was used to visualize metabolite levels in shPC vs. shNC samples under positive (left) and negative modes (right). The horizontal and vertical axes represent the samples and metabolites, respectively. Green represents down-regulated levels while red refers to up-regulated levels. (B) Volcano map of differential metabolites under positive (left) and negative (right) modes in shPC vs. shNC samples. The points in pink, blue, and grey refer to metabolites with up, down, and normal regulated levels, respectively. The dot size represents the VIP value. (C) Functional enrichment analysis for differential metabolites. The vertical and horizontal axes represent pathway categories and count number, respectively. The dot size represents the metabolite ratio of pathway enrichment. The color changes from blue to red refer to decreasing p-values. A dot with a smaller p-value represents a more significant difference for the pathway category. (D) Integrated analysis of transcriptomic and metabolomics data to identify crucial pathways regulated by PC.
Dysregulated pathways at both transcriptomics and metabolomics levels
Using an integrated analysis of transcriptomic and metabolomics data, eight pathways were found to be dysregulated by DEGs and differential metabolites in SKOV3 cells induced by PC knockdown, such as neuroactive ligands receptor interaction; glycine, serine, and threonine metabolism; and ABC transporter pathways (Figure 7D).
Ovarian cancer is one of the most common cancers in women and is a major cause of cancer-related death worldwide . Although patients with ovarian cancer benefit from first-line treatments during early stages of cancer, most cases are diagnosed with a high-stage carcinoma. Despite advances in treatment management for ovarian cancer patients, the prognosis remains poor. The screening of potential biomarkers and a better understanding of the pathogenesis of ovarian cancer could contribute to the development of novel target therapies . Thus, in this study, we attempted to discover prognosis-related biomarkers and explore the related mechanisms underlying the development and progression of ovarian cancer.
In our study, 66 prognosis-related genes in ovarian cancer samples were identified by combining TCGA and GSE49997 datasets. A metabolic pathway was the most significant pathway enriched by GGPS1, ME1, DSE, NTPCR, PPOX, and PC.ME1 and DSE have been reported to be associated with the development of ovarian cancer in previous studies [18, 19]. Among the remaining metabolic pathway-related genes, the high expression of PC was closely related with poor survival of patients with ovarian cancer and PC was prominently overexpressed in the four ovarian cancer cell lines, making it a candidate gene for further analysis. To our knowledge, this is the first study exploring the clinical significance of PC in ovarian cancer.
Pyruvate carboxylase is an anaplerotic enzyme that catalyzes the irreversible carboxylation of pyruvate to oxaloacetate. It plays an essential role in various metabolic processes in mammals, such as gluconeogenesis, lipogenesis, biosynthesis of neurotransmitters, and glucose-induced insulin secretion . A previous study shows that PC is significantly increased in gastric adenocarcinomas and correlates with tumor stage and pathological grade . The expression and activity of PC is significantly enhanced in non-small cell lung cancer tissues and required for tumor cell proliferation . A recent study by Phannasil et al. shows that the up-regulation of PC could promote cell growth and invasion of breast cancer cells . In addition, the increased expression of PC is found in isocitrate dehydrogenase (IDH)-mutant glioma tissues and PC has been determined to be the therapeutic target . Previous reports are consistent with findings of our in vivo experiments. Results of the present study show that cell proliferation, cell cycle progression, and cell invasion and migration were significantly inhibited in PC-depleted SKOV3 cells, which indicate, like its roles in other cancer types, that PC might function as a major tumor gene to promote the progression of ovarian cancer.
In order to further explore the changes in genes and metabolites in ovarian cancer cells induced by PC knockdown, transcriptome sequencing and metabolomics sequencing were performed. Our results showed that numerous genes and metabolites were differentially expressed in PC-knockdown SKOV3 cells.
The PPI network of DEGs induced by PC knockdown suggested that STAT1 (degree = 54) and TP53 (degree = 46) were the hub genes with high degrees, which revealed the significant role of these two proteins in the alteration of gene expression. STAT1 is a member of the STAT family and functions as a signal messenger involved in many biological processes, including cell proliferation, differentiation, and apoptosis. The activation of STAT1 is triggered by JAK ligands and increased expression levels of STAT1 have been reported in patients with ovarian cancer . Overexpression or knockdown of STAT1 can directly induce or suppress ovarian cancer cell proliferation, migration, and invasion. One of the potential mechanisms might be related to crosstalk with TGF-β signaling pathways, which is a major factor in EOC progression . In addition, TP53, also known as p53, is well-established as an important tumor suppressor. Mutations in p53 have been frequently found in human tumors including ovarian cancer. Novel functions of p53 in metabolic regulation have also been reported in recent studies [27, 28]. p53 is involved in the regulation of many processes such as glycolysis, mitochondrial metabolism, and fatty acid synthesis . The analysis of the PPI network in the present study also demonstrated that STAT1 and p53 are candidate genes in ovarian cancer development, which is consistent with results from previous studies.
Additionally, metabolites are the final products in biological processes and can be influenced by genetic or environmental factors. Metabolomics analyses have demonstrated that pyruvate carboxylation acts as a key process, providing carbon backbones for downstream metabolites for the biosynthesis of cellular components such as membrane lipids, nucleotides, and amino acids . By integrating the analysis of transcriptomic and metabolomics data, we demonstrated that these DEGs and differential metabolites in PC-knockdown tumor cells were closely associated with neuroactive ligands receptor interactions; glycine, serine, and threonine metabolism; and ABC transporter pathways. As described previously in the results, p53 plays an important role in the metabolic mechanisms of tumors. Mutations in p53 can enhance the glycolytic rate of fibroblasts and disturb biosynthetic processes of serine and glycine . Furthermore, various ABC transporters were increased after EOC chemotherapy, which decreases the accumulation of chemotherapy drugs and finally induce a poor prognosis of ovarian cancer by enhancing chemo-resistance . Considering this, we concluded that the disturbance of PC could significantly influence gene expression and metabolomics in ovarian cancer cells and lead to the dysregulation of serine, glycine biosynthesis, and ABC transporter pathways in ovarian cancer.
There are some limitations in our study. First, because of heterogeneity in the samples and technology limitations, the metabolites analyzed may be incomplete and not available for all ovarian cancer cases. Second, the sample size for the survival analysis was relatively small, and more samples should be included. Lastly, the function of PC was explored at the cellular level and functional validation in clinical samples is urgently needed.
In conclusion, PC plays a significant role in the pathogenesis of ovarian cancer. The reduction of PC expression can significantly inhibit cell proliferation, cell cycle progression, and cell migration and invasion. Significantly decreased PC expression may impact ovarian tumor biology through the dysregulation of STAT1 and TP53 expression, serine and glycine biosynthesis, and ABC transporter pathways. The identified pathways in this study might provide new diagnostic and treatment options for patients with ovarian cancer.
Materials and Methods
The mRNA expression profiles associated with ovarian cancer were downloaded from The Cancer Genome Atlas (TCGA) database, which were derived from 299 tissue samples. The microarray dataset GSE49997 was downloaded from the Gene Expression Omnibus (GEO) repository (https://www.ncbi.nlm.nih.gov/geo/) based on the platform of GPL2986 ABI Human Genome Survey Microarray Version 2. The GSE49997 dataset includes 204 EOC tissue samples, of which 194 samples contain corresponding clinical information.
Identification of survival-related genes
To identify crucial genes in ovarian cancer development, the gene expression profiles of samples with clinical information from TCGA and GSE49997 datasets were further analyzed. The raw count data were downloaded and normalized using the count per million (CPM) method. Genes with low expression levels were removed, and the remaining 15,683 genes with high expression levels were subjected to survival analyses. Prognostic risk assessment was based on survival information and normalized gene expression profiles. Patients were divided into high-expression and low-expression groups based on the median value of gene expression. Kaplan-Meier (KM) curves were generated using the Survival package in R and statistically tested using the log-rank test. Genes with a p-value of less than 0.05 were considered as potential genes associated with prognosis in EOC patients. Similarly, after normalization, the gene expression datasets from EOC tissues in the GEO dataset were subjected to survival analyses. The genes closely related with prognosis based on the two databases were analyzed using Venn diagrams. The overlapping genes were then selected for further analysis.
Pathway enrichment analysis and literature search
Text mining was performed for the overlapping genes using Perl code. The published genes that are closely related with ovarian cancer were searched in the PubMed database. In addition, the overlapping genes were subjected to the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis using the Database for Annotation, Visualization and Integrated Discovery (DAVID) online tool. Pathways with p < 0.05 and counts ≥ 2 were considered significant.
Four human ovarian cancer cell lines (SKOV3, Caov-3, OV-1063, OVCAR-3), normal epithelial ovarian cells (IOSE80), and human embryonic kidney (293T) cells were purchased from the cell bank of China Academic of Science. The SKOV3 cells were cultured in McCoy’s 5A Media (modified with tricine) supplemented with 10% fetal bovine serum (FBS). The Caov-3 and 293T cell lines were maintained in 90% DMEM with 10% FBS. The OV-1063 and IOSE80 cell lines were maintained in 90% RPMI 1640 with 10% FBS, and the OVCAR-3 cells were cultured in 80% RPMI 1640 with 20% FBS, sodium pyruvate, and 0.01 mg/ml bovine insulin. All the cell lines were cultured in an atmosphere of 5% CO2 and 95% air at 37 °C.
Real-time qPCR analysis
Total RNA was extracted from cells using the RNAiso Plus (Trizol) reagent (TaKaRa, Japan) and cDNA was synthesized using the PrimeScript™RT Master Mix (Perfect Real Time) kit (TAKARA, Japan) according to the manufacturer’s instructions. Real-time PCR was performed to evaluate the expression levels of GGPS1, NTPCR, PPOX, PC, PRICKLE2, TCF7L1, and PPP3CA in tumor cells. A total of 8 μl of cDNA was used as template in a final 20 μl PCR volume containing 1 μl forward primer, 1 μl reverse primer, and 10 μl SYBR Premix EX Taq (2x). PCRs were run as follows: 50.0 °C for 3 min, 95.0 °C for 3 min, followed by 40 cycles of 95.0 °C for 10 s and 60.0 °C for 30 s. Following PCR, a melting curve was obtained at temperatures from 60 °C to 95 °C, at increments of 0.5 °C for 10 s. Primer sequences are listed in Table 3.
Table 3. The primer sequences in PCR analysis.
Based on designing sigma software (https://www.sigmaaldrich.com/life-science/functional-genomics-and-rnai/sirna/mission-predesigned-sirna.html), the shRNA sequence of PC was obtained. Three shRNAs sequences in CDS region for PC, GCCCAGTTTATGGTGCAGAAT (shPC-1), GCCAAGGAGAACAACGTAGAT (shPC-2), and ATGGGCATCCGCCTGGATAAT (shPC-3) were selected to knock down endogenous PC, and messy sequence (NC- GTTCTCCGAACGTGTCACGTC) was used as control sequence. Then, the primers were designed and synthesized by Suzhou Jinweizhi Biological Technology Co. LTD. (Table 3). Inducible shRNA construction was performed by ligating annealed oligonucleotides into the pLKO.1 Puro vector digested with EcoRI and AgeI. The ligated products were transformed into chemically competent Stbl3 E. coli using Lipofectamine 3000 reagent (Invitrogen) according to the manufacturer’s instructions and grown in LB supplemented with ampicillin. After validation of shRNA insertion using gene sequencing, the plasmids were packaged into lentiviral particles using the third-generation lentiviral production system. SKOV3 cells were divided into three groups: black (without any treatment), shNC (transfected with negative control plasmid), shPC (transfected with shRNA-PC). Stable SKOV3 cells were selected with 0.5 μg/ml of puromycin over 3 days. Lastly, the efficiency of PC knockdown in SKOV3 cells was confirmed using real-time qPCR analysis, and the degree of knockdown [(1- the relative mean expression of PC in shPC group/shNC group) × 100%] for three PC shRNAs were calculated. The sable shRNA-PC transfected SKOV3 cells with highest knockdown degree were used for following experiments.
Cell viability and cell cycle analyses
Cell viability assays were performed using the cell counting kit-8 (CCK-8; Beyotime Biotechnology, China). Briefly, the SKOV3 cells in three groups (3 × 104 cells/ml) were seeded in 96-well plates at a density of 3,000 cells/well, respectively. Cells were incubated at 37 °C in 5% CO2 atmosphere conditions for 24, 48, 72, and 96 hours, and then 10 μl of CCK-8 solution was added into each well for 1.5 hours. The OD value of each well sample was measured at 450 nm using the Multiscan Spectrum System (BD Biosciences, USA). All assays were repeated in triplicate.
Flow cytometry was used for cell cycle analyses. shPC and control SKOV3 cells were digested with trypsin and collected using centrifugation at 300 × g, 4 °C for 5 min. Then, cells were fixed in 70% ethanol overnight at 4 °C. Subsequently, the samples were treated with RNase A for 30 min and stained with propidium iodide (BD Biosciences, USA) for 15 min. Flow cytometry was immediately conducted using a FACSCAN flow cytometer (BD Biosciences, USA) following the manufacturer’s instructions. All assays were repeated in triplicate.
Transwell migration and invasion assays
Invasion and migration assays were performed using BD Matrigel culture inserts. First, 8 μm pore 24-well Transwell inserts (Corning, USA) were coated with 20 μl of Matrigel (BD Biosciences) and incubated for 30 min at 37 °C in a 5% CO2 incubator to allow gel formation. Control shNC and shPC SKOV3 cells were suspended at equal cell densities in serum-free medium. A total of 100 μl of cells (2 × 105 cells/ml) were seeded in the upper chamber of untreated Transwells for cell migration analyses and in the upper chamber of Transwells treated with Matrigel for cell invasion analyses. The lower chambers were supplemented with 500 μl of 10% FBS-containing medium. Then, cells were incubated at 37 °C in 5% CO2 atmosphere conditions for 12 and 24 h, respectively. The cells in Transwells were treated with 4% paraformaldehyde for 30 min and stained with crystal violet for 10 min. Subsequently, the upper chambers were washed three times in PBS and non-invading cells on the inner surface were carefully removed using cotton swabs. Finally, the invading cells on the reverse side of the upper chamber were counted using a microscope at three different fields to assess cell migration and invasion abilities. All assays were repeated in triplicate.
Metabolomics data acquisition and identification
Metabolomics data acquisition was performed using an ultra-performance liquid chromatography-quadrupole-time of flight liquid chromatography/mass spectrometry (UPLC-Q-TOF LC/MS) system . The samples included 10 shPC-infected SKOV3 cell samples (1 × 107 cells for each sample), 10 shNC SKOV3 cell samples, and 5 quality control (QC) samples. Cells in each sample were mixed with 500 μl of methyl alcohol/acetonitrile/distilled water (2:2:1, v/v/v) solution and disrupted using ultrasonication for 30 min. After incubation at -20 °C for 1 h, the supernatant was collected using centrifugation at 13,000 rpm at 4 °C for 15 min followed by freeze-drying. For metabolomic analyses, samples were re-dissolved in 100 μl of acetonitrile solution (1:1 ratio of acetonitrile and water, by volume) and centrifuged at 14,000 × g at 4 °C for 15 min. The supernatant was collected for LC/MS analysis.
The metabolic products acquisition and identification were both performed in positive-ion (ESI+) and negative-ion (ESI-) modes and analyzed using a Triple-TOF 5600 mass spectrometer. Data normalizations were performed on LC/MS data based on internal standards using ESI+ and ESI- modes.
Principal component analysis (PCA), extensively used in the statistical learning field, and the partial least squares discriminant analysis (PLS-DA) are commonly applied to evaluate differences between groups [33, 34]. In the present study, we used two-dimensional PCA to ensure detection stability and three-dimensional PLS-DA for cross validation. In addition, univariate nonparametric analyses and multivariate PLS-DA were performed for metabolic profiling using the ropls software, version 1.6.2 (http://bioconductor.org/packages/release/bioc/html/ropls.html). The thresholds for potential biomarker selection were set as follows: a false discovery rate (FDR) of <0.01, |log2 fold change (log2FC)| > 1 in univariate analysis, and Variable Importance in the Projection (VIP) > 1 in multivariate analysis. Bidirectional hierarchical clustering analyses were performed to assess data classification ability and concentration levels of the screened metabolites. For differential metabolites obtained from tumor cells, we transformed these data into the KEGG ID format using the MetaboAnalyst online tool  and performed KEGG pathway analyses using MBROLE 2.0 (http://csbg.cnb.csic.es/mbrole2/analysis.php) with a p-value < 0.05.
Total mRNA extractions and cDNA library preparations were performed for shPC and control shNC cell samples. Then, two-paired end sequencing was performed using the Illumina platform (Illumina, San Diego, CA, USA).
Clean reads data were obtained using the Trimmomatic tool (version 3.6) . Then, the reads were mapped to the human reference genome (GRCH38, Gencode)  with the Hisat 2 software, version 2.05 . Gene expression levels were evaluated by counting reads mapped to protein-coding regions using FeatureCounts tools (v1.6.0) , and expression values were normalized using Fragments Per Kilobase of Exon Per Million Fragments Mapped (FPKM) method. Genes with an FPKM value ≥ 0.1 in at least three samples were further analyzed. According to filtered gene abundance expression profiles, the correlation of gene expression levels between samples was analyzed using PCA with the ggord package (https://zenodo.org/badge/latestdoi/35334615) .
DEGs between the shPC and shNC groups were identified using the quasi-likelihood F-tests method of the edgeR software . Genes with |log2FC| > 1.585 and FDR < 0.05 were considered to be significant, and two-dimensional clustering heatmaps were used to visualize the gene expression profiles of DEGs. Subsequently, gene ontology (GO) function and KEGG pathway enrichments were performed for genes of interest using the Clusterprofiler package . A p-value < 0.05 was considered as a significant difference.
Ovarian carcinoma-related genes were retrieved from the DisGeNET database  (http://www.disgenet.org/web/DisGeNET/menu/home). Differentially expressed genes were subjected to gene set enrichment analysis (GSEA) using the DOSE package .
The PPI pairs with required confidence (combined score > 0.7) were obtained with the STRING online tool  (https://string-db.org/). Cytoscape software  was utilized to construct a PPI network and the network topological properties were analyzed based on degree , betweenness , and closeness  using the cytoscape CytoNCA plugin .
Moreover, the significant modules in the PPI network were mined using the MCODE tool  with a screening score > 10. For the intriguing DEGs, we performed microRNA and transcription factor (TF) predictions by using Enrichr tools . The corresponding miRNA-gene pairs and TF-gene pairs were also identified from the miRTarBase and ENCODE databases . Finally, we integrated these miRNAs, TFs, and candidate genes to construct a systemic regulatory network.
Integrated pathway analysis of transcriptomics and metabolomics data
Integrated analyses are conducted in order to understand the biological function of post-genomic data at a higher level than individual biomolecules. The “IMPaLA” web tool has been used to integrate more than one type of omics data for pathway analysis . In this study, the joint pathway analysis was performed for DEGs and differential metabolites based on the KEGG database. Lastly, pathways with number_of_overlapping_metabolites/gene > 0 and metabolite with p < 0.05 were considered to be significant.
TCGA: The Cancer Genome Atlas; GEO: Gene Expression Omnibus; PC: Pyruvate carboxylase; DEGs: Differentially expressed genes; shPC: PC-shRNA; EOC: Epithelial ovarian carcinomas; PPI: Protein-protein interaction; BCAAs: Branched-chain amino acids; T2D: Type 2 diabetes; FFAs: Free fatty acids; CPM: Count per million; KM: Kaplan-Meier; KEGG: Kyoto Encyclopedia of Genes and Genomes; DAVID: Database for Annotation, Visualization and Integrated Discovery; FBS: Fetal bovine serum; QC: Quality control; PCA: Principal component analysis; PLS-DA: Partial least squares discriminant analysis; FDR: False discovery rate; FPKM: Fragments Per Kilobase of Exon Per Million Fragments Mapped; GO: Gene ontology; GSEA: Gene set enrichment analysis; TF: Transcription factor; IDH: Isocitrate dehydrogenase.
Conception and design of the research: JT; acquisition of data: HS, JZ; analysis and interpretation of data: HS, JT; obtaining funding: JT; drafting the manuscript: HS; revision of manuscript for important intellectual content: JT; conducting experiments: JZ. All authors read and approved the final manuscript.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
This work was supported by Medical Scientific Research Foundation of Zhejiang Province, China (grant number 2019KY495).
- 1. Bray F, Ferlay J, Soerjomataram I, Siegel RL, Torre LA, Jemal A. Global cancer statistics 2018: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA Cancer J Clin. 2018; 68:394–424. https://doi.org/10.3322/caac.21492 [PubMed]
- 2. Matulonis UA, Sood AK, Fallowfield L, Howitt BE, Sehouli J, Karlan BY. Ovarian cancer. Nat Rev Dis Primers. 2016; 2:16061. https://doi.org/10.1038/nrdp.2016.61 [PubMed]
- 3. Kashyap S, Davis OK. Ovarian cancer and fertility medications: a critical appraisal. Semin Reprod Med. 2003; 21:65–71. https://doi.org/10.1055/s-2003-39996 [PubMed]
- 4. Olsen CM, Green AC, Whiteman DC, Sadeghi S, Kolahdooz F, Webb PM. Obesity and the risk of epithelial ovarian cancer: a systematic review and meta-analysis. Eur J Cancer. 2007; 43:690–709. https://doi.org/10.1016/j.ejca.2006.11.010 [PubMed]
- 5. Desai A, Xu J, Aysola K, Qin Y, Okoli C, Hariprasad R, Chinemerem U, Gates C, Reddy A, Danner O, Franklin G, Ngozi A, Cantuaria G, et al. Epithelial ovarian cancer: an overview. World J Transl Med. 2014; 3:1–8. https://doi.org/10.5528/wjtm.v3.i1.1 [PubMed]
- 6. Gurler Main H, Xie J, Muralidhar GG, Elfituri O, Xu H, Kajdacsy-Balla AA, Barbolina MV. Emergent role of the fractalkine axis in dissemination of peritoneal metastasis from epithelial ovarian carcinoma. Oncogene. 2017; 36:3025–36. https://doi.org/10.1038/onc.2016.456 [PubMed]
- 7. Berkenblit A, Cannistra SA. Advances in the management of epithelial ovarian cancer. J Reprod Med. 2005; 50:426–38. [PubMed]
Emmings E, Mullany S, Chang Z, Landen CN
Jr, Linder S, Bazzaro M. Targeting mitochondria for treatment of chemoresistant ovarian cancer. Int J Mol Sci. 2019; 20:229. https://doi.org/10.3390/ijms20010229 [PubMed]
- 9. Yin P, Xu G. Metabolomics for tumor marker discovery and identification based on chromatography-mass spectrometry. Expert Rev Mol Diagn. 2013; 13:339–48. https://doi.org/10.1586/erm.13.23 [PubMed]
- 10. Xia J, Broadhurst DI, Wilson M, Wishart DS. Translational biomarker discovery in clinical metabolomics: an introductory tutorial. Metabolomics. 2013; 9:280–99. https://doi.org/10.1007/s11306-012-0482-9 [PubMed]
- 11. Ren S, Shao Y, Zhao X, Hong CS, Wang F, Lu X, Li J, Ye G, Yan M, Zhuang Z, Xu C, Xu G, Sun Y. Integration of metabolomics and transcriptomics reveals major metabolic pathways and potential biomarker involved in prostate cancer. Mol Cell Proteomics. 2016; 15:154–63. https://doi.org/10.1074/mcp.M115.052381 [PubMed]
- 12. Auslander N, Yizhak K, Weinstock A, Budhu A, Tang W, Wang XW, Ambs S, Ruppin E. A joint analysis of transcriptomic and metabolomic data uncovers enhanced enzyme-metabolite coupling in breast cancer. Sci Rep. 2016; 6:29662. https://doi.org/10.1038/srep29662 [PubMed]
- 13. Fan T, Rong Z, Dong J, Li J, Wang K, Wang X, Li H, Chen J, Wang F, Wang J, Wang A. Metabolomic and transcriptomic profiling of hepatocellular carcinomas in Hras12V transgenic mice. Cancer Med. 2017; 6:2370–84. https://doi.org/10.1002/cam4.1177 [PubMed]
- 14. Connor SC, Hansen MK, Corner A, Smith RF, Ryan TE. Integration of metabolomics and transcriptomics data to aid biomarker discovery in type 2 diabetes. Mol Biosyst. 2010; 6:909–21. https://doi.org/10.1039/b914182k [PubMed]
- 15. Zhang G, He P, Tan H, Budhu A, Gaedcke J, Ghadimi BM, Ried T, Yfantis HG, Lee DH, Maitra A, Hanna N, Alexander HR, Hussain SP. Integration of metabolomics and transcriptomics revealed a fatty acid network exerting growth inhibitory effects in human pancreatic cancer. Clin Cancer Res. 2013; 19:4983–93. https://doi.org/10.1158/1078-0432.CCR-13-0209 [PubMed]
- 16. Wild CP, Stewart BW, Wild C. (2014). World cancer report 2014: World Health Organization, Geneva, Switzerland.
- 17. Hennessy BT, Coleman RL, Markman M. Ovarian cancer. Lancet. 2009; 374:1371–82. https://doi.org/10.1016/S0140-6736(09)61338-6 [PubMed]
- 18. Radford H, Wilson AP. A comparison of immunohistochemical staining of human cultured mesothelial cells and ovarian tumour cells using epithelial and mesothelial cell markers. Anal Cell Pathol. 1996; 11:173–82. [PubMed]
- 19. Tanaka S, Tsuda N, Kawano K, Sakamoto M, Nishida T, Hashimoto T, Shichijo S, Kamura T, Itoh K. Expression of tumor-rejection antigens in gynecologic cancers. Jpn J Cancer Res. 2000; 91:1177–84. https://doi.org/10.1111/j.1349-7006.2000.tb00902.x [PubMed]
- 20. Lao-On U, Attwood PV, Jitrapakdee S. Roles of pyruvate carboxylase in human diseases: from diabetes to cancers and infection. J Mol Med (Berl). 2018; 96:237–47. https://doi.org/10.1007/s00109-018-1622-0 [PubMed]
- 21. Kang J, Tang S, Zhao G, Xu G, Wang B, Liang S, Ding J. [Expression and clinical significance of pyruvate carboxylase in gastric adenocarcinoma]. Xi Bao Yu Fen Zi Mian Yi Xue Za Zhi. 2013; 29:85–88. [PubMed]
Sellers K, Fox MP, Bousamra M
2nd, Slone SP, Higashi RM, Miller DM, Wang Y, Yan J, Yuneva MO, Deshpande R, Lane AN, Fan TW. Pyruvate carboxylase is critical for non-small-cell lung cancer proliferation. J Clin Invest. 2015; 125:687–98. https://doi.org/10.1172/JCI72873 [PubMed]
- 23. Phannasil P, Thuwajit C, Warnnissorn M, Wallace JC, MacDonald MJ, Jitrapakdee S. Pyruvate carboxylase is up-regulated in breast cancer and essential to support growth and invasion of MDA-MB-231 cells. PLoS One. 2015; 10:e0129848. https://doi.org/10.1371/journal.pone.0129848 [PubMed]
- 24. Izquierdo-Garcia JL, Cai LM, Chaumeil MM, Eriksson P, Robinson AE, Pieper RO, Phillips JJ, Ronen SM. Glioma cells with the IDH1 mutation modulate metabolic fractional flux through pyruvate carboxylase. PLoS One. 2014; 9:e108289. https://doi.org/10.1371/journal.pone.0108289 [PubMed]
- 25. Josahkian JA, Saggioro FP, Vidotto T, Ventura HT, Candido Dos Reis FJ, de Sousa CB, Tiezzi DG, de Andrade JM, Koti M, Squire JA. Increased STAT1 expression in high grade serous ovarian cancer is associated with a better outcome. Int J Gynecol Cancer. 2018; 28:459–65. [PubMed]
- 26. Tian X, Guan W, Zhang L, Sun W, Zhou D, Lin Q, Ren W, Nadeem L, Xu G. Physical interaction of STAT1 isoforms with TGF-β receptors leads to functional crosstalk between two signaling pathways in epithelial ovarian cancer. J Exp Clin Cancer Res. 2018; 37:103. https://doi.org/10.1186/s13046-018-0773-8 [PubMed]
- 27. Flöter J, Kaymak I, Schulze A. Regulation of metabolic activity by p53. Metabolites. 2017; 7:21. https://doi.org/10.3390/metabo7020021 [PubMed]
- 28. Li T, Kon N, Jiang L, Tan M, Ludwig T, Zhao Y, Baer R, Gu W. Tumor suppression in the absence of p53-mediated cell-cycle arrest, apoptosis, and senescence. Cell. 2012; 149:1269–83. https://doi.org/10.1016/j.cell.2012.04.026 [PubMed]
- 29. Liu J, Zhang C, Hu W, Feng Z. Tumor suppressor p53 and its mutants in cancer metabolism. Cancer Lett. 2015; 356:197–203. https://doi.org/10.1016/j.canlet.2013.12.025 [PubMed]
- 30. Altman BJ, Stine ZE, Dang CV. From krebs to clinic: glutamine metabolism to cancer therapy. Nat Rev Cancer. 2016; 16:619–34. https://doi.org/10.1038/nrc.2016.71 [PubMed]
- 31. Ween MP, Armstrong MA, Oehler MK, Ricciardelli C. The role of ABC transporters in ovarian cancer progression and chemoresistance. Crit Rev Oncol Hematol. 2015; 96:220–56. https://doi.org/10.1016/j.critrevonc.2015.05.012 [PubMed]
- 32. Zoerner AA, Schroeder C, Kayacelebi AA, Suchy MT, Gutzki FM, Stichtenoth DO, Tank J, Jordan J, Tsikas D. A validated, rapid UPLC-MS/MS method for simultaneous ivabradine, reboxetine, and metoprolol analysis in human plasma and its application to clinical trial samples. J Chromatogr B Analyt Technol Biomed Life Sci. 2013; 927:105–11. https://doi.org/10.1016/j.jchromb.2013.01.016 [PubMed]
- 33. Trygg J, Holmes E, Lundstedt T. Chemometrics in metabonomics. J Proteome Res. 2007; 6:469–79. https://doi.org/10.1021/pr060594q [PubMed]
- 34. Westerhuis JA. Assessment of PLSDA cross validation. Metabolomics. 2008; 4:81–89. https://doi.org/10.1007/s11306-007-0099-6
- 35. Chong J, Soufan O, Li C, Caraus I, Li S, Bourque G, Wishart DS, Xia J. MetaboAnalyst 4.0: towards more transparent and integrative metabolomics analysis. Nucleic Acids Res. 2018; 46:W486–94. https://doi.org/10.1093/nar/gky310 [PubMed]
- 36. Bolger AM, Lohse M, Usadel B. Trimmomatic: a flexible trimmer for illumina sequence data. Bioinformatics. 2014; 30:2114–20. https://doi.org/10.1093/bioinformatics/btu170 [PubMed]
- 37. Harrow J, Frankish A, Gonzalez JM, Tapanari E, Diekhans M, Kokocinski F, Aken BL, Barrell D, Zadissa A, Searle S, Barnes I, Bignell A, Boychenko V, et al. GENCODE: the reference human genome annotation for the ENCODE project. Genome Res. 2012; 22:1760–74. https://doi.org/10.1101/gr.135350.111 [PubMed]
- 38. Nakamura Y, de Paiva Alves E, Veenstra GJ, Hoppler S. Tissue- and stage-specific Wnt target gene expression is controlled subsequent to β-catenin recruitment to cis-regulatory modules. Development. 2016; 143:1914–25. https://doi.org/10.1242/dev.131664 [PubMed]
- 39. Liao Y, Smyth GK, Shi W. featureCounts: an efficient general purpose program for assigning sequence reads to genomic features. Bioinformatics. 2014; 30:923–30. https://doi.org/10.1093/bioinformatics/btt656 [PubMed]
- 40. Beck M. ggord: Ordination Plots with ggplot2. R package version 011. 2016; 9000. (https://github.com/fawda123/ggord.
- 41. Lun AT, Chen Y, Smyth GK. It’s DE-licious: a recipe for differential expression analyses of RNA-seq experiments using quasi-likelihood methods in edgeR. Methods Mol Biol. 2016; 1418:391–416. https://doi.org/10.1007/978-1-4939-3578-9_19 [PubMed]
- 42. Yu G, Wang LG, Han Y, He QY. clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS. 2012; 16:284–87. https://doi.org/10.1089/omi.2011.0118 [PubMed]
- 43. Piñero J, Bravo À, Queralt-Rosinach N, Gutiérrez-Sacristán A, Deu-Pons J, Centeno E, García-García J, Sanz F, Furlong LI. DisGeNET: a comprehensive platform integrating information on human disease-associated genes and variants. Nucleic Acids Res. 2017; 45:D833–39. https://doi.org/10.1093/nar/gkw943 [PubMed]
- 44. Yu G, Wang LG, Yan GR, He QY. DOSE: an r/bioconductor package for disease ontology semantic and enrichment analysis. Bioinformatics. 2015; 31:608–09. https://doi.org/10.1093/bioinformatics/btu684 [PubMed]
- 45. Szklarczyk D, Franceschini A, Wyder S, Forslund K, Heller D, Huerta-Cepas J, Simonovic M, Roth A, Santos A, Tsafou KP, Kuhn M, Bork P, Jensen LJ, von Mering C. STRING v10: protein-protein interaction networks, integrated over the tree of life. Nucleic Acids Res. 2015; 43:D447–52. https://doi.org/10.1093/nar/gku1003 [PubMed]
- 46. Shannon P, Markiel A, Ozier O, Baliga NS, Wang JT, Ramage D, Amin N, Schwikowski B, Ideker T. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 2003; 13:2498–504. https://doi.org/10.1101/gr.1239303 [PubMed]
- 47. Opsahl T, Agneessens F, Skvoretz J, Socnet J, Networks W. Node centrality in weighted networks: Generalizing degree and shortest paths. Social Networks. 2010; 32:245–251. https://doi.org/10.1016/j.socnet.2010.03.006
- 48. Cukierski WJ, Foran DJ. Using betweenness centrality to identify manifold shortcuts. Proc IEEE Int Conf Data Min. 2008; 2008:949–58. https://doi.org/10.1109/ICDMW.2008.39 [PubMed]
- 49. Du Y, Gao C, Chen X, Hu Y, Sadiq R, Deng Y. A new closeness centrality measure via effective distance in complex networks. Chaos. 2015; 25:033112. https://doi.org/10.1063/1.4916215 [PubMed]
- 50. Tang Y, Li M, Wang J, Pan Y, Wu FX. CytoNCA: a cytoscape plugin for centrality analysis and evaluation of protein interaction networks. Biosystems. 2015; 127:67–72. https://doi.org/10.1016/j.biosystems.2014.11.005 [PubMed]
- 51. Bader GD, Hogue CW. An automated method for finding molecular complexes in large protein interaction networks. BMC Bioinformatics. 2003; 4:2. https://doi.org/10.1186/1471-2105-4-2 [PubMed]
- 52. Chen EY, Tan CM, Kou Y, Duan Q, Wang Z, Meirelles GV, Clark NR, Ma’ayan A. Enrichr: interactive and collaborative HTML5 gene list enrichment analysis tool. BMC Bioinformatics. 2013; 14:128. https://doi.org/10.1186/1471-2105-14-128 [PubMed]
- 53. Rosenbloom KR, Dreszer TR, Long JC, Malladi VS, Sloan CA, Raney BJ, Cline MS, Karolchik D, Barber GP, Clawson H, Diekhans M, Fujita PA, Goldman M, et al. ENCODE whole-genome data in the UCSC genome browser: update 2012. Nucleic Acids Res. 2012; 40:D912–17. https://doi.org/10.1093/nar/gkr1012 [PubMed]
- 54. Kamburov A, Cavill R, Ebbels TM, Herwig R, Keun HC. Integrated pathway-level analysis of transcriptomics and metabolomics data with IMPaLA. Bioinformatics. 2011; 27:2917–18. https://doi.org/10.1093/bioinformatics/btr499 [PubMed]