MOMIC: A multi-omics pipeline for data analysis, integration and interpretation

  Citation: Madrid-Márquez, L.; Rubio-Escudero, C.; Pontes, B.; González-Pérez, A.; Riquelme, J.C.; Sáez, M.E. MOMIC: A Multi-Omics Pipeline for Data Analysis, Integration and Interpretation. Appl. Sci. 2022,12, 3987. https://doi.org/ 10.3390/app12083987 Academic Editors: Konstantinos E. Psannis and Christos L. Stergiou Received: 20 January 2022 Accepted: 2 April 2022 Published: 14 April 2022 Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations. Copyright: © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/). applied sciences Article MOMIC: A Multi-Omics Pipeline for Data Analysis, Integration and Interpretation Laura Madrid-Márquez 1,2 , Cristina Rubio-Escudero 2,* , Beatriz Pontes 2, Antonio González-Pérez 1, José C. Riquelme 2and Maria E. Sáez 1 1Andalusian Bioiformatics Research Centre (CAEBi), 41013 Sevilla, Spain; [email protected] (L.M.-M.); [email protected] (A.G.-P.); [email protected] (M.E.S.) 2Department of Computer Languages and Systems, University of Sevilla, 41013 Sevilla, Spain; [email protected] (B.P.); [email protected] (J.C.R.) *Correspondence: cr[email protected] Abstract: Background and Objectives: The burst of high-throughput omics technologies has given rise to a new era in systems biology, offering an unprecedented scenario for deriving meaningful biological knowledge through the integration of different layers of information. Methods: We have developed a new software tool, MOMIC, that guides the user through the application of different analysis on a wide range of omic data, from the independent single-omics analysis to the combination of heterogeneous data at different molecular levels. Results: The proposed pipeline is developed as a collection of Jupyter notebooks, easily editable, reproducible and well documented. It can be modified to accommodate new analysis workflows and data types. It is accessible via momic.us.es, and as a docker project available at github that can be locally installed. Conclusions: MOMIC offers a complete analysis environment for analysing and integrating multi-omics data in a single, easy-to-use platform. Keywords: multi-omics; integration; analysis; pipeline; reproducibility 1. Introduction Recent technological advances have allowed the generation of vast amounts of data from different omics layers, offering information from diverse aspects of cell biology ranging from genetic variation to the quantification of the proteins encoded by these genes [1]. This fact has given rise to a new era in systems biology, generating the need of integrating different data types. However, the methodological spectrum for the analysis of biological data, and the theoretical framework required to interpret the obtained information is lagging behind [ 2 ], and consequently, few standardized workflows for integrative data analyses have been proposed. Integrative methods can collectively mine multiple types of biological data and provide richer insights. A variety of software applications and analytical frameworks [ 3 , 4 ] are already available for handling diverse type of data. However, the lack of a central repository which catalogs, links, rates and summarizes these tools is one of the factors limiting their use by the research community. Furthermore, some of these resources present notable inconveniences: difficult installation, moderate learning curve, limited customisation or paid license. Increasing researcher awareness and accessibility to tools will improve their use and uptake [5]. This article presents MOMIC, a comprehensive pipeline capable of analyzing and summarising single-omics data by means of meta-analysis, as well as performing integrative analysis of different molecular levels. Enrichment analysis and visualization methods have also been implemented to aid data interpretation. The tool comprises a set of Jupyter notebooks [ 6 ] designed for each molecular level, based on publications and best practices available in the literature. An important aspect of MOMIC is its reproducibility; all steps Appl. Sci. 2022,12, 3987. https://doi.org/10.3390/app12083987 https://www.mdpi.com/journal/applsci Appl. Sci. 2022,12, 3987 2 of 17 are kept in the same notebook, including commands that are traditionally executed in a terminal, so the user can always see past outputs and exact parameters. It is presented as a web tool and as a Docker project to be installed locally. The article is structured as follows: Section 2presents the materials and methods used in this work, detailing all the protocols applied. Section 3explores existing alternatives and compares them to MOMIC. In Section 4the authors show diverse analysis outcomes on synthetic and public datasets and how MOMIC has been successfully used in the data analysis required for the ADAPTED IMI project [ 7 ]. Finally, Section 6presents the final conclusions and remarks related to this work. 2. Materials and Methods MOMIC currently compiles protocols for whole genome SNP data (GWAS), mRNA expression (both from arrays and RNAseq experiments) and protein data. Along with enrichment analysis and methods for combining distinct data. The data analysis workflows starts with the pre-processing and quality control of the individual datasets. Each of the datasets are then independently analyzed following different protocols designed for each type of data. To consolidate these independent results for the same omic, a meta-analysis can be performed, generating single lists of SNPs or genes. After this, one can take a step forward and combine diverse omic results, along with functional genomics analysis, and try to elucidate the potential causative changes that lead to disease. With this integrative approach the user can get a single ranked list of candidate genes summarising the heterogeneous data at different molecular levels. Finally, visualization and pathway analysis tools are available to further explore analysis results. This whole procedure has been summarized in Figure 1. Figure 1. MOMIC pipeline steps. Appl. Sci. 2022,12, 3987 3 of 17 2.1. Core Software MOMIC is presented as a collection of Jupyter notebooks using JupyterLab and Jupyter- Hub on top, written mainly in R language and containerized in docker [ 8 ]. JupyterLab is an open-source web-based interactive development environment that allows users to create and share documents, called notebooks, that contain live code, equations, visualizations and narrative text. JupyterHub brings the power of notebooks to multiple users. MOMIC is distributed as a docker-compose project that contains the instructions needed to automatically create a fully working pipeline server with JupyterLab. These instructions assemble commands that enable convenient extensions in JupyterLab, like git for notebook version controlling and table of content, set up docker volumes for data persistence, arrange the pipeline source code and all the necessary libraries and third-party software. An alternative to a local installation is to use the pipeline hosted at momic.us.es. Note this alternative is intended for light analysis or quick testing. Login details are specified in the user manual and upon request. The collection of notebooks are provided as read-only, serving as templates. The user can replicate the same analysis duplicating the template and modifying it according to his needs. Notebooks can be downloaded as ipython notebooks, HTML, PDF or scripts. The later exports the notebook in scripting language that can be reused at the user convenience. 2.2. Implemented Protocols We describe below the protocols developed. These are executed on publicly available datasets and synthetic data related to Alzheimer Disease and serves as an illustrative example of how to use the pipeline to carry out customised analysis. Figures S1–S4 from the Supplementary Material precisely describe the collection of notebooks designed for each protocol, defining the individual tasks, the flow of information and expected inputs and outputs. The user manual and the jupyter notebooks contain more detailed explanations and documented source code. Genome Wide Expression Studies (GWES) from microarrays. RNA microarrays can simultaneously measure the expression level of thousands of genes within a particular mRNA sample [ 9 ]. This can be used to compare the level of gene transcription in clinical or biological conditions in order to find differences in expression levels between predefined groups of samples. This pipeline starts from raw expression data and ends with a set of differentially expressed genes, as shown in Figure 2. In order to capture these significant genes, the standard protocol to follow comprises data pre-processing, differential expression analysis (DE) and annotation steps. Pre-processing involves transforming raw data into a matrix of normalized intensities, using the RMA [ 10 ], limma [ 11 ] or preprocessCore [ 12 ] R packages depending on the data platform. Batch effect removal is also available through the combat function of the sva package [ 13 ]. DE is performed with limma, an R package based on linear models for microarray data, including optional correction for covariates, followed by annotation to translate probes to genes, using Entrez and Human Genome Nomenclature Committee (HGNC) identifiers. Additionally, heatmaps, PCA and volcano plots can be easily generated to assist the interpretation of results. Genome Wide Expression Studies (GWES) from RNASeq. RNA-Seq is a particular technology-based sequencing technique which uses next-generation sequencing (NGS) to reveal the presence and quantity of RNA in a biological sample at a given moment. This pipeline starts from raw sequence reads, and ends with a set of differentially expressed genes (see Figure 3). The steps to complete this analysis are: quality check of the raw reads using fastqc [ 14 ], alignment of reads against the reference genome using STAR [ 15 ], quality control of aligned reads (by inspection of QC plots) and reads quantification with STAR. Differential expression can be then performed using DESeq2 [ 16 ], an R package based on a model using the negative binomial distribution, followed by annotation with biomaRt library [17]. As for microarrays, plots like MA, heatmap and PCA can be easily obtained. Appl. Sci. 2022,12, 3987 4 of 17 Figure 2. Microarray analysis workflow. Figure 3. RNA-Seq analysis workflow. Appl. Sci. 2022,12, 3987 5 of 17 Genome Wide Association Studies (GWAS): GWAS is an observational study of a genome-wide set of genetic variants in different individuals to examine if any variant is associated with a trait. The GWAS protocol is available for genome builds GRCh37/hg19 and GRCh38/hg38, including the UCSC liftOver tools for updating genetic positions. This pipeline starts from genomic data in PLINK format, and ends with Single Nucleotide Polymorphisms (SNPs) annotated for surrounding genes and rs identifiers (dbSNP v150 for the GRCh37 pipeline and dbSNP v151 for the GRCh38 pipeline), which are then aggregated to the level of whole genes, testing the joint association of all markers in the gene with the phenotype, as depicted in Figure 4. Procedure followed here is a composition from Anderson et al. protocol [ 18 ] and AT Marees et al. [ 19 ]. Steps are: initial quality control excluding bad quality individuals and SNPs, genotype imputation using the Michigan [ 20 ] or TopMed [ 21 ] imputation servers, case control association study using PLINK [ 22 ] and gene-wise statistics computed using MAGMA [23]. SNP level data can be visualised with Manhattan and QQ plots [24]. Figure 4. GWAS analysis workflow. Proteomics: Differential expression between samples can be also explored in MOMIC at the protein level. This pipeline starts from an intensity matrix, generated by specific software such as MaxQuant/Perseus, and ends with a set of differentially expressed proteins (see Figure 5). The first step is data processing, which includes removing decoy matches and matches to contaminant, extracting the LFQ intensities columns, filtering on missing values, log transformation, normalisation, unique peptide counting and visualisation. Appl. Sci. 2022,12, 3987 6 of 17 Differential expression is performed similarly to microarray DE analysis using DEqMS R package [ 25 ], which is built on top of limma and accounts for variance dependence on the number of quantified peptides. Heatmap and volcano plots are also generated to assist the interpretation of results. Figure 5. Proteomics analysis workflow. Meta-analysis: Meta-analysis is the statistical procedure for synthesising data across studies addressing the same question at the same molecular level. Two different one-step protocols are provided here, one for combining multiple GWAS results using the software METAL [ 26 ] and another for transcriptomics and proteomics results utilising MetaDE R library [ 27 ]. The former starts from SNP level GWAS results and ends with a single list of SNPs and associated summary statistics. The later starts from the results of the DE analysis of various studies, and ends with an unified list of differentially expressed genes or proteins across studies with their correspondent average effect sizes and p-values. Integrative analysis: Integrative analysis aims at consolidating heterogeneous data at different omics levels to understand their interrelation and combined influence on the disease processes. This one-step protocol starts from the results of the meta-analysis from different studies, and ends with a unified ranked list of differentially expressed genes. It is performed using the Robust Rank Aggregation method [ 28 ]. It detects genes that are ranked consistently better than expected under the null hypothesis of uncorrelated inputs and assigns a significance score for each gene. Enrichment: Enrichment analysis can identify gene onthology (GO) terms and pathways which are statistically over or under-represented within the list of interest, by systematically mapping genes and proteins to their associated biological annotations. The enrichment is performed using WebGestalt [ 29 ] in R and the visualization with GOplot [ 30 ] and pheatmap from R CRAN. 2.3. Availability and Requirements MOMIC is free to use and it is distributed under the MIT License. It is hosted at the University of Seville data center at momic.us.es. Sources, documentation and guides for users and administrators are available on GitHub: MOMIC server: https://github.com/laumadmar/MOMIC_server.git (accessed on 10 January 2022) Collection of notebooks: https://github.com/laumadmar/MOMIC_notebooks.git (accessed on 10 January 2022) Appl. Sci. 2022,12, 3987 7 of 17 Documentation: https://laumadmar.github.io/MOMIC_server (accessed on 10 January 2022) If hosting MOMIC locally, the minimum recommendation is 15 GB of RAM which may need to be increased with large datasets, specially for some processes during an RNASeq Analysis or GWAS. The most computationally intensive tasks are: the preparation of the 1 KG [ 31 ] data prior the population stratification on a GWAS analysis, which takes up to 8 GB of virtual memory and takes 26 min to complete, and a GWAS association analysis taking up to 45 MB and 45 min. Disk space can vary depending on the size of the user data, 260 GB is the actual size of the data volume for momic.us.es. Suplementary Table S1 summarises the data types used as inputs and their sizes. MOMIC has been tested on Ubuntu 18.04, CentOS 8, Rocky Linux 8.5, Windows 10 and macOS High Sierra. 3. Motivation and Existing Alternatives MOMIC emerged out of the necessity of exploring data at different molecular levels as part of the ADAPTED project, reviewed in the Results section. The main requirements regarding the bioinformatics analysis were to perform them in a comprehensive way, easy to share and fully reproducible while keeping good evidence of the work done. After examining existing tools none of them quite suited the above requirements and therefore MOMIC was developed. Table 1shows some of these tools and the main differences with MOMIC, being Galaxy and GenePattern the ones that serve similar functionality. Despite Galaxy being a powerful and useful tool, it has some limitations compared to MOMIC in terms of accessibility to the source code and ease of data manipulation and customisation of modules. GenePattern is a great tool for users with no programming experience but not that well suited for those who need full control over the source code and straight forward customisation. It offers a good variety of analysis and modules but lacks GWAS and integration of different omics. Appl. Sci. 2022,12, 3987 8 of 17 Table 1. Comparison of existing tools for omic analyses. Name Platform Friendliness Functionality Comparison with MOMIC Availability Offers a platform for Similar to MOMIC in functionality https://www.genepattern.org/ GenePattern Web Easy reproducible bioinformatics but the code is closed source, (accessed 1 April 2022) making it impossible to full customise the pipelines, however changing the input parameters does allow for minor changes. It focuses mainly on transcriptomics and lacks GWAS analysis and integration of omics. GenePattern notebook also extends JupyterHub and it is available via web and as a local server Enables researchers without informatics Serves the same aim but unlike MOMIC, https://usegalaxy.org/ Galaxy Web Medium expertise to perform computational it does not offer real time code (accessed 1 April 2022) analyses through the web and visualisation, easy data manipulation or customisation. Developing new tools in Galaxy is not straight forward and requires XML skills Offers a wide range of novel This R package can be installed http://mixomics.org/ mixOmics R Difficult multivariate methods for the exploration in MOMIC to extend the proposed (accessed 1 April 2022) and integration of biological datasets transcriptomics pipelines. It does not with a particular focus on offer protocols for GWAS analysis on variable selection Appl. Sci. 2022,12, 3987 9 of 17 Table 1. Cont. Name Platform Friendliness Functionality Comparison with MOMIC Availability Offers integrative visualization Requires specific input format http://www.paintomics.org/ Paintomics 3 Web Easy of multiple omic datasets and data wrangling. This just (accessed 1 April 2022) onto KEGG pathways provides visualisation, as oppose to MOMIC which also provides a platform for data analysis. Results obtained from MOMIC could be used to feed this tool if desired. It does not offer protocols for GWAS analysis Offers interactive NGS analysis Good alternative for NGS analysis https://www.basepairtech.com/ Basepair Web Easy pipelines for users with analysis but there is a fee (accessed 1 April 2022) no programming experience to pay per sample. The source code is not exposed so the customisation is limited Appl. Sci. 2022,12, 3987 16 of 17 6. Kluyver, T.; Ragan-Kelley, B.; Pérez, F.; Granger, B.; Bussonnier, M.; Frederic, J.; Kelley, K.; Hamrick, J.; Grout, J.; Corlay, S.; et al. Jupyter Notebooks—A publishing format for reproducible computational workflows. In Positioning and Power in Academic Publishing: Players, Agents and Agendas; Loizides, F., Schmidt, B., Eds.; IOS Press: Amsterdam, The Netherlands, 2016; pp. 87–90. 7. Madrid, L.; Moreno-Grau, S.; Ahmad, S.; González-Pérez, A.; de Rojas, I.; Xia, R.; Adami, P.V.M.; García-González, P.; Kleineidam, L. ; Yang, Q.; et al. Multiomics integrative analysis identifies APOE allele-specific blood biomarkers associated to Alzheimer’s disease etiopathogenesis. Aging 2021,13, 9277. [CrossRef] 8. Merkel, D. Docker: Lightweight linux containers for consistent development and deployment. Linux J. 2014,2014, 2. 9. Baldi, P.; Hatfield, G.W. DNA Microarrays and Gene Expression: From Experiments to Data Analysis and Modeling; Cambridge University Press: Cambridge, MA, USA, 2011. 10. Gautier, L.; Cope, L.; Bolstad, B.M.; Irizarry, R.A. affy—Analysis of Affymetrix GeneChip data at the probe level. Bioinformatics 2004,20, 307–315. [CrossRef] 11. Smyth, G.K. Limma: Linear models for microarray data. In Bioinformatics and Computational Biology Solutions Using R and Bioconductor; Gentleman, R., Carey, V.J., Huber, W., Irizarry, R.A., Dudoit, S., Eds.; Springer: New York, NY, USA, 2005; pp. 397– 420. 12. Bolstad, B. preprocessCore: A Collection of Pre-Processing Functions; R Package Version 1.50.0; Bioconductor: Santo Domingo, Dominican Republic, 2020 . 13. Leek, J.T.; Johnson, W.E.; Parker, H.S.; Jaffe, A.E.; Storey, J.D. The Sva Package for Removing Batch Effects and Other Unwanted Variation in High-Throughput Experiments. Bioinformatics 2012,28, 882–883. [CrossRef] 14. Andrews, S.; Krueger, F.; Segonds-Pichon, A.; Biggins, L.; Krueger, C.; Wingett, S. FastQC; Babraham Institute: Babraham, UK, 2010. 15. Dobin, A.; Davis, C.A.; Schlesinger, F.; Drenkow, J.; Zaleski, C.; Jha, S.; Batut, P.; Chaisson, M.; Gingeras, T.R. STAR: Ultrafast universal RNA-seq aligner. Bioinformatics 2013,29, 15–21. [CrossRef] 16. Love, M.I.; Huber, W.; Anders, S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014,15, 550. [CrossRef] [PubMed] 17. Durinck, S.; Spellman, P.T.; Birney, E.; Huber, W. Mapping identifiers for the integration of genomic datasets with the R/Bioconductor package biomaRt. Nat. Protoc. 2009,4, 1184–1191. [CrossRef] [PubMed] 18. Anderson, C.A.; Pettersson, F.H.; Clarke, G.M.; Cardon, L.R.; Morris, A.P.; Zondervan, K.T. Data quality control in genetic case-control association studies. Nat. Protoc. 2010,5, 1564–1573. [CrossRef] [PubMed] 19. Marees, A.T.; de Kluiver, H.; Stringer, S.; Vorspan, F.; Curis, E.; Marie-Claire, C.; Derks, E.M. A tutorial on conducting genomewide association studies: Quality control and statistical analysis. Int. J. Methods Psychiatr. Res. 2018 ,27, e1608. [CrossRef] [PubMed] 20. Das, S.; Forer, L.; Schönherr, S.; Sidore, C.; Locke, A.E.; Kwong, A.; Vrieze, S.I.; Chew, E.Y.; Levy, S.; McGue, M.; et al. Nextgeneration genotype imputation service and methods. Nat. Genet. 2016,48, 1284–1287. [CrossRef] [PubMed] 21. Taliun, D.; Harris, D.N.; Kessler, M.D.; Carlson, J.; Szpiech, Z.A.; Torres, R.; Taliun, S.A.G.; Corvelo, A.; Gogarten, S.M.; Kang, H.M.; et al. Sequencing of 53,831 diverse genomes from the NHLBI TOPMed Program. Nature 2021 ,590, 290–299. [CrossRef] 22. Purcell, S.; Neale, B.; Todd-Brown, K.; Thomas, L.; Ferreira, M.A.; Bender, D.; Maller, J.; Sklar, P.; De Bakker, P.I.; Daly, M.J.; et al. PLINK: A tool set for whole-genome association and population-based linkage analyses. Am. J. Hum. Genet. 2007 ,81, 559–575. [CrossRef] 23. de Leeuw, C.A.; Mooij, J.M.; Heskes, T.; Posthuma, D. MAGMA: Generalized gene-set analysis of GWAS data. PLoS Comput. Biol. 2015,11, e1004219. [CrossRef] 24. Turner, S.D. qqman: An R package for visualizing GWAS results using QQ and manhattan plots. Biorxiv 2014 , 005165. [CrossRef] 25. Zhu, Y. Bioconductor-DEqMS: A Tool to Perform Statistical Analysis of Differential Protein Expression for Quantitative Proteomics Data. R Package Version 2019,1, 10-18129. 26. Willer, C.J.; Li, Y.; Abecasis, G.R. METAL: Fast and efficient meta-analysis of genomewide association scans. Bioinformatics 2010 , 26, 2190–2191. [CrossRef] [PubMed] 27. Wang, X.; Kang, D.D.; Shen, K.; Song, C.; Lu, S.; Chang, L.C.; Liao, S.G.; Huo, Z.; Tang, S.; Ding, Y.; et al. An R package suite for microarray meta-analysis in quality control, differentially expressed gene analysis and pathway enrichment detection. Bioinformatics 2012,28, 2534–2536. [CrossRef] [PubMed] 28. Kolde, R.; Laur, S.; Adler, P.; Vilo, J. Robust rank aggregation for gene list integration and meta-analysis. Bioinformatics 2012 , 28, 573–580. [CrossRef] [PubMed] 29. Zhang, B.; Kirov, S.; Snoddy, J. WebGestalt: An integrated system for exploring gene sets in various biological contexts. Nucleic Acids Res. 2005,33, W741–W748. [CrossRef] 30. Walter, W.; Sánchez-Cabo, F.; Ricote, M. GOplot: An R package for visually combining expression data with functional analysis. Bioinformatics 2015,31, 2912–2914. [CrossRef] 31. The 1000 Genomes Project Consortium. A global reference for human genetic variation. Nature 2015.526, 68–74. [CrossRef] 32. Clough, E.; Barrett, T. The gene expression omnibus database. In Statistical Genomics; Springer: Berlin/Heidelberg, Germany, 2016; pp. 93–110. 33. Shock, N.W. Normal human aging: The Baltimore longitudinal study of aging. JAMA 1986,255, 960. Appl. Sci. 2022,12, 3987 17 of 17 34. Madrid, L.E.A. Integrated Genomic, Transcriptomic and Proteomic Analysis for Identifying Markers of Alzheimer’s Disease. Diagnosticsl 2021,11, 2303. [CrossRef] 35. de Rojas, I.; Moreno-Grau, S.; Tesi, N.; Grenier-Boley, B.; Andrade, V.; Jansen, I.E.; Pedersen, N.L.; Stringa, N.; Zettergren, A.; Hernández, I.; et al. Common variants in Alzheimer’s disease and risk stratification by polygenic risk scores. Nat. Commun. 2021 , 12, 3417. [CrossRef] 36. Ochoa, D.E.A. Open Targets Platform: Supporting systematic drug–target identification and prioritisation. Nucleic Acids Res. 2021,49, D1302–D1310. [CrossRef] 37. Lambert, J.C.; Ibrahim-Verbaas, C.A.; Harold, D.; Naj, A.C.; Sims, R.; Bellenguez, C.; Jun, G.; DeStefano, A.L.; Bis, J.C.; Beecham, G.W.; et al. Meta-analysis of 74,046 individuals identifies 11 new susceptibility loci for Alzheimer’s disease. Nat. Genet. 2013 ,45, 1452–1458. [CrossRef] [PubMed] 38. Gibbs, R.A.; Belmont, J.W.; Hardenbol, P.; Willis, T.D.; Yu, F.L.; Yang, H.M.; Ch’ang, L.Y.; Huang, W.; Liu, B.; Shen, Y.; et al. The International Hapmap Project; Nature Publishing Group: London, UK, 2003.