Machine learning and bioinformatics models to identify gene expression patterns of ovarian cancer associated with disease progression and mortality

Journal Pre-proofs Machine learning and Bioinformatics Models to Identify Gene Expression Patterns of Ovarian Cancer Associated with Disease Progression and Mortality Md. Ali Hossain, Sheikh Muhammad Saiful Islam, Julian M.W. Quinn, Fazlul Huq, Mohammad Ali Moni PII: DOI: Reference:

S1532-0464(19)30232-1 https://doi.org/10.1016/j.jbi.2019.103313 YJBIN 103313

To appear in:

Journal of Biomedical Informatics

Received Date: Revised Date: Accepted Date:

16 February 2019 20 September 2019 13 October 2019

Please cite this article as: Hossain, d.A., Saiful Islam, S.M., Quinn, J.M.W., Huq, F., Moni, M.A., Machine learning and Bioinformatics Models to Identify Gene Expression Patterns of Ovarian Cancer Associated with Disease Progression and Mortality, Journal of Biomedical Informatics (2019), doi: https://doi.org/10.1016/j.jbi. 2019.103313

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

© 2019 Published by Elsevier Inc.

Machine learning and Bioinformatics Models to Identify Gene Expression Patterns of Ovarian Cancer Associated with Disease Progression and Mortality Md. Ali Hossaina,b,∗, Sheikh Muhammad Saiful Islamc,∗, Julian M.W. Quinnd , Fazlul Huqe , Mohammad Ali Monid,e,∗∗ a

Dept of CSE, Manarat International University, Dhaka-1212, Bangladesh b Dept of CSE, Jahangirnagar University, Savar, Dhaka, Bangladesh c Dept. of Pharmacy, Manarat International University, Dhaka-1212, Bangladesh d Bone biology divisions, Garvan Institute of Medical Research, Sydney, NSW 2010, Australia e The University of Sydney, School of Medical Sciences, NSW 2006, Australia

Abstract Ovarian cancer (OC) is a common cause of cancer death among women worldwide, so there is a pressing need to identify factors influencing OC mortality. Much OC patient clinical data is publicly accessible via the Broad Institute Cancer Genome Atlas (TCGA) datasets which include patient age, cancer site, stage and subtype and patient survival, as well as OC gene transcription profiles. These allow studies correlation of OC patient survival (and other clinical variables) with gene expression to identify new OC biomarkers to predict patient mortality. We integrated clinical and tissue transcriptome data from patients available from the TCGA portal. We determined OC mRNA expression levels (compared to normal ovarian tissue) of 41 genes already implicated in OC progression, and assessed their ability to use their OC tissue expression levels predicts patient survival. We employed Cox Proportional Hazard regression models to analyse clinical factors and transcriptomic information to determine relative effects on survival that is associated with each factor. Multivariate analysis of combined data (clinical and gene mRNA expression) found age and ovary tumour site significantly correlated with pa∗ ∗∗

Joint First Author Corresponding Author Email address: [email protected] (Mohammad Ali Moni)

Preprint submitted to BioRxiv

October 23, 2019

tient survival. The univariate analysis also confirmed significant differences in patient survival time when altered transcription levels of TLR4, BSCL2, CDH1, ERBB2, and SCGB2A1 were evident, while multivariate analysis that considered the 41 genes simultaneously revealed a significant relationship of mortality with TLR4, BSCL2, CDH1, ERBB2 and PTPRE genes. However, analyses that considered all 41 genes with clinical variables together identified genes TLR4, BSCL2, CDH1, ERBB2, BRCA2 and SCGB2A1 as independently related to mortality in OC. These studies indicate that the latter genes influence OC patient survival, i.e., expression levels of these genes provide mechanistic and predictive information in addition to that of the clinical traits. Our study provides strong evidence that these genes are important prognostic indicators of patient survival that give clues to biological processes that underlie OC progression and mortality. Keywords: Ovarian cancer, Clinical factors, Gene expression, Survival analysis, RNA seq, Molecular pathways

1

2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20

1. INTRODUCTION Ovarian carcinoma (OC) accounts for the great majority of cases of cancer affecting the ovaries. It is the fifth most common cause of cancer deaths in women with an estimated 239,000 new cases and 152,000 deaths worldwide annually [1] and remains the leading cause of death from gynecologic malignancy [2] with five year mortality rate in the United States of 35% in 14,000 annual deaths [3]. Current standard treatment mainly comprises platinum-based chemotherapeutics, however the most prevalent OC type, high grade serous carcinoma (HGSC), often shows high chemotherapy resistance as well as greatly altered genome and transcriptome. High throughput gene sequence and gene expression data from patient tissues emerging as an important research resource to identify the key intrinsic factors that affect cancer behaviour, with numerous genetic mutations shown to be associated with OC development and progression [4]. Such genes are of wide interest for their potential prognostic power and as possible drug targets. While tumour-associated coding gene mutations can affect both gene product function and levels, it is their gene transcript levels that give more direct information regarding levels of their respective protein gene products [5, 6]. The latter is important in determining invasiveness, metastasis establishment, immune evasion and growth at metastatic sites, as well as indirect 2

21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58

pathological effects (via secreted factors) on other tissues that impair health. Clearly, OC gene expression patterns will influence patient fate and their quantification may also yield clues to important cellular pathways underlying disease processes. Tumour tissue transcript levels determined by RNA sequencing (RNAseq ) are thus becoming commonly used to study tumour clinical features [7]. Earlier work on OC has focused on influence of clinical features of the patient on outcomes [8, 9, 10, 11], which have enabled identification of significant risk factors such as age at diagnosis, age at menopause, disease stage and histological grade, tumour size, type of therapy received and family history of OC. In addition, new transcriptomic data from tumour tissue has identified a number of potential gene expression markers for patient progression [12, 13, 14]. Besides using survival analysis, some other methods have been applied for identifying prognostic biomarker genes in various cancers. These include prediction analysis for microarrays (PAM) [15], gene co-expression networks (GCN) [16] from microarray data, logistic regression along with odds ratio (OR) [17, 18] from DNA methylation, RNA expression and microarray data [19]. However, generally less well studied until recently [20] is how such tumour transcript levels may interact with clinical factors to affect cancer outcomes, which is a crucial consideration in determining how to use and interpret the transcript data [21]. Publicly available tumour transcript datasets are increasing in number along with clinical datasets, and this makes feasible large scale investigations of interaction between transcript levels and clinical trait. In particular, studies on combined datasets (containing both types of data from the same patient) can be particularly powerful tools for revealing important features of gene function in the tumour context, and informing the clinical use of such data. We thus employed the rich database available from The Cancer Genome Atlas (TCGA) project, a collaboration between the National Cancer Institute (NCI) and National Human Genome Research Institute (NHGRI). The TCGA project is supported by the Broad Institute to make available patient survival, clinical, and gene expression data so that researchers can analyse these either separately or in combination, to clarify the influences important to OC patient survival and which can be used as independent predictors. As part of this analysis we selected 41 genes of interest that have previously been validated as important in OC or OC models. The influence of clinical factors and disease marker gene expression on OC patient survival can thus be analysed by a number of methods related to 3

83

Cox Proportional Hazard (PH) modelling. Here, we employed standard Cox PH models [22] for univariate and multivariate analysis [23] as well as the penalized Cox PH model [24] for combined analysis [20]. Important clinical features with associations to the disease can be selected by consulting OC literature, then employed in such modelling [25]. In this study, we analysed post-diagnosis survival time. Clinical variables employed were patient age, ethnicity, anatomical site of cancer, histological grade of cancer, primary tumour site and neoplasm status with tumour. These were chosen from published studies that investigated associations between gene expression, clinical factors and OC patient survival. We gathered survival, clinical, and gene expression data from the TCGA datasets separately and combined them in order to assess the joint role of genetic and clinical factors. Our analysis followed the methodology of Ahmed et al. [20] and Xu and Moni (2015) [26] who used Cox PH regression modelling for the analysis of post-diagnosis survival time. We performed an analysis to identify most significant genes among selected 41 genes selected from literature searches as well as clinical variables, affecting patients survival in OC. To do so, we used product-limit estimators to estimate survival function for each gene separately, then used log rank test to detect which genes differ in expression levels significantly between altered and not altered group. Finally we performed univariate Cox PH regression analysis to determine each gene’s likelihood of contribution to deaths and two multivariate Cox PH regression analyses, the first taking consideration of 41 genes simultaneously and with the second considering all 41 genes and clinical variables to determine most significant genes and clinical variables in context of likelihood of risk of death.

84

2. METHODS & MATERIALS

59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82

85 86 87 88 89 90 91 92 93

In this study, our multi-step analysis procedure that we developed and applied is shown in Fig.1. We integrated transcriptomic and clinical data of patients with OC for selected biomarker genes obtained from literature and OMIM Database and applied univariate, multivariate and combined Cox PH regression analysis to identify significant biomarker genes that affect survival of OC patients. In order to identify important relevant protein protein interactions (PPI), biological pathways and ontological terms, we further performed PPI network analysis, GO ontology analysis (Biological pathway) as well as KEGG pathway analysis.

4

Figure 1: The multi-stage analysis methodology was used in this study. TCGA data (transcriptomic and clinical data) of OC were downloaded from cBioPortal and then combined transcriptomic and clinical data of patients with OC for selected biomarker genes obtained from literature and OMIM Database. Then we applied univariate, multivariate and combined Cox Proportional Hazard regression analysis to identify significant biomarker genes and clinical variables that affect survival of OC patients. We constructed protein-protein interaction network around selected 41 genes to identify any significant interaction among them which may play a pivotal role in survival of OC patients as well cancer progression, also we executed GO annotation analysis for enrichment relevant to biological processes. Finally we assessed the OC cancer related important pathways in which the prognostic genes that we found, are implicated.

94 95 96 97

2.1. Data We collected the RNAseq data for this study from TCGA genome data analysis centre (http://gdac.broadinstitute.org/) which is an interactive data system for researchers to search, upload, download, and analyse cancer ge5

98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135

nomic data [27]. Since our goal was to explore a particular point of interest, survival analysis of OC on clinical and genetic factors, we retrieved the anonymized clinical data and RNAseq data for OC (Ovarian Serous Cystadeno Carcinoma TCGA, Provisional) from the cBioPortal [28]. In the Clinical dataset there are 577 cases with 87 features. Cases that had RNAseq gene expression data included 535 cases with 11885 genes. We employed six clinical factors (ethnicity, anatomical site of cancer, histological grade of cancer, primary tumour site, and neoplasm status with tumour) along with 41 genes found in above 11885 genes. These genes are part of gene set of the OMIM database[29] as among the database repositories that delineate all currently known disease genes and their associated disorders, the OMIM database is the best-curated [30], and is frequently cited as significant factors in the OC literature (Table 1) [31, 32, 33, 34, 29]. Using these data we investigated a single outcome variable, namely OC-specific survival. 48 OC patient records that did not include clinical information were excluded from this analysis. We matched patient ID in both clinical and RNAseq dataset and identified 529 patients with data available for both. Among the clinical variables six clinical variables given above were considered [39, 32]. Tumour histological subtype and age at initial diagnosis were collected from pathology reporting. OC stage was recorded according to American Joint Committee on Cancer (AJCC) staging classification [40]. Patient survival data was taken from overall number of months of patient survival but converted to days of survival. Normal and tumour samples were identified using the TCGA barcode; two digits at position 14-15 of the barcode denote the sample type. Tumour type spans from 01 to 09, normal type from 10 to 19 and control samples from 20 to 29. Fold change (FC) was used as a common approach in differential gene expression analysis between two conditions [41]. For example, if gene expression read counts (determined by RNAseq) for a patient was 60 while that of a normal patient was 30, patient gene FC was 2. FC is most suitable when the gene expression distribution is symmetric. However, in RNAseq analyses, expression levels are modelled by the discrete counts of reads mapped to a known gene in a reference genome. Poisson and Negative binomial distribution assumptions are taken for reading counts. When the abundance rate of a particular gene expression level is very low, read count distributions modelled by the Poisson or Negative binomial are skewed to the right. Thus, as using FC as a measure of differential expression may not be appropriate in such cases we transformed the gene expression value using a standardizing transformation and calculated z-scores 6

Table 1: Ovarian cancer associated significant 41 genes and their sources.

Gene Name BMP15, FMR1, FOXL2, CDH1, NR5A1, OPCML, DIAPH2, RAD51C, RAD51D, MSH2, FSHR, PARK2, POF1B, RRAS2, SF1, BRCA1, BRCA2, CTNNB1, ERBB2, PIK3CA and MSH2 BRCA1, BRCA2, BSCL2, CTNNB1, ERBB2, MSH2, MUC16, PIK3CA, KRAS, TP53 KRAS, PIK3CA, and TP53

136 137 138 139 140

Resource Ref. Significant gene of OC from OMIM database [29] https://www.omim.org/

http://www.cancer-genetics.org/X1003.htm

[33]

Mutation profile and clinical outcome of mixed endometrioid-serous endometrial carcinomas are different from that of pure endometrioid or serous carcinomas KLK6, CP, SLPI , EZH2, Predicting biomarkers for ovarian cancer usHPN, SCGB2A1 ing gene-expression microarrays TLR4 The inflammatory microenvironment in epithelial ovarian cancer: a role for TLR4 and MyD88 and related proteins WFDC2 HE4 (WFDC2) gene overexpression promotes ovarian tumor growth RXRG, DVL1, MEF2B, The impact on high-grade serous ovarian MEF2C, BSCL2, CD36, cancer of obesity and lipid metabolismRXRG, E2F1 and TLR4 related gene expression patterns: the underestimated driving force affecting prognosis

[35]

[31] [36]

[37] [38]

for each expression value. For the expression from RNAseq experiments, the standard rule is to compute the relative expression of an individual gene in tumour samples using gene expression distribution in a reference population. A reference population was considered either all diploid tumours for the gene in question or, when available, normal adjacent tissue. The resulted value

7

141 142 143 144

(z-score) indicates the number of standard deviations away from the mean expression in the reference population. Here, we considered genes showing z > ±1.96 to be differentially expressed. We computed z-scores for RNAseq data using following formula: (Expression for gene X in tumor A - Mean expression for gene X in normal) (Standard deviation of expression for gene X in reference population) (1) We thus used the z-score values to define samples with ”altered”(Overexpress or Under-express) and ”normal” (unaltered) expression of a gene. We have assumed a sample to be altered if the z-score for that sample is equal to or higher than a specific threshold value such as z=2, as noted. We therefore define altered versus normal as follows: z=

145 146 147 148 149

z ≥ 2 => Over − express

(2)

z ≤ −2 => U nder − express

(3)

z < 2||z > −2 => normal

(4)

150 151

152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171

2.2. Methods We performed the following analysis related to survival of patients with OC. First we used product-limit estimator for estimation of survival function, then we performed log rank test to determine whether the survival functions of two different groups (patients with altered gene expression and patients with unaltered gene expression) exhibit statistically significant differences, then we used Cox PH regression models to determine the significance of genes and clinical factors on comparative risk of death. Finally, we performed functional analyses of our most significant genes found from our analysis. We selected important clinical and demographic variables affecting OC by literature review and similarly selected 41 genes having experimentally determined significance in context of OC. The expression z-score for each gene was identified as being in ”altered” or ”normal” categories based on a significance threshold value (z < 2) as noted in the Data section above. We performed Cox PH regression for every gene individually known as univariate regression as well as multivariate analysis taking all 41 genes simultaneously and finally multivariate regression on combined set of clinical and gene expression data. Survival analysis is used to estimate expected duration of time until one or more events happen, such as death in cancer and failure in mechanical systems [42]. Normally, survival analysis is carried in three steps: determining 8

172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191

192 193 194 195 196 197 198 199 200 201

202

203 204

time to event, formulating a censoring indicator for subject inclusion in the analysis, and the time to occurrence of event. Censoring in survival analysis is usually done in two ways, right censoring and left censoring. Right censoring occurs when a subject leaves the study before an event occurs, or the study ends before the event has occurred. Left censoring happens when the event of interest has already occurred before enrollment. This is very rarely encountered. Right censoring is again of two types. First one is Type I right censoring results from completely random dropout (e.g emigration) and/or end of study with no event occurrence and the second one, Type II right censoring, occurs with end of study after fixed number of events amongst the subjects has occurred. In survival analysis, one is interested in estimation of survival function of samples divided into subgroups or as whole. There are two types of estimation of it, one is non parametric in which no prior assumption of distribution of survival time is made and other is parametric estimation in which there are some assumptions along with assuming predefined distribution of survival time. As stated above, we used non parametric technique for estimation of survival function. There are several such techniques, we used product-limit estimator. In short product-limit (PL) estimator of the survival function is defined as follows: j Y dj ˆ S (tj ) = (1 − ) (5) n j i=1 Here Sˆ(tj ) is estimated survival function at time tj , dj is the number of events occurred at tj , and nj is the number of subjects available at tj . After estimating survival function, two or more groups can be compared using logrank test. For example, we used Log-rank test to detect the most significant genes in the case of patients survival time in altered versus unaltered groups in context of gene expression. The null hypothesis is following: H0 : Survival function for patients with altered gene expression is not different from the patients with normal (unaltered) gene expression HA : Survival functions are different for these two groups. Symbolically these can be written H0 : Saltered (t) = Snot altered (t) HA : Saltered (t) 6= Snot altered (t)

(6) (7)

Survival analysis methods can also be extended to assess several risk factors simultaneously similar to multiple linear and multiple logistic regression 9

205 206 207 208 209 210 211 212 213 214 215 216 217

analysis. One of the most popular regression techniques for survival analysis is Cox PH regression, which is used to relate one or several risk factors or exposures, considered simultaneously, to survival time. In a Cox PH regression model, the measure of effect is the hazard rate, which is the risk of failure (i.e., the risk or probability of suffering the event of interest), given that the participant has survived up to a specific time. Using Cox PH regression, first we performed univariate survival analysis by selecting each gene separately, then conducted a multivariate survival analysis taking all 41 genes simultaneously and finally fit a Cox PH model on all of selected six clinical factors and 41 genes combined, modelling the hazards of having the event under investigation (Ovarian Serous Cystadeno carcinoma, in this case), using an undetermined baseline hazard function and an exponential form of a set of covariates. Mathematically we can write the model as following: h(t|Xi ) = h0 (t)exp(β T Xi )

(8)

230

Whereas h(t|X) is the hazard function conditioned on a subject i with covariate information given as the vector xi , h0 (t) is the baseline hazard function which is independent of covariate information, and represents vector of regression coefficients to the covariates correspondingly. We have calculated the hazard ratio (HR) based on the estimated regression coefficients from the fitted Cox PH model to determine whether a specific covariate affects patient survival. The HR for a covariate xr can be expressed by the following simple formula exp (βr ). Thus, the HR for any covariate can be calculated by applying an exponential function to the corresponding (βr ) coefficient. Pathways and gene ontology (GO) for these genes were analysed using KEGG pathway database (http : //www.genome.jp/kegg/pathway.html) and enriched using (http : //amp.pharm.mssm.edu/Enrichr/enrich), a web based software tool.

231

3. Results and Discussion

218 219 220 221 222 223 224 225 226 227 228 229

232 233 234 235 236 237

6 clinical factors (including age) and mRNA expression data for 41 genes were employed in our analysis. All clinical factors were categorical except age. Age distribution of patients is shown in Figure 3. Average age of the patients at the time of their diagnoses was 59.7 years, with a range between 41 and 89 years old. The descriptive summary statistics of these factors shown in Table 2. From Table 2, we observe that the OC stage variable 10

Diagnosis Age 150 100 50 0 <=30 30

40

50

60

70

80

90

>90 NA

Figure 2: Distribution of Age at Onset of Diagnosis Table 2: Descriptive Statistics of Clinical Predictors

238 239 240 241 242 243 244 245

has the highest number (10) of categories and 71.03% of all patients are of stage IIIC. The next largest group is stage IV with 15.36% cases. There are 8 categories for histology type with the highest percentage of patients (84.67%) from was G3, i.e., poorly differentiated, while the second largest group (12.02%) of patients had moderately differentiated grade. Regarding ethnicity, most patients (85.27%) are from the European-caucasian descent population, while African American is next largest group with 5.89%. In the case of OC anatomical site, most (72.53%) women had bilateral cancers, while

11

246 247 248

249 250 251 252 253 254 255 256 257

cancer patients with left and right ovary are 14.65% and 12.82% respectively. Most women had cancer in the ovary itself (99.13%), the remaining minor percentage of women had tumours located in omentum and/or peritoneum. 3.1. Survival pattern for gene expression data We have estimated survival function for altered and unaltered group for each of the 41 genes by applying product limit (PL) estimator. We then compared estimated survival function for altered and unaltered group using log rank test. Those genes for which there is statistically significant difference are shown below. The significant role for these genes is indicated by their pvalues in differential survival pattern when comparing their expression level in two categories (altered and unaltered). Surprisingly, the p53 gene (TP53), which is known for its implication in cancer in general, has p-value greater than 0.05. From Figure 3, we can see that patients having altered expression

Figure 3: Survival pattern of altered (Over-express and Under-express) and normal(nonaltered) groups for BSCL2, CDH1, ERBB2, SCGB2A1 and TLR4 258 259

of TLR4, BSCL2, CDH1, ERBB2, and SCGB2A1 genes are less likely to 12

260 261 262 263

264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280

281 282 283 284 285 286 287 288 289 290 291 292 293 294 295

survive compared to the non-altered group. Note that, the red line in the graphs indicates normal gene expression, the blue line indicates the underexpression of gene expression group and green line denotes over-expression of gene expression group. 3.2. Modeling the hazard risk on the RNA-Seq data One can measure the relative likelihood of risk of death in OC for each gene separately as well as for all genes simultaneously. Consequently one can determine which genes are most significant in case of survival of patients in OC. For these purposes, Cox PH regression model is used. We thus considered both univariate (separately for every gene) and multivariate (incorporating rest of the genes) for Cox PH model on each of the 41 genes. Table 3 shows the estimated Coefficients (β), with corresponding hazard ratios (HR), and p-values from those analyses. We noted that the p-values for genes TLR4, BSCL2, CDH1 and ERBB2 showed that their expression profiles have statistically significant association with OC patient survival in both univariate and multivariate analyses. In contrast, gene SCGB2A1 showed significant association only in univariate analysis, while PTPRE showed significance only in multivariate analysis. Surprisingly, even though the KRAS gene has previously been shown by experimentation to have a significant role in OC [43], we did not find its significance in our analysis, both in univariate and multivariate analysis. 3.3. Modeling hazard on the combined model containing both clinical and RNAseq data We have performed multivariate Cox PH regression analysis with both clinical and RNA-Seq data simultaneously. We have presented estimated regression coefficients in Table 4 and 5 along with hazard ratios (HR), z-values and corresponding p-values in our combined Cox PH model. Here, n= 497, number of events (deaths)= 178 after deleting 32 observations due to missing data values. Table 4 and table 5 shows the summary of the Cox Proportional Hazard Model result for Clinical data and RNA-Seq data, respectively. We observe from Table 4 that if the patient’s age increases one-year, likelihood of hazard increases 1.03 times. We also observed that AMERICAN INDIAN and BLACK OR AFRICAN AMERICAN patients had a greater risk of death. Notably, we found no evidence that ethnicity, histologic grade and anatomical OC site have any statistically significant association with hazard. Surprisingly, our analysis indicates that patients with a right side ovary 13

Table 3: Summary of Univariate and Multivariate Cox Proportional Hazard Model for mRNA Seq data. Univariate Gene

Overexpressed

ᵦ

HR

PValue

Multivariate

Underexpressed

ᵦ

HR

PValue

Overexpressed

ᵦ

HR

PValue

Underexpressed

ᵦ

HR

PValue

AKT1

-1.62E-01 8.51E-01

0.464 2.13E-01 1.24E+00 0.4058 -1.45E-01 8.65E-01 0.5859 1.57E-01 1.17E+00 0.6080

BMP15

-3.88E-02 9.62E-01

0.908 3.70E-01 1.45E+00 0.3696 1.64E-01 1.18E+00 0.6685 3.89E-01 1.48E+00 0.4379

BRCA1

-1.19E-02 9.88E-01

0.973 6.28E-01 1.87E+00 0.5317 -3.23E-01 7.24E-01 0.4897 5.18E-01 1.68E+00 0.6204

BRCA2

-8.56E-02 9.18E-01

0.754 -4.02E-01 6.69E-01 0.4247 2.28E-01 1.26E+00 0.4627 -1.16E-01 8.90E-01 0.8259

BSCL2

2.57E-02 1.03E+00

0.946 9.40E-01 2.56E+00 0.0000 -5.06E-02 9.51E-01 0.9096 9.12E-01 2.49E+00 0.0004

CD36

1.68E-01 1.18E+00

0.813 1.88E-01 1.21E+00 0.4042 3.34E-01 1.40E+00 0.6574 8.19E-02 1.09E+00 0.7727

CDH1

4.87E-03 1.00E+00

0.989 5.93E-01 1.81E+00 0.0027 -2.25E-01 7.99E-01 0.5580 8.30E-01 2.29E+00 0.0008

CLDN3 CP

4.08E-01 1.50E+00 … …

0.566 2.56E-01 1.29E+00 0.2794 7.76E-02 1.08E+00 0.9253 1.46E-01 1.16E+00 0.6614 … … … … -1.08E-01 8.98E-01 0.6489 -2.74E-01 7.60E-01 0.4658

CTNNB1

-2.44E-01 7.84E-01

0.427 -1.25E-01 8.83E-01 0.6063 -1.69E-01 8.44E-01 0.6422 -3.84E-01 6.81E-01 0.1907

DIAPH2

-2.31E-01 7.94E-01

0.451 -4.71E-01 6.24E-01 0.5077 -2.52E-01 7.77E-01 0.4891 1.33E-02 1.01E+00 0.9871

DVL1

-3.21E-02 9.68E-01

E2F1

6.27E-02 1.06E+00

0.892 2.11E-02 1.02E+00 0.9164 -6.55E-03 9.93E-01 0.9812 8.35E-02 1.09E+00 0.7266 … … … … … … 0.763 2.68E-01 1.31E+00 0.3093

ERBB2

6.79E-01 1.97E+00

0.017 -3.26E-01 7.22E-01 0.0753 5.90E-01 1.80E+00 0.1138 -7.16E-01 4.89E-01 0.0030

EZH2

-3.46E-01 7.08E-01

0.260 -2.17E-03 9.98E-01 0.9949 -3.94E-01 6.74E-01 0.3278 9.28E-02 1.10E+00 0.8149

FMR1

-1.27E-01 8.81E-01

FOXL2

-1.47E-01 8.63E-01

0.565 5.22E-01 1.69E+00 0.0792 -1.48E-01 8.62E-01 0.5522 3.33E-01 1.40E+00 0.3826 … … … -3.99E-02 9.61E-01 0.8750 … … … 0.505

FSHR

-8.56E-02 9.18E-01

HPN KLK6

1.95E-01 1.22E+00 … …

0.754 -4.02E-01 6.69E-01 0.4247 -1.47E-01 8.63E-01 0.6508 -1.93E-01 8.24E-01 0.6494 … … … … … … 0.348 2.36E-01 1.27E+00 0.3187 … … … … 2.89E-01 1.34E+00 0.1732 6.91E-02 1.07E+00 0.8303

KRAS

-4.35E-02 9.57E-01

0.777 -1.64E-01 8.49E-01 0.5368 2.01E-01 1.22E+00 0.2547 -1.40E-01 8.70E-01 0.6405

MEF2B

-7.57E-02 9.27E-01

0.732 1.53E-02 1.02E+00 0.9758 -7.60E-02 9.27E-01 0.7793 3.65E-01 1.44E+00 0.5332

MEF2C

4.92E-01 1.64E+00

0.094 -5.52E-01 5.76E-01 0.4364 -5.18E-02 9.50E-01 0.8800 3.13E-03 1.00E+00 0.9970

MSH2

-2.70E-01 7.64E-01 … …

OPCML

1.41E-01 1.15E+00

0.380 1.73E-02 1.02E+00 0.9513 2.75E-01 1.32E+00 0.5023 -1.17E-01 8.90E-01 0.7401 … … … … 1.37E-01 1.15E+00 0.4727 3.23E-02 1.03E+00 0.9156 … … … … 1.37E-01 1.15E+00 0.4727 -2.99E-02 9.70E-01 0.9303 … … … … … … 0.619 5.26E-02 1.05E+00 0.9178

MUC16 NR5A1

…

…

PARK2

1.41E-01 1.15E+00

0.619 8.69E-02 1.09E+00 0.8472 2.00E-01 1.22E+00 0.4993 -9.02E-02 9.14E-01 0.8377

PIK3CA

3.84E-01 1.47E+00

POF1B

-4.64E-02 9.55E-01

0.134 -1.26E-02 9.87E-01 0.9720 -6.93E-02 9.33E-01 0.6677 -1.59E+01 1.30E-07 0.9918 … … … 0.733 2.55E-02 1.03E+00 0.9241

PTPRE

-1.50E+01 2.99E-07

0.990 -5.18E-02 9.50E-01 0.8145 7.52E-01 2.12E+00 0.0337 -7.36E-01 4.79E-01 0.1320

RAD51D

4.12E-01 1.51E+00

0.109 -3.55E-01 7.01E-01 0.3527 -9.66E-02 9.08E-01 0.8287 1.13E-01 1.12E+00 0.6519

RAD51C

2.39E-01 1.27E+00

0.479 4.56E-02 1.05E+00 0.8263 3.56E-01 1.43E+00 0.4809 -5.55E-02 9.46E-01 0.8261

RRAS2

-5.58E-01 5.72E-01

RXRG

-4.72E-01 6.24E-01

0.218 3.02E-01 1.35E+00 0.0548 -8.64E-01 4.21E-01 0.1120 1.86E-01 1.20E+00 0.3343 … … … -5.85E-01 5.57E-01 0.0955 … … … 0.125

SCGB2A1 2.39E+00 1.09E+01 SF1 -3.85E-01 6.80E-01

0.019 2.88E-01 1.33E+00 0.2492 2.21E+00 9.12E+00 0.0594 4.61E-01 1.59E+00 0.2157

TLR4

6.79E-01 1.97E+00

TP53

6.13E-02 1.06E+00

0.282 -6.52E-01 5.21E-01 0.5153 -7.02E-01 4.95E-01 0.1025 -4.62E-01 6.30E-01 0.6658 … … … … 2.25E-01 1.25E+00 0.3424 -3.90E-02 9.62E-01 0.9221 … … … … … … 0.011 7.16E-01 2.05E+00 0.0138 … … … … … … 0.951 2.34E-01 1.26E+00 0.8199

WFDC2

1.68E-01 1.18E+00

0.813 1.88E-01 1.21E+00 0.4042 1.59E+00 4.90E+00 0.2130 1.26E-01 1.13E+00 0.7582

SLPI

…

…

14

Table 4: Summary of the Cox Proportional Hazard Model for combined mRNA seq and Clinical data, showing Clinical data

ᵦ

Variables

HR

P-Value

2.78E+00 1.61E+01 0.00074 AMERICAN INDIAN OR ALASKA NATIVE -2.73E-02 9.73E-01 0.95317 ASIAN 7.37E-01 2.09E+00 0.01637 BLACK OR AFRICAN AMERICAN NATIVE HAWAIIAN OR OTHER PACIFIC ISLANDER 1.54E+00 4.68E+00 0.16029 2.60E-02 1.03E+00 3.30E-05 age 3.83E-02 1.04E+00 0.97248 anatomic_site Peritoneum ovary 6.95E-01 2.00E+00 0.40454 histologic_grade G2 8.01E-01 2.23E+00 0.33189 G3 2.65E+00 1.42E+01 0.08085 G4 3.36E-01 1.40E+00 0.83559 GB -1.91E-01 8.26E-01 0.84474 GX -3.71E-01 6.90E-01 0.05113 tumour_site Left 2.69E-01 1.31E+00 0.19474 tumour_site Right 1.22E+00 3.40E+00 0.25561 Stage IB -7.38E-01 4.78E-01 0.49448 Stage IC -1.58E+01 1.43E-07 0.99185 Stage IIA 6.35E-02 1.07E+00 0.93521 Stage IIB -1.18E+00 3.09E-01 0.00483 Stage IIC -9.63E-02 9.08E-01 0.90753 Stage IIIA -2.79E-01 7.57E-01 0.42995 Stage IIIB -3.10E-01 7.33E-01 0.08836 Stage IIIC NA NA NA Stage IV

296 297 298 299 300 301 302 303 304 305 306 307 308 309 310

tumour have a more likely to survive compared to patients with a bilateral tumour, but patients having left side and bilateral tumour exhibit same likelihood of risk of death. Tumour stage is an important factor affecting risk of death in many cancers, and consistent with this, we found that patients with Stage IIC OC are 0.31 times less likely to survive compared to stage IA, although there was no evidence that other OC stages exerted a significant influence on the likelihood of death compared to stage IA. We then conducted our hazard analysis using the 41 selected genes. Among these, only TLR4, BSCL2, CDH1, ERBB2, SCGB2A1 and BRCA2 expression levels showed a significant association to patient survival. Patients with under-expressed ERBB2 gene expression have respectively a 0.5 times lower risk of death compared to those with no such alterations. In contrast, patients with overexpressed TLR4, SCGB2A1 and BRCA2 gene expressions have respectively 1.93, 12.2, and 1.91 times higher rates of likelihood of death than those with no alteration in expression in these genes. Notable, under-expressed genes

15

Table 5: Summary of the Cox Proportional Hazard Model for combined mRNA seq and Clinical data, showing genes

311

BSCL2, CDH1 have 2.01, 2.55 times higher rates of likelihood of death than

16

312 313 314 315 316 317 318 319 320 321 322 323 324 325

normal groups and ERBB2 has 0.5 time less likelihood of death than normal groups. A Venn diagram of the significant genes summarises the relationships between altered gene expression and risk of death identified using different methodologies and is shown in figure 4. From univariate, multivariate analysis and Cox PH regression analysis on combined data, we found that 7 genes among the 41 studied showed significant association with risk of death; these were TLR4, BSCL2, CDH1, ERBB2, PTPRE, SCGB2A1 and BRCA2. Gene expression of SCGB2A1 was significant in both univariate and combined hazard models, that of BRCA2 only in combined hazard models, and PTPRE expression only in multivariate analysis. Notably, TLR4, BSCL2, CDH1, and ERBB2 gene expression were found significant in all of the three types of models. The latter suggests that knowledge of these gene expression levels does not add to information from the clinical data.

Figure 4: Venn diagram of significant genes found in univariate, multivariate and combined hazard model analysis

326 327 328 329 330 331 332 333 334

3.4. Pathway and functional correlation analysis of the significant genes We observed that nineteen significant pathways including MicroRNAs in cancer, endometrial cancer, bladder cancer, non-small cell lung cancer, and pancreatic cancer, prostate cancer are associated with the significantly regulated genes for Ovarian Cancer. Genes associated with these pathways and corresponding p-values are presented in the table (see Table 6 ). We also performed the biological process ontology enrichment analysis (see table 7) of these identified significant genes using the Enrichr software tool. We found 26 biological pathways associated with these significant genes 17

Table 6: Pathway associated with the selected 7 significantly associated genes with the Ovarian Cancer. Nineteen KEGG pathways were found using Enrichr for these genes. Genes associated with these pathways and corresponding Adjusted p-values are presented in this table. KEGG ID

Pathways Name

Adjusted P-Value

Genes AKT1, KRAS, ERBB2, TP53, PIK3CA, CDH1, CTNNB1 E2F1, AKT1, KRAS, ERBB2, TP53, RXRG, PIK3CA E2F1, AKT1, KRAS, MSH2, ERBB2, TP53, RXRG, BRCA2, PIK3CA, CDH1, CTNNB1, DVL1 E2F1, AKT1, KRAS, ERBB2, TP53, BRCA2, PIK3CA E2F1, AKT1, KRAS, ERBB2, TP53, PIK3CA, CTNNB1

hsa05213

Endometrial cancer

4.60E-08

hsa05223

Non-small cell lung cancer

7.26E-08

hsa05200

Pathways in cancer

1.64E-07

hsa05212

Pancreatic cancer

1.80E-07

hsa05215

Prostate cancer

1.11E-06

hsa05216

Thyroid cancer

4.86E-06 KRAS, TP53, RXRG, CDH1, CTNNB1

hsa05218

Melanoma

hsa05205

Proteoglycans in cancer

hsa05219

hsa04012 hsa05206 hsa03440 hsa04066

Bladder cancer Central carbon metabolism in cancer ErbB signaling pathway MicroRNAs in cancer Homologous recombination HIF-1 signaling pathway

8.30E-06 E2F1, AKT1, KRAS, TP53, PIK3CA, CDH1 AKT1, KRAS, RRAS2, ERBB2, TP53, PIK3CA, 1.11E-05 TLR4, CTNNB1 2.01E-05 E2F1, KRAS, ERBB2, TP53, CDH1

hsa04015

Rap1 signaling pathway

9.72E-03 AKT1, KRAS, PIK3CA, CDH1, CTNNB1

hsa04151

PI3K-Akt signaling pathway Pathogenic Escherichia coli infection Fanconi anemia pathway Adherens junction

1.17E-02 AKT1, KRAS, TP53, PIK3CA, TLR4, BRCA1

hsa05230

hsa05130 hsa03460 hsa04520

1.18E-04 AKT1, KRAS, ERBB2, TP53, PIK3CA 5.09E-03 5.36E-03 6.06E-03 6.70E-03

AKT1, KRAS, ERBB2, PIK3CA E2F1, KRAS, ERBB2, EZH2, TP53, BRCA1 RAD51D, RAD51C, BRCA2 AKT1, ERBB2, PIK3CA, TLR4

1.80E-02 CDH1, TLR4, CTNNB1 1.94E-02 RAD51C, BRCA2, BRCA1 3.34E-02 ERBB2, CDH1, CTNNB1

338

as shown in Table 7. We have investigated protein protein interaction network generated using STRING [44], a web-based visualization software resource. Most of the key genes are connected with each other through the PPI network. However some genes are not connected with the main network.

339

4. DISCUSSION

335 336 337

340 341 342 343 344 345 346

We used univariate, multivariate Cox PH ratio analysis using mRNA expression data, estimating the survival curve using product limit procedure and determining whether there is any statistically significant difference between the altered and un-altered groups using log rank test for each gene. This identified five significant genes (TLR4, BSCL2, CDH1, ERBB2, and SCGB2A1) in univariate, five (TLR4, BSCL2, CDH1, ERBB2, and PTPRE) in multivariate and six genes (TLR4, BRCA2, BSCL2, CDH1, ERBB2, and 18

Table 7: Gene Ontolgy terms with underlying Biological pathways associated with 7 significant genes in context of the Ovarian Cancer along. 26 Gene Ontology terms along with biological pathways, found significantly associated with these significant genes with adjusted p-values shown in this table. Go Term

Adjusted P-Value

Pathways Name

GO:0010628 positive regulation of gene expression

1.28E-07

GO:0045893 positive regulation of transcription, DNA-templated

2.09E-06

GO:0000732 strand displacement GO:0071681 cellular response to indole-3-methanol positive regulation of transcription from RNA GO:0045944 polymerase II promoter GO:0000731 DNA synthesis involved in DNA repair GO:0038128 ERBB2 signaling pathway GO:0043406 positive regulation of MAP kinase activity DNA damage response, signal transduction by p53 GO:0006978 class mediator resulting in transcription of p21 class mediator double-strand break repair via homologous GO:0000724 recombination intrinsic apoptotic signaling pathway in response to GO:0042771 DNA damage by p53 class mediator GO:0001934 positive regulation of protein phosphorylation

2.90E-05 5.23E-05

GO:0006351 transcription, DNA-templated GO:0030521 GO:0006289 GO:0006302 GO:0031052 GO:1900227 GO:0007283 GO:0014066 GO:0060907 GO:0030307 GO:0050702 GO:0000165 GO:0048015 GO:0008585

347

androgen receptor signaling pathway nucleotide-excision repair double-strand break repair chromosome breakage positive regulation of NLRP3 inflammasome complex assembly spermatogenesis regulation of phosphatidylinositol 3-kinase signaling positive regulation of macrophage cytokine production positive regulation of cell growth interleukin-1 beta secretion MAPK cascade phosphatidylinositol-mediated signaling female gonad development

5.94E-05 7.20E-05 9.24E-05 3.44E-04

Genes E2F1, MEF2C, HPN, KRAS, ERBB2, TP53, TLR4, PARK2, BRCA1 E2F1, MEF2C, FOXL2, TP53, BRCA2, CDH1, BMP15, BRCA1, CTNNB1, DVL1 RAD51D, RAD51C, BRCA2, BRCA1 CDH1, BRCA1, CTNNB1 E2F1, MEF2C, AKT1, MEF2B, FOXL2, TP53, TLR4, PARK2, BRCA1, CTNNB1, NR5A1 RAD51D, RAD51C, BRCA2, BRCA1 AKT1, KRAS, ERBB2, PIK3CA MEF2C, KRAS, ERBB2, EZH2

6.18E-04 TP53, BRCA2, BRCA1 6.70E-04 RAD51D, RAD51C, BRCA2, BRCA1 2.34E-03 MSH2, TP53, BRCA2 3.17E-03 AKT1, KRAS, ERBB2, DVL1 E2F1, MEF2C, MEF2B, ERBB2, EZH2, SF1, 3.91E-03 TP53, RXRG, PARK2, BRCA1, CTNNB1, NR5A1 4.07E-03 SCGB2A1, BRCA1, CTNNB1 4.07E-03 RAD51D, TP53, BRCA2 1.03E-02 MSH2, BRCA2, BRCA1 1.16E-02 BRCA2, BRCA1 1.16E-02 CD36, TLR4 1.19E-02 1.41E-02 1.61E-02 1.63E-02 1.84E-02 2.27E-02 2.51E-02 3.43E-02

E2F1, RAD51C, BRCA2, FSHR, WFDC2 AKT1, ERBB2, PIK3CA CD36, TLR4 AKT1, HPN, ERBB2 CD36, TLR4 MEF2C, MEF2B, KRAS, ERBB2 AKT1, ERBB2, PIK3CA BRCA2, FSHR

SCGB2A1) in combined analysis.

348 349 350 351 352 353 354 355 356 357 358 359

One study described the TLR4 gene as a novel prognostic biomarker in solid tumors [45]. It was found that the expression of the DAMP(damageassociated molecular patterns)-derived molecules were upregulated in the tumor or tumor microenvironment and that caused TLR4-related chronic inflammation, leading to carcinogenesis, cancer progression, and metastasis [45]. TLR4 signaling was reported to be related to numerous cancers, such as lung, ovarian, liver, pancreatic, colon and gastric cancer [45]. Many studies found TLR4 to be strongly associated with shorted poor prognosis [46], survival time [47, 36, 45], advance stage of high grade serous ovarian cancer [47], increased risk of developing OC [48]. It is also reported to be a drug candidate for OC by inhibiting cell growth and induction of apop-

19

Figure 5: PPI network of 41 selected genes of OC

360 361 362 363

tosis [49]. From KEGG pathway analysis [50] , TLR4 was found to be involved in PI3K/AKT/mTOR pathway, HIF-1 signaling pathway, Pathogenic Escherichia coli infection and Measles.TLR4 actively participates in HIF-1 signalling pathway which promotes hypoxia induced metastasis of OC [51].

364 365 366 367 368 369 370 371 372 373 374 375

BRCA2, one of the two most thoroughly characterized genetic risk factors for ovarian cancer, functions in the double strand DNA break repair pathway, and BRCA2 mutations are present in 5-15% of all ovarian cancer cases [52]. Several studies [53, 54, 55] identified mutated BRCA2 as posing significantly high risk for OC. Mutations in both BRCA1 and BRCA2 genes have been experimentally implicated in OC [55]. In addition, the PI3K/AKT/mTOR pathway is activated in approximately 70% of OC cases [56]. From a recent study, it was found that this pathway is also actively involved in Thyroid cancer [57]. We found from KEGG pathway analysis that the BRCA2 gene is principally associated with the Pathways in cancer, Pancreatic cancer, Homologous recombination, and Fanconi anemia pathways.

376 377

BSCL2 is a gene encoding an endoplasmic reticulum protein that believed

20

378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415

to be involved in lipid droplet morphology. It has been reported to be associated with progression-free and overall survival in high-grade ovarian serous carcinoma (HGOSC) [38]. One study has reported that the CDH1 gene is a significant biomarker of OC [58]. CDH1, a tumor suppressor gene and a member of the cadherin family, plays an important role in epithelial cell-cell adhesion and in maintaining normal tissue architecture and thus reduction of CDH1 expression may influence invasion and metastasis of cancer cells [59, 58]. It has also been reported that CDH1 promoter methylation status is significantly associated with the risk of OC incidence [58]. CDH1 is also expressed significantly as a bio-marker for suppressing gastric cancer [60]. Our KEGG pathways analysis suggested that this gene is may be involved in Adherens junction, Bacterial invasion of epithelial cells and Melanoma pathways, especially in Rap1 signalling pathway that plays dominant role in tumor cell migration and invasion [61] ERBB2 gene encodes the tyrosine kinase-type cell surface receptor HER2, which binds epidermal growth factor. ERBB2/HER2 expression is a cancer prognostic biomarker in OC patients [62]. From KEGG pathway analysis, this gene is found actively involved in Endometrial cancer, Non-small cell lung cancer, Pathways in cancer, Pancreatic cancer, Prostate cancer, Proteoglycans in cancer, Bladder cancer, Central carbon metabolism in cancer and ErbB signaling pathway. ErbB is a notable member of the extracellular growth factor receptor family (EGFR) and is reported to be implicated in many cancers including OC [63]. PTPRE is a member of protein kinase family is a signaling molecule regulating numerous cellular processes such as cell growth, differentiation, mitotic cycle, and oncogenic transformation. PTPRE activates Src family kinases through dephosphorylation of C-terminal inhibitory sites of Src, Yes, and Fyn which stimulates cell proliferation [64]. In a study, assessing the impact of obesity and lipid metabolism associated gene on serous ovarian cancer found PTPRE is significantly implicated in shorter-term survival [38]. Another study evaluating contribution of individual genes in a very large gene set identified PTPRE as one of 56 of the most significant genes implicated in OC [65]. SCGB2A1 (Secretoglobin Family 2A Member 1) is a highly differentially expressed gene in all grades of OC [66, 67]. This gene has also been significantly associated with breast cancer incidence [68] and hypotrichosis [69]. Gene Ontology (GO) annotations related to this gene include androgen recep21

426

tor signaling pathway. This pathway is critically involved in OC and can play a more important role in the pathogenesis of OC [70]. Among the 41 genes, we found 7 significant genes that affect the survival of OC whereas only 4 genes (TLR4, CDH1, ERBB2 and BRCA2) are strongly connected with the PPI network. Because the other three important genes namely BSCL2, SCGB2A1, PTPRE are not directly connected with other four significant genes and their directly connected neighbour proteins. But the former three genes interact with the later ones if we consider more than one intermediate proteins among the significant 7 genes. Along with connected proteins, these seven proteins participate in important OC related pathways and biological processes, thus exerting dominant role in survival of OC patients.

427

5. CONCLUSIONS

416 417 418 419 420 421 422 423 424 425

428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451

In this study, survival analysis of 577 OC patients revealed that out of 41 genes chosen for their previously idedntified involvement in OC, altered transcript levels of ten of these predicted reduced OC survival. Using product limit or Kaplan-Meier analyses, we found a significant difference in survival time between altered and non-altered genes. In cases of TLR4, BSCL2, CDH1, ERBB2, and SCGB2A1 genes, patients with these altered expressions of these five genes had significantly lower survival time than patients with non-alteration of these five genes. When taking consideration of mutual effect of all selected 41 genes in multivariate analysis, TLR4, BSCL2, CDH1, ERBB2, and PTPRE exhibited significant influence on survival of OC, suggesting that they may explain any survival alteration predicted by SCGB2A1 levels. TLR4, BSCL2, CDH1 and ERBB2 genes were notable for their high statistical significance in the Kaplan-Meier survival analysis, making these genes of particular interest, at least in this OC patient cohort. When we extended our analysis to evaluate which clinical factors and genes play dominant roles in determining survival time in OC patient, and considering interaction between clinical variables and genes, five genes including TLR4, BSCL2, CDH1, ERBB2, SCGB2A1, and BRCA2 were found to be the genes having most influence on survival, i.e., their influence was independent of clinical features so would be useful in patient survival prediction. Again, TLR4, BSCL2, CDH1 and ERBB2, seems the most promising candidate for further investigation since these genes are significant in all three types of analysis we performed. These genes may be candidates for therapeutic drug discovery itself or they may point the way to a crucial cell pathway that 22

458

influences patient survival Among the clinical variables, AMERICAN INDIAN, AFRICAN AMERICAN, patient age, and cancer stage status are associated with significantly more risk of death within 5 years in OC. In addition, we found that patients with left site tumour in ovary had less risk of death. Our approach can be used in case of other types of cancers to identify key genetic and clinical factors in patient survival.

459

Acknowledgements

460

Conflicts of interest

452 453 454 455 456 457

461

462

463 464 465 466

467 468

469 470 471

472 473 474 475

476 477 478

479 480 481 482

We have no conflicts of interest. References [1] J. Ferlay, I. Soerjomataram, M. Ervik, R. Dikshit, S. Eser, C. Mathers, M. Rebelo, D. Parkin, D. Forman, F. Bray, Cancer incidence and mortality worldwide: Iarc cancerbase no. 11. 2013, Lyon, France: International Agency for Research on Cancer (2014). [2] R. Seigel, D. Naishadham, A. Jemal, Cancer statistics, 2014, Ca Cancer J Clin 64 (2014) 9–29. [3] E. Gov, K. Y. Arga, Differential co-expression analysis reveals a novel prognostic gene module in ovarian cancer, Scientific reports 7 (2017) 4996. [4] R. Rahman, T. Islam, E. Gov, B. Turanli, G. Gulfidan, M. Shahjaman, N. A. Banu, N. H. Mollah, K. Y. Arga, M. A. Moni, Identification of prognostic biomarker signatures and candidate drugs in colorectal cancer: Insights from systems biology analysis (2018). [5] M. A. Hossain, T. A. Asa, J. M. Quinn, M. M. Rahman, F. Huq, M. A. Moni, Network-based genetic profiling, and therapeutic target identification of thyroid cancer, bioRxiv (2018) 480632. [6] M. A. Hossain, T. A. Asa, F. Huq, J. M. Quinn, M. A. Moni, A networkbased approach to identify molecular signatures and comorbidities of thyroid cancer, in: Proceedings of International Joint Conference on Computational Intelligence, Springer, pp. 235–246. 23

483 484 485

486 487 488

489 490 491

492 493

494 495

496 497 498 499

500 501

502 503 504 505

506 507 508 509

510 511 512 513 514

[7] M. A. Moni, P. Li`o, Network-based analysis of comorbidities risk during an infection: Sars and hiv case studies, BMC bioinformatics 15 (2014) 333. [8] Q.-Q. Zheng, P. Wang, R. Hui, A.-M. Yao, Prognostic analysis of ovarian cancer patients using the cox regression model, Chinese Journal of Cancer 28 (2009) 170–2. [9] M. A. Moni, H. Xu, P. Lio, Cytocom: a cytoscape app to visualize, query and analyse disease comorbidity networks, Bioinformatics 31 (2014) 969–971. [10] Z. Zhang, Semi-parametric regression model for survival data: graphical visualization with r, Annals of translational medicine 4 (2016). [11] M. A. Moni, P. Lio, comor: a software for disease comorbidity risk assessment, Journal of clinical bioinformatics 4 (2014) 1. [12] W. Zhang, T. Ota, V. Shridhar, J. Chien, B. Wu, R. Kuang, Networkbased survival analysis reveals subnetwork signatures for predicting outcomes of ovarian cancer treatment, PLoS computational biology 9 (2013) e1002975. [13] B. Friedenson, Brca1 and brca2 pathways and the risk of cancers other than breast or ovarian, Medscape General Medicine 7 (2005) 60. [14] K. Tecza, J. Pamula-Pilat, Z. Kolosza, N. Radlak, E. Grzybowska, Genetic polymorphisms and gene-dosage effect in ovarian cancer risk and response to paclitaxel/cisplatin chemotherapy, Journal of Experimental & Clinical Cancer Research 34 (2015) 2. [15] W.-C. Shangkuan, H.-C. Lin, Y.-T. Chang, C.-E. Jian, H.-C. Fan, K.-H. Chen, Y.-F. Liu, H.-M. Hsu, H.-L. Chou, C.-T. Yao, et al., Risk analysis of colorectal cancer incidence by gene expression analysis, PeerJ 5 (2017) e3003. [16] H.-M. Hsu, C.-M. Chu, Y.-J. Chang, J.-C. Yu, C.-T. Chen, C.-E. Jian, C.-Y. Lee, Y.-T. Chiang, C.-W. Chang, Y.-T. Chang, Six novel immunoglobulin genes as biomarkers for better prognosis in triple-negative breast cancer by gene co-expression network analysis, Scientific reports 9 (2019) 4484. 24

515 516 517 518 519

520 521 522 523 524

525 526 527 528

529 530 531

532 533 534

535 536

537 538 539

540 541

542 543 544 545

[17] H.-F. Chang, C.-C. Wu, C.-A. Sun, C.-M. Chu, F.-G. Lin, J.-F. Hsieh, C.-H. Hsu, C.-H. Huang, T. Yang, Y.-M. Tsai, et al., Clinical stage and risk of recurrence and mortality: interaction of dna methylation factors in patients with colorectal cancer, Journal of Investigative Medicine 64 (2016) 1200–1207. [18] Y.-T. Chang, C.-T. Yao, S.-L. Su, Y.-C. Chou, C.-M. Chu, C.-S. Huang, H.-J. Terng, H.-L. Chou, T. Wetter, K.-H. Chen, et al., Verification of gene expression profiles for colorectal cancer using 12 internet public microarray datasets, World Journal of Gastroenterology: WJG 20 (2014) 17476. [19] M. A. Hossain, T. A. Asa, M. R. Rahman, M. A. Moni, Networkbased approach to identify key candidate genes and pathways shared by thyroid cancer and chronic kidney disease, Informatics in Medicine Unlocked (2019) 100240. [20] T. Ahmed, M. Begum, Association between gene expression, clinical factors and survival in patients with breast cancer, Journal of Biomedical Analytics 1 (2018) 1–14. [21] L. Arzuman, M. A. Moni, P. Beale, Q. Y. Jun, M. Molloy, J. M. Quinn, F. Huq, Protein expression patterns in ovarian cancer cells associated with monofunctional platinums treatment, bioRxiv (2019) 628958. [22] D. R. Cox, Regression models and life-tables, in: Breakthroughs in statistics, Springer, 1992, pp. 527–541. [23] R. M. Gabriel, G. B. Glavin, Multivariate data analysis in empirical research, The Pavlovian Journal of Biological Science: Official Journal of the Pavlovian 13 (1978) 93–112. [24] G. Heinze, M. Schemper, A solution to the problem of monotone likelihood in cox regression, Biometrics 57 (2001) 114–119. [25] M. R. Rahman, T. Islam, M. A. Al-Mamun, T. Zaman, M. R. Karim, M. A. Moni, The influence of depression on ovarian cancer: Discovering molecular pathways that identify novel biomarkers and therapeutic targets, Informatics in Medicine Unlocked 16 (2019) 100207.

25

546 547 548

549 550 551

552 553 554 555

556 557 558 559

560 561 562

563 564 565

566 567 568

569 570

571 572 573 574 575

576 577

[26] H. Xu, M. A. Moni, P. Li`o, Network regularised cox regression and multiplex network models to predict disease comorbidities and survival of cancer, Computational biology and chemistry 59 (2015) 15–31. [27] K. Tomczak, P. Czerwi´ nska, M. Wiznerowicz, The cancer genome atlas (tcga): an immeasurable source of knowledge, Contemporary oncology 19 (2015) A68. [28] E. Cerami, J. Gao, U. Dogrusoz, B. E. Gross, S. O. Sumer, B. A. Aksoy, A. Jacobsen, C. J. Byrne, M. L. Heuer, E. Larsson, et al., The cbio cancer genomics portal: an open platform for exploring multidimensional cancer genomics data, 2012. [29] A. Hamosh, A. F. Scott, J. S. Amberger, C. A. Bocchini, V. A. McKusick, Online mendelian inheritance in man (omim), a knowledgebase of human genes and genetic disorders, Nucleic acids research 33 (2005) D514–D517. [30] J. Amberger, C. Bocchini, A. Hamosh, A new face and new challenges R for online mendelian inheritance in man (omim ), Human mutation 32 (2011) 564–567. [31] T. Adib, S. Henderson, C. Perrett, D. Hewitt, D. Bourmpoulia, J. Ledermann, C. Boshoff, Predicting biomarkers for ovarian cancer using gene-expression microarrays, British journal of cancer 90 (2004) 686. [32] M. R. McLemore, C. Miaskowski, B. E. Aouizerat, L.-m. Chen, M. J. Dodd, Epidemiologic and genetic factors associated with ovarian cancer, Cancer nursing 32 (2009) 281. [33] C. SJ., Ovarian cance— Cancer, genetics web, 2017. [Online; accessed 05-Aprial-2018]. [34] Which genes are most likely to have a mutation that increases the risk of ovarian cancer?....Target ovarian cancer, https://www.targetovariancancer.org.uk/ information-and-support/hereditary-ovarian-cancer/ which-genes-are-most-likely-have-mutation, 2014. [35] L. Coenegrachts, D. Garcia-Dios, J. Depreeuw, M. Santacana, S. Gatius, M. Zikan, P. Moerman, L. Verbist, D. Lambrechts, X. Matias-Guiu, 26

578 579 580

581 582 583 584

585 586 587 588

589 590 591 592

593 594 595

596 597 598

599 600 601

602 603 604 605

606 607 608 609

et al., Mutation profile and clinical outcome of mixed endometrioidserous endometrial carcinomas are different from that of pure endometrioid or serous carcinomas, Virchows Archiv 466 (2015) 415–422. [36] Z. Li, M. S. Block, R. A. Vierkant, Z. C. Fogarty, S. J. Winham, D. W. Visscher, K. R. Kalli, C. Wang, E. L. Goode, The inflammatory microenvironment in epithelial ovarian cancer: a role for tlr4 and myd88 and related proteins, Tumor Biology 37 (2016) 13279–13286. [37] R. G. Moore, E. K. Hill, T. Horan, N. Yano, K. Kim, S. MacLaughlan, G. Lambert-Messerlian, Y. D. Tseng, J. F. Padbury, M. C. Miller, et al., He4 (wfdc2) gene overexpression promotes ovarian tumor growth, Scientific reports 4 (2014) srep03574. [38] M. A. Cuello, S. Kato, F. Liberona, The impact on high-grade serous ovarian cancer of obesity and lipid metabolism-related gene expression patterns: the underestimated driving force affecting prognosis, Journal of cellular and molecular medicine 22 (2018) 1805–1815. [39] P. G. Moorman, R. T. Palmieri, L. Akushevich, A. Berchuck, J. M. Schildkraut, Ovarian cancer risk factors in african-american and white women, American journal of epidemiology 170 (2009) 598–606. [40] Ovarian cancer stages....American cancer https://www.cancer.org/cancer/ovarian-cancer/ detection-diagnosis-staging/staging.html, 2017.

society,

[41] Tutorial: Survival analysis of tcga patients integrating gene expression (rnaseq) data.....Biostars, https://www.biostars.org/p/153013/ #179081, 2017. [42] M. E. Hossain, A. Khan, M. A. Moni, S. Uddin, Use of electronic health data for disease prediction: A comprehensive literature review, IEEE/ACM transactions on computational biology and bioinformatics (2019). [43] E. Ratner, L. Lu, M. Boeke, R. Barnett, S. Nallur, L. J. Chin, C. Pelletier, R. Blitzblau, R. Tassi, T. Paranjape, et al., A kras-variant in ovarian cancer acts as a genetic marker of cancer risk, Cancer research (2010) 0008–5472. 27

610 611 612 613

614 615 616

617 618 619

620 621 622 623

624 625 626 627

628 629 630 631

632 633

634 635 636 637

638 639 640

[44] D. Szklarczyk, J. H. Morris, H. Cook, M. Kuhn, S. Wyder, M. Simonovic, A. Santos, N. T. Doncheva, A. Roth, P. Bork, et al., The string database in 2017: quality-controlled protein–protein association networks, made broadly accessible, Nucleic acids research (2016) gkw937. [45] B. Hao, Z. Chen, B. Bi, M. Yu, S. Yao, Y. Feng, Y. Yu, L. Pan, D. Di, G. Luo, et al., Role of tlr4 as a prognostic factor for survival in various cancers: a meta-analysis, Oncotarget 9 (2018) 13088. [46] K. Sawada, C. Ohyagi-Hara, T. Kimura, K.-i. Morishige, Integrin inhibitors as a therapeutic agent for ovarian cancer, Journal of oncology 2012 (2012). [47] M. S. Block, R. A. Vierkant, P. F. Rambau, S. J. Winham, P. Wagner, N. Traficante, A. Toloczko, D. G. Tiezzi, F. A. Taran, P. Sinn, et al., Myd88 and tlr4 expression in epithelial ovarian cancer, in: Mayo Clinic Proceedings, volume 93, Elsevier, pp. 307–320. [48] A. G. Kutikhin, A. E. Yuzhalin, A. N. Volkov, A. S. Zhivotovskiy, E. B. Brusina, Correlation between genetic polymorphisms within il-1b and tlr4 genes and cancer risk in a russian population: a case-control study, Tumor Biology 35 (2014) 4821–4830. [49] C. Brignole, D. Marimpietri, D. Di Paolo, P. Perri, F. Morandi, F. Pastorino, A. Zorzoli, G. Pagnan, M. Loi, I. Caffa, et al., Therapeutic targeting of tlr9 inhibits cell growth and induces apoptosis in neuroblastoma, Cancer research 70 (2010) 9816–9826. [50] M. A. Moni, P. Lio, How to build personalized multi-omics comorbidity profiles, Frontiers in cell and developmental biology 3 (2015). [51] X. Yu, Y. Wang, H. Qiu, H. Song, D. Feng, Y. Jiang, S. Deng, H. Meng, J. Geng, Aeg-1 contributes to metastasis in hypoxia-related ovarian cancer by modulating the hif-1alpha/nf-kappab/vegf pathway, BioMed research international 2018 (2018). [52] K. Yoshida, Y. Miki, Role of brca1 and brca2 as regulators of dna repair, transcription, and cell cycle in response to dna damage, Cancer science 95 (2004) 866–871.

28

641 642

643 644 645 646

647 648

649 650 651 652

653 654 655 656

657 658 659

660 661 662 663

664 665 666 667 668

669 670 671

[53] S. J. Ramus, S. A. Gayther, The contribution of brca1 and brca2 to ovarian cancer, Molecular oncology 3 (2009) 138–150. [54] N. Biglia, P. Sgandurra, V. E. Bounous, F. Maggiorotto, E. Piva, E. Pivetta, R. Ponzone, B. Pasini, Ovarian cancer in brca1 and brca2 gene mutation carriers: analysis of prognostic factors and survival, ecancermedicalscience 10 (2016). [55] H. M. Sowter, A. Ashworth, Brca1 and brca2 as ovarian cancer susceptibility genes, Carcinogenesis 26 (2005) 1651–1656. [56] M. L. Gasparri, E. Bardhi, I. Ruscito, A. Papadia, A. A. Farooqi, C. Marchetti, G. Bogani, I. Ceccacci, M. D. Mueller, P. B. Panici, Pi3k/akt/mtor pathway in ovarian cancer treatment: Are we on the right track?, Geburtshilfe und Frauenheilkunde 77 (2017) 1095. [57] M. A. Hossain, T. A. Asa, M. R. Rahman, M. A. Moni, Networkbased approach to identify key candidate genes and pathways shared by thyroid cancer and chronic kidney disease, Informatics in Medicine Unlocked (2019) 100240. [58] Q. Wang, B. Wang, Y.-m. Zhang, W. Wang, The association between cdh1 promoter methylation and patients with ovarian cancer: a systematic meta-analysis, Journal of ovarian research 9 (2016) 23. [59] H. Oka, H. Shiozaki, K. Kobayashi, M. Inoue, H. Tahara, T. Kobayashi, Y. Takatsuka, N. Matsuyoshi, S. Hirano, M. Takeichi, et al., Expression of e-cadherin cell adhesion molecules in human breast cancer tissues and its relationship to metastasis, Cancer research 53 (1993) 1696–1701. [60] C.-M. Chu, C.-J. Chen, D.-C. Chan, H.-S. Wu, Y.-C. Liu, C.-Y. Shen, T.-M. Chang, J.-c. Yu, H.-J. Harn, C.-P. Yu, et al., Cdh1 polymorphisms and haplotypes in sporadic diffuse and intestinal gastric cancer: A case– control study based on direct sequencing analysis, World journal of surgical oncology 12 (2014) 80. [61] Y.-L. Zhang, R.-C. Wang, K. Cheng, B. Z. Ring, L. Su, Roles of rap1 signaling in tumor cell migration and invasion, Cancer biology & medicine 14 (2017) 90.

29

672 673 674

675 676

677 678 679

680 681 682

683 684 685 686 687

688 689 690 691 692

693 694 695 696 697

698 699 700 701

702 703

[62] H. Luo, X. Xu, M. Ye, B. Sheng, X. Zhu, The prognostic value of her2 in ovarian cancer: a meta-analysis of observational studies, PloS one 13 (2018) e0191972. [63] Q. Sheng, J. Liu, The therapeutic potential of targeting the egfr family in epithelial ovarian cancer, British journal of cancer 104 (2011) 1241. [64] S. Granot-Attas, A. Elson, Protein tyrosine phosphatase epsilon activates yes and fyn in neu-induced mammary tumor cells, Experimental cell research 294 (2004) 236–243. [65] C. Coveney, D. Boocock, R. Rees, S. Deen, G. Ball, Data mining of gene arrays for biomarkers of survival in ovarian cancer, Microarrays 4 (2015) 324–338. [66] S. Bellone, R. Tassi, M. Betti, D. English, E. Cocco, S. Gasparrini, I. Bortolomai, J. Black, P. Todeschini, C. Romani, et al., Mammaglobin b (scgb2a1) is a novel tumour antigen highly differentially expressed in all major histological types of ovarian cancer: implications for ovarian cancer immunotherapy, British journal of cancer 109 (2013) 462. [67] K. Fischer, A.-C. von Br¨ unneck, D. Hornung, C. Denkert, C. Ufer, H. Schiebel, H. Kuhn, A. Borchert, Differential expression of secretoglobins in normal ovary and in ovarian carcinoma–overexpression of mammaglobin-1 is linked to tumor progression, Archives of biochemistry and biophysics 547 (2014) 27–36. [68] M. Zafrakas, B. Petschke, A. Donner, F. Fritzsche, G. Kristiansen, R. Kn¨ uchel, E. Dahl, Expression analysis of mammaglobin a (scgb2a2) and lipophilin b (scgb1d2) in more than 300 human tumors and matching normal tissues reveals their co-expression in gynecologic malignancies, BMC cancer 6 (2006) 88. [69] K. Tanahashi, K. Sugiura, M. Kono, H. Takama, N. Hamajima, M. Akiyama, Highly prevalent liph founder mutations causing autosomal recessive woolly hair/hypotrichosis in japan and the genotype/phenotype correlations, PloS one 9 (2014) e89261. [70] T. Mizushima, H. Miyamoto, The role of androgen receptor signaling in ovarian cancer, Cells 8 (2019) 176. 30

Highlights: 

TLR4, BSCL2, CDH1, ERBB2, and SCGB2A1 gene expression affects patient survival.



5 genes acted in concert to affect OC patient survival, so may be drug targets.



TLR4, BSCL2, CDH1, ERBB2, SCGB2A1 and BRCA2 most affected clinical outcomes.



After taking clinical factors into account, ethnicity affected patient survival.



7 of our identified genes linked to pathways associated with OC progression.