Study population and treatment protocol

The CC and MGC were collected at the Fertility Clinic, Copenhagen University Hospital, Rigshospitalet, Denmark between February 2009 and June 2010 as a part of a prospective, randomised study (0 IU hCG; 50 IU hCG, 100 IU hCG and 150 IU hCG daily during controlled ovarian stimulation (COS)) [17]. The Danish National Committee on Biomedical Research Ethics (HB-2008-146) and The Danish Medicines Agency (2612–3928, EudraCT number 2008-008355-42) approved the study. All participants provided both written and verbal informed consent. The consent was documented with participant’s signature and the consent form was approved by the Ethics Committes.

A total of 60 women, scheduled for IVF, aged 25–37 years and with BMI >18 and < 30 kg/m², regular menstrual cycle, early follicular phase serum FSH levels within 1–12 IU/l and early follicular phase total antral follicle (2–10mm) count ≥ 6 donated CC and MGC. The patients underwent controlled ovarian stimulation in a GnRH agonist protocol as described in [17]. Oocyte retrieval was performed 36 h (±2 h) after hCG administration. Single embryo transfer was performed Day 3 and 16 pregnancies resulted in live births. All clinical (week 7) and ongoing pregnancies (week 10–12) were confirmed with an ultrasound scan and followed up to delivery. The sixteen CC and corresponding MGC representing the competent oocytes and follicles giving rise to the babies were included in the study. The competent oocytes were evenly distributed in the 4 study groups [17]. To represent incompetent oocytes, another 16 CC and corresponding MGC were chosen among the samples with follicle cells connected to the oocytes that gave rise to an embryo that failed to implant. These samples were chosen to match the competent oocytes in relation to COS treatment, female age and embryo morphology.

For demographic data t-test and chi-squared tests were used and p<0.05 considered significant.

Cumulus and granulosa cell isolation

Follicular fluids from all follicles were collected individually in tubes (Nunc Centrifuge 11 ml, Nunc, Denmark). When an oocyte was obtained, the number of the follicle was recorded and the oocytes were numbered consecutively. Thus, each oocyte could be linked to an individual follicle and hence the correspondent MGC. Immediately after isolation of the COC, the MGC were isolated from the follicular fluid by centrifugation at 400 g in 8 min at room temperature. The pellet was re-suspended in 1 ml HEPES-buffered solution (Sydney IVF Follicle Flush Buffer, COOK Medical, Australia) and placed on a density gradient column (Histopaque 1077, Sigma-Aldirch, St-Louis, USA) and centrifuged at 400 g for 30 minutes. After centrifugation, the inter-phase containing the MGC was aspirated and the MGC washed twice by mixing with1 ml of the buffer solution followed by centrifugation at 400 g for 8 minutes at room temperature. After the second wash, pellet was re-suspended in 500 μl in a cryo-vial (CryoVial, 2 ml Greiner bio-one, Germany) and stored at -80°C until RNA extraction.

Within 1 hour after the oocyte retrieval a part of the cumulus mass was mechanically removed (18G needles) from the COC and transferred to a cryo-vial (CryoVial, 2 ml Greiner bio-one, Germany) with as little fluid as possible and stored at -80°C until RNA extraction.

Embryo Evaluation

IVF insemination was done with standard IVF procedure at 3 h (±1 h) after oocyte retrieval. All oocytes were followed individually and zygote and embryo quality assessed at 20 h (±1 h), 28 h (±1 h), 44 h (±1 h) and 68 h (±1 h) after oocyte retrieval. The embryo quality evaluation consisted of assessment of cell number and three parameters of embryo morphology: degree of fragmentation, cell stage specific blastomere uniformity and multinucleation (ALPHA and ESHRE consensus parameters, 2011).

RNA purification

RNA from individual samples was purified using Arcturus PicoPure RNA Isolation kit (Applied Biosystems®, Naerum, Denmark) according to manufactures instructions but with small modifications. In short: 100 μl extraction buffer was used for initial lysis of each sample, except for samples were sample volume excided 70 μl. For these samples 2 x sample volume of extraction buffer was added. Following lysis and centrifugation, the supernatant was transferred to a 2 ml Eppendorf tube from the kit (Arcturus PicoPure RNA Isolation kit, Applied Biosystems®, Naerum, Denmark) where it was mixed 1:1 with 70% ethanol. For some MGC samples with large volume, several transfers and spins were needed to collect all the RNA in the spin tubes. In the end, the RNA was extracted in 12 μl nuclease free water, giving a ~10 μl end volume. The integrity of the RNA was measured using bioanalyzer pico chip (RNA 6000 Pico Assay for 2100 bioanalyzer, Agilent Technologies) and only samples with a RNA integrity number (RIN) above 7 as analyzed by the 2100 bioanalyzer (Agilent Technologies, Glostrup, Denmark) were included in the analysis.

Microarray analysis

RNA was amplified and labeled using a pico amplification kit according to manufactures instructions. In short, approximately half of the total RNA from each sample (5 μl of 10 μl) was amplified using the Ovation Pico WTA v.2 RNA Amplification System from NuGEN® Inc. (NuGEN®, San Carlos, CA, USA) and biotin labeling was performed with the Encore Biotin Module (NuGEN®). The labeled samples were hybridized to the Human Gene 1.0 ST GeneChip array (Affymetrix, Santa Clara, CA, USA). The arrays were washed and stained with phycoerytrin conjugated streptavidin (SAPE) using the Affymetrix Fluidics Station® 450, and the arrays were scanned in the Affymetrix GeneArray® 3000 7G scanner to generate fluorescent images, as described in the Affymetrix GeneChip® protocol. Cell intensity files (CEL files) were generated in the GeneChip® Command Console® Software (AGCC) (Affymetrix, Santa Clara, CA, USA).

All samples are MIAME compliant and were handled according to SOP in the Microarray Center. A total of 27 CC and 19 MGC samples were amplified and hybridized to arrays. The 46 samples were submitted to ArrayExpress at EMBL using MIAMExpress. The experiment accession number is E-MTAB-4012.

Formulation of Classifiers

Classifiers were constructed based on expression matrices obtained by Plier normalization [18]. Following the normalization, the data was initially filtered to include only probe sets with functional annotation (both protein coding and non-protein coding); reducing the number of probe sets included in the analysis to 28,054. The formulation and performance assessment of the classifier models was performed on this data set, with and without, an additional pre-filtering step to exclude probe sets expressed below back-ground level (average expression above 10). Back-ground filtering resulted in the inclusion of 17,279 probe sets in the CC data set and 15,146 probe sets in the GC data set, respectively.

Based on these data sets, three different classification algorithms (diagonal linear discriminant analysis (dLDA) [19,20], k-nearest neighbors (K-NN) [20,21] and support vector machines (SVM) with linear kernel [22,23] were generated and tested. In order to formulate and select the best model for classification of LB, the misclassification rate for each classifier was estimated using leave-one-out cross-validation (LOOCV), during which we applied uni-variate t-tests with a grit of p-values (p < 0.001, p < 0.005, p < 0.01, p < 0.05) for selection of the optimal number of probe sets to include in each model. The overall accuracy, as well as PPV and NPV were derived for both the cumulus and granulosa data separately. The output of these analyses is presented for the optimal classification models and model parameters (S2 and S3 Tables). The overall best performing SVM classifier, based on the back-ground filtered data set with a p-value cut-off of 0.01, resulting in a gene signature of 34 probe sets, representing 30 unique genes, was selected as the final model for further validation (Table 1, Cumulus and granulosa Classifiers: PLIER annot and PLIER annot above bg).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. Performance of cumulus and granulosa cell support vector machine (SVM).

https://doi.org/10.1371/journal.pone.0153562.t001

Relative weight of the genes in the CC signature

Following the cross validation procedure, we calculated the cross-validation support for each gene showing the percentage of the cross-validation training sets in which each particular gene was selected. The cross-validation support shows the percentage of the cross-validation training sets in which each particular gene was selected, giving the relative strength (weight) of each gene in the signature.

Predictive probability of LB classification

For the optimal CC classifier, the predictive probability (estimate of the probability of getting a particular class label) of each samples, was derived from the trained SVM model. This was achieved by transforming the SMV classifier into a probabilistic classifier by calculating the predictive probability of a sample being one or the other type using logit estimates as implemented in the e1071 R library [24]. The predictive probabilities of the final and optimal classifier are visualized in a scatter plot of the function, p(LB), live birth, by plotting the predictive probability on the y-axis for each consecutive sample shown on the x-axis.

Statistical significance of the Error rate

A permutation test was performed in order to determine if the cross-validated misclassification rate of the final CC classification model was lower than expected by chance (Tusher et al. 2001; Simon et al. 2007). In 1000 random permutations of the class label, the entire cross-validation was repeated for classifying the random grouping of the samples. The proportion of the 1000 random permutations that resulted in a smaller or similar cross-validation misclassification rate as obtained with the real data determine the permutation p-value. The statistical significance of the error rate was determined both for the linear SVM, 1-NN, 3-NN and dLDA classifier.

Gene function enrichment and Biological Networks

The cumulus gene signature was subjected to Ingenuity pathway analysis using IPA Ingenuity software Ingenuity®Systems (www.ingenuity.com). The list of the signature probe sets with p value and fold change was imported into Ingenuity. By application of the Ingenuity Pathway Knowledge Base each gene identifier was overlaid onto a global molecular network linking their functionality to the function of other genes and biological enriched functionality and mechanistic networks of these genes were then generated based on their connectivity and enrichment statistics. A network score was calculated based on the hyper-geometric distribution and calculated with the right-tailed Fisher’s exact test. The score is the negative log of this p-value. The score takes into account, the number of network eligible molecules in the network and its size, as well as the total number of network eligible molecules analyzed and the total number of molecules in the knowledge base that could potentially be included in the network. The score represents the chance of getting a network containing at least the same number of network eligible molecules by chance when randomly picking the number of genes that can be in networks from the knowledge base. The top molecular and cellular functions and networks are shown.

Downstream effect analysis was performed to predict the effect of directional gene expression change on biological processes based on the expected causal effects derived from the literature. The predicted effect is based on a value calculated by the IPA z-score algorithm. The z-score predicts the direction of change for the function. An absolute z-score of ≥ 2 is considered significant. A biological function is predicted to be Increased if the z-score is ≥ 2 and decreased if the z-score ≤ -2.

External validation data set

The optimal CC classifier-signature of 34 probe sets (30 unique gene) was validated on an external data set generated on the Agilent platform using the Whole Human Genome Microarray 4x44K microarray [9]. The full data set comprising three subseries of CC microarray expression profiles was downloaded from the GEO repository at NCBI. The three samples sets used were part 1: GEO accession GSE37110, including 10 embryos of poor quality (EP) and 11 blastocyst (B); part 2: GEO accession GSE37117, including 18 EP and 18 blastocyst B, and part3: GEO accession GSE37116, including 13 EP and 11 B.

The microarray expression profiles were generated in a two-color array with dye-swap of replicate samples to avoid systematic noise in the data related to the choice of dye (cye-5 or cye-3). We downloaded the preprocessed data of normalized unlogged expression values representing one-color channel intensities. The data was used for validation of the gene signature after applying a log2 transformation. Gene symbols were used to select the probes on the 4x44K array that overlap with the unique genes in the 30-gene signature. Due to redundant representation of some of the genes on the Agilent array, the gene signature was represented by 68 redundant mRNA transcripts (representing 25 unique genes) on the Agilent 4x44 platform. For assessing the performance of this gene signature, represented by 25 unique genes, the probe with the intermediate expression value for each of the genes with redundant probes was used to define the 25-gene signature in the external data sets.

Performance of signature on external data

The external data sub sets (GSE37110, GSE37116 and GSE37117) were tested one by one. For each sub set, a test set data object was constructed and the overlapping gene signature was tested using SVM classification with linear kernels. The predictive probabilities were derived and visualized along with the predictive probability of the training set samples.

Receiver operating characteristic (ROC) curves

ROC curve characteristics analysis was performed on all data set (training and validation data) with the end point LB or NP in the training data, and blastocyst or embryo of poor quality used in the validation data set. A ROC curve was generated for each partition of the validation data set. The ROC curve was created by plotting the true positive rate (TPR, or sensitivity) versus the false positive rate (FPR, or 1—specificity). The area under the curve (AUC) indicates the degree of predictive ability of the gene expression, where 0.5 is random chance and 1.0 is perfect predictive ability. The ROC and AUC were produced using the functionality of the R library ROCR.

Clinical presentation of patients in the study

A total of 46 expression profiles were generated, encompassing 27 CC samples and 19 MGC samples, of which 12 CC sample and 9 MGC samples, respectively, were obtained from follicles where the oocyte after implantation gave rise to a live birth of a healthy baby (LB) and 15 CC and 10 MGC samples, respectively, which did not lead to pregnancy (NP). The expression profiles of the CC and MGC samples were analyzed separately. The basic demographic parameters for the women which CC and MGC were analyzed are presented in S1 Table. No differences were present between the two groups.

Construction of classifier for oocyte viability

We constructed a classifier to discriminate between patients where IVF treatment lead to live birth compared with no pregnancy, based on expression profiles on CC or MGC from the pre-ovulatory follicles.

We examined different methods for the construction of an optimal classifier of oocyte competence and tested several cut-off values of the test statistics with or with-out fold change and expression above back-ground filtering. Three different classifier models (diagonal linear discriminant analysis (dLDA), k-nearest neighbors and SVM with linear kernel function were generated and tested by LOOCV. Overall accuracy, as well as PPV and NPV were derived for both the cumulus and granulosa data separately, as shown in Table 1, S2 Table and S3 Table. A grit of p-values (p < 0.001, p < 0.005, p < 0.01, p < 0.05) was tested on all algorithms with LOOCV and the p-value cut-off (p < 0.01), which lead to the optimal performance was used to generate the gene signature of 34 probe sets, representing 30 unique genes which were used in the final classifier (S4 Table).

The overall best performance of the classifier models applied on the CC was achieved with PLIER normalized unlogged expression values by the linear SVM model and 3-NN and resulted in a LOOCV accuracy of 0.81 and 0.85, respectively (Table 1 and S2 Table). On the contrary, the dLDA algorithm showed low cross validation accuracy of 57% (S2 Table).

The performance range of the classifier models applied on PLIER normalized MGC with a p-value cut-off of 0.01 was 11% - 37% for dLDA, 3-NN (k = 3) linear SVM (S3 Table). Evaluation of the signature obtained using RMA normalized data showed slightly better results, although all analyses resulted in overall high error level and poor performance with no signature showing significant discriminative power to distinguish LB samples from NP (data not shown), showing that the low performance is probably independent of normalization strategy.

The CC data show high levels of variance between the samples within each group and relatively small overall expression differences between the LB and NP groups. To increase the power of the analysis we applied a filter to firstly exclude un-annotated probe sets and further exclude genes with an average expression below back-ground level (average expression below 10 across all CC samples) [25]. Training of the linear SVM classifier after the annotation filtering step, resulted in a signature of 82 probe sets with unambiguous annotation, 60 of which were annotated with gene symbol (protein coding) and 22 were annotated as other non-coding RNAs, hence excluding probe sets with unclear function or lacking annotation.

Rerunning the training of the of the linear SVM after filtering for genes expressed above background, the overall accuracy of the SVM with radial kernel was unchanged (81%), when using a signature representing 30 unique genes (Table 1). The full list of the 30 gene signature with p-value cut-off and fold-change is shown in S4 Table.

ROC curve analysis to distinguish between true positive rate and false positive rate of the final linear SVM classifier (annotated probes filtered to exclude probes expressed below background level) resulted in an AUC value of 0.86 indicating a good performance of the CC classifier (Fig 1).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 1. ROC curve of the training set.

Receiver operating characteristic (ROC) curve for the binary classifier built to distinguish between live birth (LB) and no pregnancy (NP). The curve shows the true positive rate versus false positive rate, i.e. the tradeoff between sensitivity and specificity. The area under the curve (AUC), which captures the ability of the classifier to correctly group the patients with follicular adenoma and those with follicular carcinoma, is equal to 0.86. A perfect classifier will have an AUC of 1.0, whereas an AUC value of 0.5 indicates that the classification is random.

https://doi.org/10.1371/journal.pone.0153562.g001

We furthermore assessed the strength of each of the 30 genes in the signature by the cross-validation support, which shows the percentage of the cross-validation training sets in which each particular gene was selected (S4 Table). Analysis of the cross-validation support revealed that 12 genes had CV support of 100 percent; and more than half of the genes in the CC signature were included in 89 percent or more of the cross validation runs, hence receiving high weight in the classifier (see S4 Table). Visualization of the 30 gene signature by hierarchical clustering showed that samples with low expression of the top 100 percent CV support genes were the samples, which received low predictive probabilities of LB (Fig 2, S2 Fig).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. Predictive probability of cumulus training set samples and validation set samples.

The predictive probability of the 30-gene signature is shown for the training set and after translation to the validation data set, part 1, 2 and 3, respectively. Each dot represents a sample and the color indicates the true (blinded) class. If a sample has a predictive value above 0.5 (p(LB) > 0.5), it is classified as predictive of leading to Live birth (LB) in the training set or reaching blastocyst (B) stage in the validation sub sets, validation-part1, validation-part2 and validation-part3, respectively. Samples with p(LB) below 0.5 are classified as predictive of no pregnancy (NP) in the training set and embryos of poor quality (EP) in the validation set. Samples which received an erroneously prediction according to their true class are indicated with a black circle and sample name.

https://doi.org/10.1371/journal.pone.0153562.g002

In order to assess whether the misclassification rate was lower than what would be expected by chance we performed 1000 random permutation of the samples. For each permutation the entire LOOCV process was repeated and the proportion of the 1000 random permutations resulting in smaller misclassification rates determined the statistical significance of the classifiers error rate. The significance of the misclassification rate after 1000 permutations was 0.284, 0.002, 0.009, for the diagonal linear discriminant analysis classifier, the 3-nearest neighbors’ classifier and the linear SVM classifier, respectively. All but the d-LDA classifier was statistically significant compared to the random setting.

Validation of the cumulus gene signature on cross-platform external data.

We tested the cumulus classifier on an external data set published by Feuerstein et al, 2012 downloaded from the public expression profile repository, gene expression omnibus (GEO) at NCBI [9]. The data set consisted of expression profiles of CCs originating from oocytes of IVF treated patients, which developed into blastocysts after 5/6 days of in-vitro culture or oocytes which were arrested at the embryo state. Three sub-sets of data were available as part of the study (part 1: GSE37110, part 2: GSE37116 and part 3: GSE37117) consisting of 21, 36 and 24 samples, respectively. Since these data were analyzed on the whole human genome oligo microarray 4x44K Agilent platform, the data sets were not directly comparable to the Affymetrix format. In order to test the CC signature of 30 unique genes, the gene signature was translated to the probe id of the Agilent platform using gene symbols, which resulted in a signature of 25 unique genes. The performance of this 25 gene signature was evaluated by linear SVM on pre-processed expression values downloaded from the GEO data base. Testing of the gene signature resulted in varying performance across the three parts of the data set, with part 3 data showing superior accuracy of 0.88 compared to 0.6 and 0.75 respectively for part 1 and part 2 (Table 2).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 2. Performance of the gene signature obtained during the formulation and training of the cumulus cell support vector machine on the external validation data set.

https://doi.org/10.1371/journal.pone.0153562.t002

The SVM predictions were translated into predictive probabilities and visualized alongside the predictive probabilities of the training data for comparison (Fig 2). In agreement with these results, ROC analysis resulted in AUC of 0.57, 0.65 and 0.88, respectively for the three sub data sets (S1 Fig).

Ingenuity pathway analysis.

Despite the fact that gene classification signatures often can be difficult to interpret with respect to the biology of the underlying disease; and missing biologic interpretation does not warrant poor clinical usefulness of well-established biomarkers [26], the probe sets list including p-value, false discovery rate (FDR) q-value, and fold change of the comparison of LB versus NP, was imported into the Ingenuity Pathway Analysis (IPA) software to investigate biological functionality. Based on the relatively short gene signature list, no upstream activator molecules were detected. However, downstream effect analysis to predict the effect of directional gene expression change on biological processes based on the expected causal effects derived from the literature, resulted in two biological functions with significant activation z-scores. Down-regulation of 10 genes in the LB group collectively pointed to the activation of apoptosis with a z-score of 3.032 and an enrichment p-value of 0.018 (Fig 3). Four related bio functions in the categories of Cell Death and Survival, and Organismal Injury and Abnormalities were predicted to be activated in the LB group, although with non-significant z-scores below 2 (S5 Table). Six of the 10 genes were furthermore overlapping with a set of eight genes, which were predictive of decreased activity of cell migration (Fig 3, S5 Table). Additional four functional categories showed z-scores indicative of decreased activity although with low significance. These were microtubule dynamics, organization of cytoplasm, quantity of neurons and growth of neurites, all part of the biological categories Cellular Assembly and Organization, and Tissue Morphology (S5 Table).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. Predicted activation states.

Ingenuity downstream effect analysis to predict the effect of directional gene expression resulted in two biological functions with significant activation z-scores indicative of predicted activation of apoptosis and decreased activity of cell migration. The mechanistic network of the implicated genes is shown along with the predicted relationship indicated by the color of the edges.

https://doi.org/10.1371/journal.pone.0153562.g003

Mechanistic networks were generated based on sets of connected upstream regulators that in conjunction would result in the gene expression changes observed in the dataset. Two mechanistic networks were generated with network scores above 6. The first network was enriched in functions related to Cellular Movement, Nervous System Development and Function, Cellular Growth and Proliferation with a network score of 24 (Fig 4A). Eleven of the 30 genes in the classifier gene signature are represented in the mechanistic network, all of which are down regulated in the LB group compared to NP. The second mechanistic network represents genes involved in Cell Death and Survival, Cell-To-Cell Signaling and Interaction, Hematological System Development and Function. Although only four of the 35 genes in the network are represented in the gene signature, they are part of the functional network involved in cell death, which is in line with the increased predicted activity of apoptosis (Fig 4B).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. Mechanistic network.

Mechanistic network based on the 30 probe set signature used to predict live birth in cumulus cells. (A) Top predicted functions represented by the network are Cellular Movement, Nervous System Development and Function, Cellular Growth and Proliferation. (B). Top predicted functions represented by the network are Cell Death and Survival, Cell-To-Cell Signaling and Interaction, Hematological System Development and Function.

https://doi.org/10.1371/journal.pone.0153562.g004