1. Systematic review and selection of studies

A systematic search of studies published in the period 2011-2019 was conducted in October 2019 following This review was conducted in October 2019, according to the Preferred Reporting Items for Systematic Reviews and Meta-Analyses (PRISMA) statement guidelines [52]. We searched breast cancer methylation data studies in GEO database [53], using the keywords: “breast cancer”, “methylation” and “Homo sapiens" in studies published in English. We applied the following exclusion criteria: (i) studies conducted in organisms other than humans; (ii) sample size less than 10 in each experimental group; (iii) different methylation profiling platform to Infinium HumanMethylation450 BeadChip from Illumina and (iv) experimental design different from case-control.

2. Bioinformatics analysis strategy

The following procedure was applied to each of the individually selected studies: i) Data acquisition. ii) Exploratory analysis and quality control of the samples. iii) Then, the studies were separated into different "case vs. control" comparisons, i.e. if a study consists of healthy samples and samples from two breast cancer subtypes we will obtain two comparisons —healthy samples vs. breast cancer subtype I and healthy samples vs. breast cancer subtype II—. iv) Analysis of differential methylation between case and control groups by comparison. v) Analysis of functional enrichment for each of the comparisons. vi) Integration of the methylation profiles and functional results in the final meta-analysis

3. Data exploration, quality control and normalization

Exploration of raw data was performed through principal component analysis (PCA) and clustering analysis. Data quality control was conducted with minfi [24] R package. In this step we check levels of the signal from the methylated and non-methylated channels, then a quantile normalization was performed on the data to compare all together

4. Individual epigenomic analysis

The differential methylation analysis for each of the different comparisons was carried out with the minfi R package. For each comparison, permutation tests were performed to obtain the methylation scores of all analyzed regions on the Illumina BeadChip 450k. A total of 7 comparisons derived from the individual studies were performed: Control vs TNBC, Control vs HER-2, Control vs Lum-A, Control vs Lum-B, Control vs Brcmutant, Control vs Invasive Ductal Carcinoma —IDC— and Control vs Invasive Lobular Carcinoma —ILC—.

Subsequently, methylation scores were annotated to gene level through the bumphunter R package [54]. In cases where more than one differentially methylated region (DMR) existed for a gene, the one with the highest absolute value of differential methylation was used. Then, a functional enrichment analysis was performed from the differential methylation results using Gene Set Analysis (GSA) [55]. To do this, the genes were ordered according to their p-values and the sign of the contrast statistic. The GSA was then performed following the logistic regression model implemented in the mdgsa R package [56], as well as its corresponding functional annotation. The p-values were corrected for each function for false discovery rate (FDR) [57]. The databases used for functional enrichment were the Gene Ontology [58] (GO) and the Kyoto PATHWAY Encyclopedia of Genes and Genomes [59] (KEGG).

Significant functions have been represented in the form of Upset plots [60]. With this type of graph we can see the number of functional elements specific and shared by each of the different breast cancer subtypes studied

5. Methylated genes meta-analysis

From the mapping of the differentially methylated regions, the genes included in them were identified. In order to obtain an integrated level of methylation in all the studies, those genes with a common differential pattern between the case and control groups were selected. For each gene, the p-values of all comparisons were then combined using the Fisher combination method (or the inverse normal/weighted method). This meta-analysis strategy provided the genes with a significant common methylation profile in all the selected studies

6. Functional DNA methylation signatures meta-analysis

Finally, once the functional enrichment study was performed for each of the individual comparisons the results were integrated into a functional meta-analysis. The methodology of the functional meta-analysis followed was the one proposed by [17]. For this task, the metafor [61] R package was used to assess the combined effect of the studies together with a random effects model. The results of the meta-analysis are more robust than the results of individual comparisons due largely to the larger sample size. The variability of individual studies was taken into account in the meta-analysis so that if a study had little variability, it would carry more weight in the meta-analysis and be more influential in the calculation of the log odds ratio (LOR). In turn, both a heterogeneity analysis to check the suitability of the different studies selected and a sensitivity analysis and assessment of bias to detect whether any of the comparisons had an excessive influence on the final meta-analysis were performed.

Each function analysed in the meta-analysis is accompanied by the combined estimate of the effect of the studies (LOR), the 95% confidence interval and the p-value adjusted by the Benjamini and Hochberg method. Thus, those functions with an adjusted p-value equal to or less than 0.05 were considered significant. For each significant function, forest and funnel plots were used to measure the contribution of each study to the meta-analysis and to assess its variability

Fig. 1. Data-analysis workflow.

© 2020 by the authors. Submitted for possible open access publication under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).