Ethics statement

Ethics committee approval for conducting the study, use of consent forms and scientific evaluation of the data were obtained from Heidelberg University Hospital’s review board. The approval was granted by the named board (S-502/2013). All patients or their legal guardian signed a written informed consent.

Tissue collection and DNA extraction

All next generation sequencing (NGS) and DNA-methylation array data were generated and obtained within the scope of the INdividualized Therapy FOr Relapsed Malignancies in Childhood (INFORM) project [22]. Here is a brief description of the INFORM methodology including OS, described in detail in [22].Only patients with relapsed or refractory osteosarcoma between the ages 2–30 were included in the study. Based on histological assessment the content of tumor in the tumor material must be higher than 40%. Fresh frozen samples of tumor material and of its matching control were quality assessed before their molecular analysis.

DNA sequencing, quality control and read alignment

For the whole-exome sequencing (WES) the paired-end libraries were prepared with the Agilent SureSelectXT Human V5 kit. Low-coverage whole-genome sequencing (lcWGS) was prepared with the exclusion of the enrichment step in SureSelectXT Human V5 kit. RNA-Seq paired-end 100 bp from poly(A)+ RNA was obtained with Illumina TruSeq RNA Kit v2 on an Illumina HiSeq2500 only tumor reverse DNA library. The mean coverage WES was 165-fold for the tumor and 117-fold for the control, WGS coverage was 3.4-fold for the tumor and 2.5-fold for the control. Average number of RNA sequencing reads was 135 million reads.

Quality control, trimming and filtering of FASTQ-files was done with PRINSEQ (ver.0.20.4) [23].Reads were aligned in SAM format with the BWA-MEM algorithm (bwa ver. 0.7.16a) [24]. A full 1000genomes phase2 reference genome sequence hs37d5 was used for alignment. Sam-files were converted to bam-files, sorted by coordinates and PCR duplicates were removed with the MarkDuplicate subprogram of the Picard software suite (ver. 2.16.0) [25].Bam-files were indexed with samtools [26]. For the final indel realignment, the Genome Analysis Toolkit module (ver.3.8) (GATK) was used [27].

Genomic scars’ method

Based on the previous studies, NGS-based assay myChoiceHRD (Myriad Genetics) was developed to estimate the degree of so called genomic scar HRD score, i.e. a score combined of counts of specific genomic events such as loss of heterozigosity (LOH), telomeric allelic imbalance (TAI), large-scale transitions (LST) and BRCA1/2 sequences. HRD-LOH score is represented by the number of LOH regions longer than 15 Mb [15].HRD-TAI is a number of allelic imbalance regions extending to subtelomeres not crossing the centromere and longer than 11 Mb [9,13]. The HRD-LST score was defined as the number of break points between regions longer than 10 Mb after filtering out regions shorter than 3 Mb [14]. Furthermore, the HRD-LST score cutoff between intact and deficient samples was modified based on the samples’ ploidy [28]. myChoiceHRD represents the sum of these three scores, i.e. HRD-LOH, HRD-TAI and HRD-LST. In order to estimate the HRD score in an initial set of 22 OS INFORM samples, genomic scarring events were counted manually, relying on human visualization via Nexus CN software (BioDiscovery, Inc) and the myChoice score was calculated by summing up HRD-LOH, HRD-TAI and HRD-LST events. XY chromosomes’ genomic scarring events were omitted.

Automation of HRD score assessment with NexusCN and Control-FREEC/HRDtools

In order to automate and remove the human interpretation of genomic scarring events which might cause significant fluctuations in HRD scoring, we employed the NexusCN feature Percent of Genome Changed (PGC) which shows what percent of the genome has aberrations. The filter data options were set as follows: remove gains smaller than 10000 kb, remove losses smaller than 10000 kb, remove all allelic imbalance calls, remove all LOH calls. Since we use only lcWGS, TAI and LOH calls were not considered.

To calculate PGC with open source software we used Control-FREEC (version 11.4.3) and HRDtools (version 0.0.1) software [29]. Control-FREEC was developed as a tool for detection and enumeration of copy-number alterations (CNA), TAI and LOH using NGS data [30]. Briefly, from tumor-control bam-files we generated pileup-files using samtools mpileup [26]. Then, we ran Control-FREEC with pileup-files as an input. We post-process Control-FREEC output to make it compatible with HRDtools. We modified HRDtools so that it can closely follow myChoice HRD genomic size filtering (LOH = 15 Mb, LST = 10 Mb, TAI = 11 Mb).Modified HRDtools provided in the output the width of CNA, TAI and LOH. Dividing summed up width by the length of the haploid genome (3,234,834,689—bases GRCh37.p13) we obtained a PGC score.XY chromosomes were not considered for the overall PGC score.

DNA-methylation arrays

The tumor DNA-methylation profile was assessed with Illumina HumanMethylation450 and Epic BeadChip (482,421/850,000 CpG sites, respectfully) according to the manufacturer’s instructions. There were 43 DNA methylation arrays of OS generated within the INFORM project. Analysis and processing of 450k and Epic arrays were done with an integrated analysis pipeline, R package ChAMP, which includes filtering low quality probes based on detection P-values, chromosomal location and presence of single nucleotide polymorphisms in the probe sequence. For filtering and processing the default values were kept [31]. In order to utilize both 450k and Epic arrays we used the function convertArray from R package minfi [32]. This function converts EPIC into 450k arrays by dropping probes that differ between two array types (Modified champ.load function from package ChAMP with incorporated convertArray method in Supplement). After obtaining beta-values, defined as the ratio of a probe fluorescence intensity over the overall intensity, beta-values were further normalized with Beta MIxture Quantile dilation (BMIQ) adjusting the beta-values of type II design probes into a statistical distribution characteristic of type I probes. BMIQ was suggested as the preferred methods for case—control, tumor—normal or, in our case, BRCAness negative—positive dataset processing. BMIQ was shown to be more effective at removing unwanted technical bias than other normalization methods [33]. Methylation probes were annotated with HumanMethylation450 v1.2 manifest (Illumina).

Setting the threshold for BRCAness positive samples

A tumor with a myChoice score above 42 was considered HRD/BRCAness-positive [9]. NexusCN PGC score above 32, FREEC-HRDtools PGC score above 28 were considered as BRCAness-positive after matching PGC scores to myChoice scores. Throughout this study we used the NexusCN PGC score, as the original calculation of the myChoice score for OS samples was done with NexusCN assistance. 43 OS INFORM samples were grouped to two classes of 26 and 17 samples, i.e. BRCAness-positive and BRCAness-negative group, respectively, based on the threshold of a NexusCN PGC score of 32 which corresponds to a myChoiceHRD score of 42.

Feature selection and testing class structure

First, variances of each probe were calculated within the BRCA positive (var_pos) and negative (var_neg) training group. var_all was computed for each sample regardless of the grouping. It can be represented as Δ-variances = var_all–(var_pos + var_neg) [34]. Probes were ranked in descending order according to Δ-variances and the top 2,000 probes (features) were selected for further random forest (RF) training and validation.

Since the number of predictors was considerably larger than the number of samples and to avoid importance bias caused by multiple highly correlated probes [35], only one probe was retained from probes correlated above 0.8. As a result, only 2.7% (54 probes) of the original probe (feature) set of 2000 probes remained. We used 41 sample set and uncorrelated predictors further in the training and testing of the RF. Correlated probes were included in the biological functional analysis but excluded from the RF training and validation stages.

To test whether we found a real class structure in the data, the corresponding null distribution was estimated by permuting the labels in the data and by enumerating the number of positive Δ-variances. Our grouping based on PGC and HRD scoring had 3,685 positive Δ-variances. The labels were permuted 1000 times and every time the number of probes with positive Δ-variances was counted. In 30 out of 1000 runs only a single probe had positive Δ-variance, the rest of the runs had zero positive Δ-variance probes. This result supports our approach to OS sample grouping and our assumption that only probes with positive Δ-variances are informative.

Random forest

To identify the predictive set of probes underlying BRCAness we utilized a machine learning approach, random forests, using the R randomForest package (version 4.6–14) in R [36–38]. For the most part, we followed a predictive modelling approach and the R script described in Valletta et al. with parameters modified for our dataset [39]. For RF training our unbalanced set of samples was stratified for class balancing. As input variables for our models we used the intensity of methylation probes. The response variable, i.e. the outcome to be predicted, was sample assignment either to BRCAness-positive (’Yes’) or to BRCAness-negative group (’No’) as it was defined by the NexusCN PGC score. The number of trees was set to 10000, 32 variables were tried at each node and the minimum node size was set to 3. The tree depth limit was not set. Uncorrelated predictors were run through a predictive modelling pipeline, where an additional statistical significance of selected features was computed with a modified mProbes/xRF feature selection method [39–41]. However, all uncorrelated probes showed high predictive power by having high statistical significance after introducing predictors with randomly shuffled values. All uncorrelated predictors were used in RF for training and testing.

Accuracy evaluation

To assess model predictive accuracy, we used RF internal error estimation, i.e. out-of-bag (OOB) error and the receiver operating characteristic (ROC) curve [37].The feature selection algorithm was used internally in a five-fold cross-validation loop. The OOB error confusion matrix generated by RF has been used for accuracy estimation.The area under the ROC curve (AUC) was computed for each fold. The average AUC was computed with combined TP/FP events generated in each fold [39]. The 5-fold CV RF models were not kept. The final RF model was trained on the entire dataset of 41 samples. Accuracy was reported as (TP+TN) / (TP + TN + FP + FN) x 100. Combined accuracy was computed when TP, FP, TN and FN events obtained in each fold of CV were combined together, e.g. if TP1 is TP of fold 1 CV, then TP1 + TP2 + TP3 + TP4 + TP5 = TP_comb. The average (combined) AUC was computed with TP_comb and FP_comb events. The heatmaps and dendrograms were generated with unsupervised hierarchical clustering using Euclidean distance and Ward.D2 clustering method.All plots and statistical tests were performed with statistical software R [38].

RNA-Seq

For alignment and generation of bam-files, a splice-aware RNA-Seq aligner, STAR (ver. 2.7.1a), has been used with default settings [42]. The number of reads mapping to genes (features) from bam-files were calculated with htseq-count which is a part of the HTSeq framework (ver. 0.9.1) [43].Non-unique reads, secondary, supplementary alignments and reads overlapping more than one feature were discarded. Gene reference GTF file was obtained from UCSC/hg19 annotation archive. Read counts were transformed using the DESeq2 framework with regularized logarithm (rlog). rlog transformation generates data on a log2 scale normalized to library size. DESeq2 was used to obtain set of DE genes [44]. Normally, several biological replicates are required for DESeq to identify genes differentially expressed above biological noise. Here, we treated BRCAness-positive and BRCAness-negative OS samples similarly to biological replicates for the purpose of capturing profoundly expressed genes between two classes. After all, most of the changes in gene expression between these two classes can be explained by the genetic heterogeneity of samples rather than random biological noise. Thus, BRCAness-positive and BRCAness-negative samples were grouped conditional on the presence or the absence of BRCAness characteristics. Gene differential expression was considered statistically significant with fdr < 0.05.

Gene set enrichment analysis

A pre-ranked list of all genes based on logFC sign and 1/p-val was analyzed by the gene set enrichment analysis tool (GSEA ver. 4.0.2) [45,46]. The enrichment was evaluated using Molecular Signatures database (MSigDB) sets, namely Hallmark set, KEGG (MSigDB ver.7.0) [47–49]. Enrichment was considered significant if the gene set family wise error rate (FWER) was less than 0.05. GSEA particular settings were set to classic enrichment statistics, no normalization mode, mean of probes. Differentially expressed upregulated and downregulated sets of genes with fdr < 0.05 were uploaded to the gene set enrichment web-tool Enrichr [50]. Enrichment was considered significant if the gene expression signature had fdr < 0.05.

HRD score evaluation with NexusCN and Control-FREEC/HRDtools

We simplified the computation of HRD scores with PGC in order to separate samples to BRCAness-positive and BRCAness-negative groups for further RF training and validation. We observed a strong correlation of PGC with canonical HRD (LST + TAI + LOH). Spearman rho correlation between the myChoice HRD score and the NexusCN PGC, and between the myChoice HRD score and the FREEC-HRDtools PGC were 0.92 (p-val = 1.3e-09) and 0.86 (p-val = 3.2e-07), respectively. Spearman rho correlation between the FREEC-HRDtools PGC and the Nexus CN PGC was 0.91 (p < 2.2e-16) (Fig A in S1 File). Since the correlation between the myChoice HRD score and the NexusCN PGC was higher and did not require human intervention, we continued with the NexusCN PGC as a BRCAness score replacing the myChoice HRD score except validating independent test set where we used the FREEC-HRDtools PGC, since the NexusCN software was no longer available. Overall, we showed that the PGC approach is suitable for separation of samples to two groups. However, using PGC as an alternative to the canonical HRD score requires testing it on a large number of tumor samples, such as ovarian cancer which was initially used to develop and test HRD scoring.

Analysis and filtering of sample set

First, we analyzed a full dataset of 43 samples where methylation probes were ranked in the order of decrease of Δ-variance and the top 2000 probes were selected for further training and testing of the RF machine learning algorithm in order to build a model discriminating between BRCAness-positive and BRCAness-negative classes of OS samples. The results of this model testing and training evaluated with ROC curves and confusion matrices are shown (Fig 1A and Fig B in S1 File and Table 1).

Download:

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

Fig 1.

ROC curves for each cross validation fold (dashed-black line) and the mean ROC curve (solid black line) for random forest models fitted (A) to 43 samples and (B) to 41 samples where two samples, HD_00J and I047_005 were removed (the dashed diagonal line represent AUC = 0.5, i.e. classification at random). (C) ROC curve for independent validation test (AUC = 0.87) (n = 20).

https://doi.org/10.1371/journal.pcbi.1009562.g001

Download:

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

Table 1. Confusion matrices of each cross validation and average fold for random forest models.

https://doi.org/10.1371/journal.pcbi.1009562.t001

The inherited noisiness of the OS samples provided negative Δ-variances for the majority of probes. In the 43 sample set we have only 1,055 probes with positive Δ-variance (~0.25%) and 426,520 probes with negative Δ-variances. With this settings taking the top 2000 probes led to misclassification of two samples at unsupervised hierarchical clustering (Fig 2A).

Download:

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

Fig 2.

Initial unsupervised hierarchical clustering of 2000 probes (A) in all available 43 samples and (C) in 41 sample set after removal of two samples (HD_00J and I047_005). Note that these two samples were clustered incorrectly. Lower part in heatmap (A) shows probes with very low intensities across all samples yielding delta variances with negative values. Brown band: BRCAness-positive, blue band: BRCAness negative samples. (B) Barplot shows the number of positive Δ-variances re-computed after excluding each sample consequently from 43 sample set. Red—BRCAness negative, green—positive samples, Y-axe—number of positive Δ-variances. Subsequent removal of HD_00J and I047_005 gives the largest increase in the number of positive Δ-variances.

https://doi.org/10.1371/journal.pcbi.1009562.g002

In order to maximize the number of informative predictors with positive Δ-variances, without significantly reducing the number of samples, we computed the number of probes with positive Δ-variances by removing one sample at a time from the 43 sample set and re-calculating the number of positive variances (Fig 2B and Table A in S1 File). Subsequent removal of two samples, I047_005 and HD_00J, in BRCAness positive and in BRCAness negative group, correspondingly showed the largest number of probes with positive Δ-variances. These two samples were the ones misclassified at unsupervised hierarchical clustering (Fig 2A).

Therefore, we assumed that the negative Δ-variances decreased the accuracy of the prediction and only probes with positive Δ-variances were informative. After removal of samples I047_005 and HD_00J, unsupervised hierarchical clustering correctly grouped 41 samples and the number of probes with positive Δ-variances increased to 3685 (Fig 2C). The accuracy estimated by AUC and confusion matrices improved for the 41 OS sample set shown in Fig 1B and Fig B in S1 File and Table 1 as well as grouping in hierarchical clustering (Fig 2C).

The OOB error evaluates model performance on data not participating in the training of RF.However, due to selection bias the OOB error computed on a model fitted on the features selected just once might lead to overfitting and yield over-optimistic error estimation [51]. To avoid this bias, a Δ-variance feature selection algorithm was used internally in a five-fold cross-validation loop [39]. Within each fold, feature selection was conducted de novo. The model fitted using the selected features anew each time and was tested against the fold set aside.

During a five-fold cross-validation loop the feature selection, i.e. selection of 2000 top probes was performed independently according to the filter method as opposed to the wrapper method. The filter method picks up inherent properties of the features based on univariate statistics [52]. Since 1/5 of the samples were excluded from training and left for testing purposes, i.e. not participating in computing of Δ-variances, the set of 2000 probes was fluctuating after each iteration. The 41 sample set showed improved accuracy evaluated with AUC and confusion matrix (Fig 1B and Fig B in S1 File and Table 1). Unsupervised hierarchical clustering of uncorrelated 246 probes/predictors, extracted after each cross-validation loop in the course of the feature selection and, subsequently, merged into a non-redundant set, channeled OS samples to corresponding BRCAness-positive and BRCAness-negative groups (Fig 3A). Also, correct grouping was observed when 4052 probes correlated to and including these 246 probes (rho > 0.8), were merged (Fig 3B) and subjected to unsupervised hierarchical clustering. A pipeline for data processing, RF classifier training and validation is shown in Fig C in S1 File.

Download:

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

Fig 3. Unsupervised hierarchical clustering of samples and probes in 41 sample set.

(A) 246 probes identified as important non-correlated predictors combined after 5-fold cross validation and (B) 4052 probes correlated to and including 246 probes discriminate between BRCAness-positive and BRCAness-negative samples. BRCAness-positive and BRCAness-negative samples indicated by brown and blue horizontal bars at the top, correspondingly. Please note, for the final RF model, used for validation with independent test set (n = 20), all 41 samples were used as one training set. This set yielded 54 uncorrelated probes out of 2000 probes.

https://doi.org/10.1371/journal.pcbi.1009562.g003

The ROC curve is a plot visualizing the true positive (TP) rate against the false positive (FP) rate at various threshold settings. The average of the area under the ROC curve (AUC) across all 5 CV folds with 43 samples was 0.88 (Fig 1A). However, after removing two samples, I047_005 and HD_00J, the average AUC increased to 0.95 (Fig 1B). Also, the combined accuracy after 5-fold CV increased from 77% to 85% in 41 sample set compared to 43 sample set (Table 1). Keeping in mind that maximum AUC is 1.0 when prediction is 100% accurate and, AUC is 0.5 with random prediction, the average AUC of the DNA-methylation based RF classifier 0.95 demonstrates accurate prediction of BRCAness in OS samples. Additionally, the RF classifier trained on 41 OS samples was validated with an independent OS test set (Fig 1C).

The arrangement of samples in order of increase of PGC, the measure of BRCAness obtained from genomic NGS data, shows a pattern of probe intensity decrease with increase of PGC (Fig D in S1 File). This suggests that gene expression is repressed through higher DNA-methylation in BRCAness negative samples and, conversely, in BRCAness positive samples expression of certain genes are upregulated. However, we cannot exclude the possibility that the reason for lower DNA methylation signal in BRCAness positive samples was simply because of massive genome alterations affecting DNA sites of methylation.

In line with a standard operating procedure for RF, we assigned a fraction of RF votes as a measure of DNA repair defect, where a fraction of RF votes > 0.5 classifies a sample as BRCAness positive.

RNA-Seq sample selection and differential expression analysis

In order to assess heterogeneity of RNA-Seq samples, Principle Component Analysis (PCA) was applied (Fig E in S1 File). The overlap of ellipses denoting the confidence interval showed poor separation of the two groups. As we expected, the two samples, HD_00J and I047_007 (already removed in RF), were far away from their corresponding clusters, and for this reason these RNA-Seq samples were excluded from subsequent consideration. Out of 41 OS samples used in methylation classifier, 31 samples had RNA-Seq data. To avoid samples in the "grey zone", i.e. around the border line separating BRCAness-negative group from BRCAness-positive group, we selected 2/3 of the top samples with the highest score from the BRCAness-positive group and 2/3 of the bottom samples with the lowest score from the BRCAness-negative group ranked according to Nexus PGC and used their corresponding RNA-Seq data if they were available. This approach excludes overlapping and border line samples. This approach yielded 13 BRCAness-positive and 7 BRCAness-negative RNA-Seq samples. To assess the heterogeneity of RNA-Seq samples, unsupervised hierarchical clustering was applied (Fig 4A). It shows each individual RNA-seq sample (leaf) joined by the branches. The height represents the distance value from the input Euclidean distance matrix. One can consider the first split in the dendrogram to be the two largest distinct groups, with more groups forming at each subsequent split down the dendrogram. In Fig 4, the demarcated boxes in red and cyan represent two groups obtained by cutting the dendrogram at even height. Moving top-down, the group in the red box and I007_018 are separated from the rest of the groups at the top level. Thus, it is expected to find most difference between sets of samples separated at the top level. Considering the right side of the first split, the outlier I007_018 is separated from the group of samples in the red box and hence it is not included in the final set of RNA-seq samples. A similar rationale was applied to I062_007, even though I062_007 belongs to the same class as the samples in the cyan box. BRCAness-positive RNA-Seq MS-00A was excluded from the stretch of BRCAness-negative samples. Unsupervised clustering and PCA of the remaining 17 RNA-Seq samples, 12 BRCAness-positive and 5 BRCAness-negative, showed that two clusters were well-separated (Fig 4B and Fig E in S1 File). Subsequently, these 17 RNA-Seq samples were used for differential expression (DE) and gene set enrichment analyses (GSEA).

Download:

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

Fig 4.

Dendrogram showing the results of our hierarchical cluster analysis of (A) 20 RNA-Seq samples selected from two thirds top and bottom of the OS sample list arranged according to PGC score; and (B) 17 remained RNA-Seq samples used for DE analysis and GSEA. BRCAness-positive samples are denoted by red color and BRCAness-negative by cyan color. Leaves representing excluded samples are denoted with ’x’. The demarcated boxes in red and cyan represent two groups obtained by cutting the dendrogram at even height and used further in DGE and GSEA. The height of a group is defined by Ward’s method.

https://doi.org/10.1371/journal.pcbi.1009562.g004

There were 1,528 significantly DE genes consisting of 449 genes upregulated in BRCAness-positive group and 1,079 genes downregulated. Among the genes upregulated in BRCAness-positive versus BRCAness-negative OS tumors were genes important in DNA replication, DNA repair, cell cycle regulation such as BLM, BRCA1, EME1, ERCC1, EYA1, FANCB, H2AX (aka H2AFX), HIST1H2BC, HIST1H2BD, HIST1H2BN, LIG3, POLD2, POLH, RPA3, SMUG1, TIPIN, TOP3A, TIMELESS. Of note, the upregulated histone H2AX was extensively studied as a DNA damage response factor. The change between γH2AX dephosphorylation and H2AX phosphorylation marks the switch between cell repair and apoptotic responses, correspondingly [53]. γH2AX foci formation correlates with a double strand break induced by irradiation exposure [54].

The majority of the gene expression signatures downregulated in BRCAness-positive versus BRCAness-negative OS tumors were those associated with the immune response indicating the level of activity of the immune cells in the OS microenvironment. The genes related to adaptive immunity such as T cell receptor genes (CD3E, CD3D, CD3G) and B cell genes (CD79A, CD83) as well as genes of the major histocompatibility complex (MHC) such as HLA-DOB, HLA-DQA1, HLA-DQB1 were significantly downregulated in BRCAness-positive OS tumors. The NF-kB, composed of a family of genes (NFKB1, NFKB2, RELB, REL), an important transcription factor related to T-cell activation and cytokines expression, was significantly downregulated or showed strong tendency towards downregulation in BRCAness-positive versus BRCAness-negative OS tumors. Expression of genes specific to cytotoxic activities mediated by T cells, such as CXCR3, PRF1, KLRF1, GZMH, GZMK was also significantly reduced in BRCAness-positive OS tumors as well as expression of the genes associated with immune response and proinflammatory cytokines (IL1A, IL18, CCL11, CCL22, CCL13) (Table A in S2 File).

Gene set enrichment analysis

The Broad Institute GSEA tool reported gene sets with a positive enrichment score, i.e. gene sets that display enrichment at the top of the pre-ranked list, were hallmark dna repair, kegg dna replication, kegg mismatch repair, kegg homologous recombination (Table 2 and Fig F in S1 File). Genes contributing most to the enrichment results in gene sets related to DNA repair and damage tolerance are listed in Table B in S1 File. The second portion of the GSEA report showed results for gene sets with a negative enrichment score, i.e. gene sets that display enrichment at the bottom of the pre-ranked list. Those signatures included hallmark allograft rejection, kegg asthma, kegg graft versus host disease, kegg allograft rejection (Table 2 and Fig G in S1 File). Genes contributing most to the enrichment results in gene sets related to immune system disease category are listed in Table C in S1 File. Asthma, allograft rejection, and graft versus host disease signatures reflect downregulation of MHC classes I and II which in turn reduces responsiveness of anti-tumor immune cells and is typical tumor evasion of the immune response [55,56] (GSEA results in the supplement).

Download:

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

Table 2. Enrichment results for gene sets.

https://doi.org/10.1371/journal.pcbi.1009562.t002

As additional GSEA we applied the web-tool Enrichr (https://amp.pharm.mssm.edu/Enrichr/) where statistically significant upregulated and downregulated genes were uploaded as separate gene sets. Enrichr with upregulated genes and KEGG database reported Fanconi anemia pathway and HR, both pathways share most of genes such as BRCA1, POLD2.Enrichr with the Reactome database reported multiple signatures related to HR and DNA-repair such as homology directed repair, DNA double-strand break repair. For the set of downregulated genes Enrichr reported MHC II genes in graft-versus-host-disease, allograft rejection and in other signatures (Enrichr results in the supplement).

Validation of classifier with independent test set

Already after training BRCAness classifier, we obtained OS samples (n = 20) through INFORM ongoing collection of samples with relapsed or refractory osteosarcoma. This test set was used for validation of BRCAness classifier. Samples were separated to BRCAness-positive (n = 14) and BRCAness-negative (n = 6) groups based on FREEC-HRDtools PGC score. New methylation data arrays were normalized together with training set. Subsequently, trained model was applied to test set for classification. Only single sample in each group was misclassified (Table 3). The validation with independent OS set resulted in AUC of 0.87, accuracy of 0.9, sensitivity and specificity of 0.93 and 0.83, correspondingly (Fig 1C). Differential expression analysis was performed with the RNA-Seq test set where all samples classified as BRCAness-positive (n = 14) and BRCAness-negative (n = 6) were grouped conditional on the presence or the absence of BRCAness characteristics, tested for DE, ranked, and subsequently analyzed for pathway enrichment. For the set of upregulated genes Reactome database reported multiple signatures related to HR and DNA-repair such as homologous DNA-pairing and strand exchange, homology directed repair through homologous recombination, DNA-double strand break repair and homology directed repair (Table D in S1 File). KEGG database reported cell cycle, dna replication signatures as well as signatures related to DNA-repair, specifically, mismatch repair (fdr < 0.05), homologous recombination repair (fdr < 0.05) for the set of upregulated genes. For the set of downregulated genes KEGG reported MHC II genes in graft-versus-host-disease, allograft rejection and in other signatures including apoptosis (Table E in S1 File). Matching signatures were observed in training set GSEA.

Download:

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

Table 3. Test set confusion matrix.

https://doi.org/10.1371/journal.pcbi.1009562.t003