3.1 Identification of DEGs using GEO datasets
GEO datasets GSE58840, GSE60335, GSE58878, GSE26465, GSE73935 and GSE54772, containing expression profiles of taxol-sensitive and resistant cell lines, were downloaded using R package “GEOquery”. The study characteristics and sizes of selected datasets were described in Table 1. 226 up-regulated genes and 214 down-regulated genes were found in taxol-resistant SKOV3 cells in microarray GSE58840, GSE60335 and GSE58878 (Figure 1A and Supplementary Table 1), while 494 up-regulated genes and 451 down-regulated genes were identified in taxol-resistant A2780 cell line in microarray GSE26465 (Figure 1B and Supplementary Table 1). 150 DEGs were screened from the GSE73935 dataset, including 71 up-regulated genes and 79 down-regulated genes in taxol-resistant OV90 cell line (Figure 1C and Supplementary Table 1). Additionally, 74 and 48 genes were up-regulated and down-regulated respectively in taxol-resistant SKOV3 cell line in dataset GSE54772 (Figure 1D and Supplementary Table 1). The overlapping up-regulated and down-regulated genes were obtained from the intersection of DEG datasets identified above (Figure 1E-F).
Table 1. mRNA sequencing datasets containing taxol-sensitive/resistant cell lines
3.2 Survival analysis
To explore if DEGs identified in taxol-sensitive and -resistant ovarian cancer cell lines are related to PFI and OS of ovarian cancer patients, samples from TCGA-OV with recurrence and therapy information are selected and analyzed (Supplementary Table 2). To each DEG identified, correlation of the expression of this gene with PFI and OS was evaluated with Kaplan-Meier method. Samples were divided into a high-expression group and a low-expression group according to the cutoff value of the certain gene, which was calculated by maximally selected rank and statistics method using R package “survminer”. For genes that show statistically significant differences in OS and PFI survival analysis, we determined whether the difference of their expression in sensitive and drug-resistant cell lines were consistent with the differences shown in survival analysis. If a specific gene had a higher expression level in taxol-resistant cell line, the survival of patients with this highly expressed gene should be poorer. Interferon Stimulated Gene 15 (ISG15), synuclein alpha (SNCA) and RIPK2 were up-regulated in taxol-resistant cell lines, while their high expression were also correlated with shorter OS and PFI in TCGA-OV dataset. Phospholipase C gamma 2 (PLCG2), ras homolog family member U (RHOU), tribbles pseudokinase 2 (TRIB2) and elongator acetyltransferase complex subunit 3 (ELP3) had low expression in taxol-resistant cell lines and their high expression was related to better survival in TCGA-OV dataset (Figure S1-6).
Dataset GSE30161, GSE32062 and GSE63885 which contained clinical information of patients with serous ovarian cancer including OS and PFI were selected to further validate effects of the expression of ISG15, SNCA, RIPK2, PLCG2, RHOU, TRIB2 and ELP2 on patients’ chemosensitivity to taxol treatment (Table 2 and Supplementary Figure 1-6). In all three datasets, high expression group and low expression group of RIPK2 showed significant difference in survival analysis with both OS and PFI, suggesting that high expression of RIPK2 was a risk factor of survival in patients with serous ovarian cancer (Figure 2).
We further validated the prediction value of RIPK2 expression by dividing TCGA-OV cohort into four groups based on patients’ RIPK2 expression level and if taxol was used during treatment. Survival analysis were carried out and we found that patients with low RIPK2 expression and taxol treatment showed the longest OS, while those with low RIPK2 but were not treated by taxol showed the shortest OS. Meanwhile, there were no significant difference of OS in patients with high expression of RIPK2, no matter they were treated by taxol or not.
Table 2. mRNA sequencing datasets containing overall survival and progress free interval of serous ovarian patients treated with taxol
3.3 Correlations of RIPK2 gene expression with taxol resistance in CCLE
CCLE contains a large panel of human cancer cell lines and their pharmacological profiles, in which the gene expression files and IC50 values to taxol of 21 ovarian cancer cell lines were included (Supplementary table 3). By dividing the expression level of RIPK2 in each cell line by the expression level of GAPDH in the same cell line, the expression of RIPK2 was normalized. Correlation of RIPK2 expression and IC50 of taxol was analyzed by R package “corrplot” with Spearman method. The correlation coefficient was 0.46 (p<0.05), indicating that higher expression of RIPK2 elevated taxol resistance of multiple ovarian cancer cell lines.
3.4 RIPK2 co-expression network in ovarian cancer
In order to gain further insight of the biological function of RIPK2 in development of taxol resistance in ovarian cancer, the co-expression mode of RIPK2 in serous ovarian cancer patients treated with taxol was analyzed. In dataset TCGA-OV, 341 genes were found to show a significant positive co-expression pattern with RIPK2, while no gene showed a negative co-expression. The expression of 704 genes had a positive correlation with RIPK2 expression in GSE30161 and 1706 had a negative correlation. 45 genes positively co-expressed while 37 negatively co-expressed with RIPK2 in GSE32063. 17 genes were found in dataset GSE63885 to have a positive co-expression pattern with RIPK2 and 4 genes had a negative co-expression pattern. A total description of the co-expressed genes was detailed in Supplementary Table 4.
Function analysis were made using the intersections of every two datasets as input. Significant GO term showed that RIPK2 co-expressed genes from multiple datasets participate mainly in cell adhesion molecule binding, positive regulation of cytokine production and focal adhesion(Figure 3A-C). KEGG pathway analysis showed an enrichment in NF-kappa B signaling pathway, NOD-like receptor signaling pathway and Ubiquitin mediated proteolysis pathway. (Figure 3D and Supplementary Table 5).
A PPI network of RIPK2 related genes were created on the basis of information from the STRING database, which further illustrates the connection of these co-expressed genes from the protein level. The average aggregation coefficient was 0.508, and the enrichment p<1.0e-3 (Figure 3E).
3.5 Genomic alterations of RIPK2 in ovarian cancer
cBioPortal tool was used to determine the alterations of RIPK2 in ovarian cancer patients who were treated with taxol in TCGA-OV database. Alterations occurred in 26 of 252 samples (10%), including 1 missense mutation (0.4%), 7 amplification (3%), 21 mRNA upregulation (8%) and 4 mRNA downregulation (2%) (Figure 4A). RIPK2 amplification results of high expression of RIPK2, which may relate to taxol resistance. Thus, AMP is the most common type of RIPK2 copy number alteration (CNA) in ovarian cancer (Figure 4B). Meanwhile there was significant difference of amplification of oxidative stress induced growth inhibitor family member 2 (OSGIN2), nibrin (NBN), Ras-Related Protein Rab-2A (RAB2A) and calbindin 1 (CALB1), etc in RIPK2-altered and -unaltered group (Figure 4C and Supplementary Table 6). Moreover, the mutation frequency of ArfGAP with SH3 domain, ankyrin repeat and PH domain 1 (ASAP1), ATP/GTP binding protein 1 (AGTPBP1), frizzled class receptor 7 (FZD7), HECT and RLD domain containing E3 ubiquitin protein ligase 5 (HERC5), KIAA0232, Mitogen-Activated Protein Kinase Kinase Kinase 10 (MAP3K10), PATJ crumbs cell polarity complex component (PATJ), PDGFA associated protein 1 (PDAP1), xin actin binding repeat containing 1 (XIRP1) was significantly associated with the alteration of RIPK2 (Figure 4D and Supplementary Table 7).
3.6 Evaluation the difference of immune cell infiltration
The immune infiltration of 64 types of immune cells, including adaptive and innate immune cells, hematopoietic progenitors, epithelial cells, and extracellular matrix cells, were evaluated by R package “xCell”, using the method of ssGSEA in ovarian cancer tissues. In dataset TCGA-OV, the infiltration of Macrophages M1, Melanocytes and plasmacytoid dendritic cells (pDC) was positively related to the expression of RIPK2, while the infiltration of Neurons was negatively related to the expression of RIPK2 (Figure 5A and Supplementary Table 8-9). Meanwhile, CD8+ naive T-cells , common lymphoid progenitors (CLP), CD4+ memory T-cells, smooth muscle and Hematopoietic stem cells (HSC) infiltrated more when RIPK2 expression levels were higher, but immature dendritic cells (iDC), Neurons, Basophils, Class-switched memory B-cells, mensenchymal stem cells(MSC), microvascular endothelial cells, natural killer T-cells(NKT), pro B-cells, Pericytes, Melanocytes, Mast cells, CD4+ T-cells, Plasma cells, MEP, lymphatic endothelial cells, Chondrocytes, pDC, endothelial cells, Myocytes and CD4+ central memory T-cells infiltrated less (Figure 5B andSupplementary Table 8-9). The infiltration of dendritic cells (DC) was positively correlated with the expression of RIPK2 and the infiltration of Mast cells was negatively correlated with the expression of RIPK2 in dataset GSE32063 (Figure 5C and Supplementary Table 8-9). In dataset GSE63885, melanocyte’s infiltration was high when RIPK2 had higher expression, while the infiltration of neurons and HSC was low (Figure 5D and Supplementary Table 8-9).