Statistical normalisation of network propagation methods for computational biology
Abstract
The advent of high-throughput technologies and their decreasing cost have fostered the creation of a rich ecosystem of public database resources. In an era of affordable data acquisition, the core challenge has shifted to improve data interpretation, in order to understand normal and disease states. To that end, leveraging the current contextual knowledge in the form of annotations and biological networks is a powerful data amplifier to elucidate novel hypotheses. Label propagation and diffusion are the linchpin of the state of the art in network algorithms. In its simplest form, label propagation predicts the labels of a given node (for instance a gene, protein or metabolite) using those of its interactors. More elaborated approaches propagate beyond direct interactors, with robust performance in many computational biology domains. It has been pointed out that the topological structure of biological networks can bias propagation algorithms. Poorly known entities are overlooked and harder to link to experimental findings, which in turn keeps them barely annotated. Some efforts try to break this circularity by statistically normalising the topological bias, but the properties of the bias and the real benefit of its removal are yet to be carefully examined. This thesis covers two blocks. First, a characterisation of the bias in diffusion-based algorithms, with the implementation of statistical normalisations. Second, the application of such normalisation in classical computational biology problems: pathway analysis for metabolomics data and target gene prediction for drug development. In the first block, the presence of the bias is confirmed and linked to the network topology, albeit dependent on which nodes have labels. Equivalences are proven between diffusion processes with variations on their definitions, thus easing its choice. Closed forms on the first and second statistical moments of the null distributions of the diffusion scores are provided and linked to the spectral features of the network. The normalisation can be detrimental if the bias favours nodes with positive labels. An ad-hoc study of the data and the expected properties of the findings is recommended for an optimal choice. To that end, this thesis contributes the diffuStats software package, easing the computation and benchmark of several normalised and unnormalised diffusion scores. The second block starts with pathway analysis for metabolomics data. This choice is driven by the relative lack of computational solutions for metabolomics, whose output still requires an effortful interpretation. Here, a knowledge graph is conceived to connect the metabolites to the biological pathways through intermediate entities, like reactions and enzymes. Given the metabolites of interest, a propagation process is run to prioritise a relevant sub-network, suitable for manual inspection. The statistical normalisation is required due to the network design and properties. The usefulness of this approach is proven not only regarding pathway findings, but also examining the metabolites and reactions within the suggested sub-networks. The knowledge network construction and the propagation algorithm are distributed in the FELLA software package. The second practical application is the prediction of plausible gene targets in disease. Besides benchmarking the effect of the statistical normalisation, particular care is put into obtaining meaningful performance estimates for practical drug development. Target data is usually known at the protein complex level, which leads to performance over-estimation if ignored. Here, this effect is corrected in a varied comparison of prioritisation algorithms, networks, performance metrics and diseases. The results support that the statistical normalisation has a small but negative impact. After correcting for the protein complex structure, network-based algorithms are still deemed useful for drug discovery.
Full text
STATISTICAL NORMALISATION OF NETWORK PROPAGATION METHODS FOR COMPUTATIONAL BIOLOGY author: sergio picart-armada advisor: alexandre perera-lluna A thesis by compendium of publications submitted in fulfilment of the requirements for the degree of Doctor of Philosophy in Biomedical Engineering in the B2SLab Centre de Recerca en Enginyeria Biomèdica Departament d’Enginyeria de Sistemes, Automàtica i Informàtica Industrial Universitat Politècnica de Catalunya February 2020
This document was written with L A TEX using the ArsClassica style by Lorenzo Pantieri ([email protected]), Copyright© 2008-2017
ABSTRACT STATISTICA L NORM ALIS ATION OF NETWORK PRO PAGATION M ETHO DS FOR C OMPU TATION AL BIO LOGY sergio picart-armada The advent of high-throughput technologies and their decreasing cost have fostered the creation of a rich ecosystem of public database resources with molecular annotations and experimental data. In an era of affordable data acquisition and abundant pre-processing tools, the core challenge has shifted to improve data interpretation through algorithms and computational tools. The understanding of normal and disease states is a fundamental piece for generating novel and valuable biological insights. To that end, leveraging the current contextual knowledge in the form of annotations and biological networks can result in a powerful data amplifier and elucidate novel patterns and hypotheses. Label propagation and diffusion are the cornerstone of the state of the art in network mining. They are driven by the guilt by association principle, which states that two interacting biological entities are prone to share functions and properties. In its simplest form, propagation algorithms predict the labels of a given node (for instance a gene, protein or metabolite) using those of its interactors. More elaborated approaches propagate beyond direct interactors, with robust performance in many areas within computational biology. It has been pointed out that the topological structure of biological networks can bias propagation algorithms in such a way that best described entities experience a systemmatic advantage. Poorly known entities are therefore overlooked and harder to link to experimental findings, which in turn keeps them barely annotated. Some efforts try to break this circularity by statistically normalising the topological bias, albeit the properties of the bias and the real benefit of its removal are yet to be carefully examined. The present thesis covers two general blocks. First of all, it seeks a proper characterisation of the bias in diffusion-based algorithms. Statistical normalisations are suggested, implemented and distributed within a scientific software package. The second block covers the application of such normalisation in classical computational biology problems that can be tackled from the network propagation standpoint. In particular, in biological pathway analysis for metabolomics data and in target gene prediction for drug development. In the first block, the presence of the bias is confirmed and linked to the network topology, albeit dependent on which nodes have labels. Some equivalences are proven between diffusion processes with variations on their definitions, therefore easing its choice. Closed forms on the first and second statistical moments of the null distribution of the diffusion scores are provided, with resemblance to the spectral features of the network. Another iii
finding is that the normalisation can be detrimental in certain scenarios, e.g. if the bias favours nodes with positive labels. An ad-hoc study of the data and the expected properties of the findings is recommended for an optimal choice. To that end, this thesis contributes the diffuStats software package to a public repository. diffuStats eases the computation and benchmark of several diffusion scores, including normalised and unnormalised ones. The second block starts with pathway analysis for metabolomics data. This choice is driven by the relative lack of computational solutions for metabolomics for being a younger discipline. The classical over representation analysis starts from a list of metabolites of interest, typically derived from an experimental study, and highlights a list of relevant biological pathways. Newer tools also use metabolic network data in their layout, but the interpretation still entails a demanding manual ad-hoc effort. This block focuses on an enhanced interpretability by building and mining a richer knowledge network. The network connects the metabolites to the biological pathways through intermediate entities, like reactions and enzymes. Given the metabolites of interest, a propagation process is run to prioritise a relevant sub-network, suitable for manual inspection. The statistical normalisation is required due to the network design and properties. The usefulness of this approach is proven not only regarding pathway findings, but also examining the metabolites and reactions within the suggested sub-networks. The knowledge network construction and the propagation algorithm are distributed in the FELLA software package, with six case studies on human and animal datasets. The second practical application is the prediction of plausible gene targets in disease by leveraging biological networks. Besides benchmarking the effect of the statistical normalisation on label propagation, particular care is put into obtaining meaningful performance estimates for practical drug development. Target data is usually known at the protein complex -or even familylevel. Studies that overlook the structure of the protein complex data report overly optimistic performance estimates. In this thesis, this effect is corrected in an exhaustive comparison of prioritisation algorithms, networks, performance metrics and diseases. The results support that the statistical normalisation has a small but negative impact. In broad terms, even after correcting for the protein complex bias, network-based algorithms are still deemed useful and encouraged for drug discovery.
NOR MALI ZACI ÓN ESTADÍST ICA SO BRE LOS ALG ORITMOS DE PROPAGACIÓN EN REDE S PARA BIO LOGÍA C OMPU TA CION AL sergio picart-armada La aparición de tecnologías experimentales de alto rendimiento ha propiciado la creación de un rico entorno de bases de datos que aglomeran todo tipo de anotaciones moleculares. Dada la creciente facilidad para la adquisición de datos en varios niveles moleculares, el reto central de la biología computacional ha virado hacia la interpretación de dicho volumen de datos. La comprensión de los procesos de normalidad y enfermedad involucrados en los cambios observados en los estudios experimentales es el motor que expande la frontera del conocimiento humano. Para ello, es fundamental aprovechar la herencia de conocimiento previo, recogido en las bases de datos en forma de anotaciones y redes biológicas, y minarlo en busca de nuevos patrones e hipótesis. Los algoritmos más extendidos para extraer conocimiento de las redes biológicas son los denominados métodos de propagación y difusión. Su trasfondo es el principio de culpa por asociación, que postula que las entidades biológicas que mantienen relación o interacción son más propensas a compartir funciones y propiedades. Dichos algoritmos aprovechan las interacciones conocidas, en formato de red, para predecir propiedades de nodos (por ejemplo genes, proteínas o metabolitos) usando las propiedades de sus interactores. Existe evidencia de que la estructura topológica de las redes sesga los algoritmos de propagación, de forma que los nodos mejor descritos gozan de una ventaja sistemática. Los nodos menos conocidos quedan en desventaja, se entorpece el descubrimiento de su implicación en los experimentos, a su vez perpetuando nuestro pobre conocimiento sobre ellos. La literatura ofrece algunos estudios donde se normaliza dicho efecto, pero las propiedades intrínsecas del sesgo y el beneficio real de dicha normalización requiere un estudio más detallado. El objeto de esta tesis tiene dos vertientes. Primero, la caracterización de la estadística del sesgo en los algoritmos de propagación, la concepción de normalizaciones estadísticas y su distribución como software científico. Segundo, la aplicación de dicha normalización en problemas clásicos de biología computacional. Concretamente, en el análisis de vías biológicas para datos de metabolómica y en la predicción de genes como dianas terapéuticas en el desarrollo de fármacos. Ambos problemas son abordables mediante técnicas de propagación y, por lo tanto, potencialmente sensibles al efecto del sesgo topológico. En el primer bloque, se corrobora la existencia del sesgo y su dependencia no sólo de la estructura de la red, sino de los nodos en los que se define la propagación. Se demuestran equivalencias matemáticas entre ciertas variaciones en la definición de la propagación, facilitando así su elección. Se proporcionan expresiones cerradas sobre los momentos estadísticos de la difusión y se halla una conexión con las propiedades espectrales de las redes. Un punto importante es que la normalización no siempre ayuda, y su aplicabilidad dependerá de cada caso particular y de las hipótesis sobre la v
topología de los nodos que deben ser descubiertos. Para ello, esta tesis deja como resultado diffuStats, un software disponible en un repositorio púlico, que permite calcular y comparar la propagación con ciertas variantes, y con presencia o ausencia de normalización. En el segundo bloque, se escoge el análisis de vías en metabolómica dada la relativa juventud de los estudios metabolómicos y, por ende, su falta de herramientas informáticas dedicadas. El análisis de vías clásico parte de una lista de metabolitos de interés, normalmente procedentes de un estudio, y reporta una lista de vías o procesos metabólicos estadísticamente relacionados con ellos. Algunas variantes usan redes de metabolitos para dar más contexto biológico, pero la interpretación de los datos sigue requiriendo un extenso esfuerzo manual. La aportación de esta tesis es la creación de una red de conocimiento que relaciona los metabolitos con las vías a través de las entidades intermedias anotadas, como reacciones y enzimas. Sobre dicha red se aplican algoritmos de propagación para identificar las entidades más relacionadas con los metabolitos de interés. La normalización estadística es necesaria, dada la estructura y las características de la red. Se demuestra no sólo la coherencia de las vías metabólicas propuestas, sino la de los metabolitos y las reacciones priorizadas. La publicación del software FELLA proporciona la construcción de la red de conocimiento y el algoritmo de difusión a la comunidad científica. FELLA va acompañado de seis casos de aplicación en estudios humanos y animales. Por otro lado, se aborda el problema de predicción de genes para dianas terapéuticas a través de redes biológicas. Además de probar el efecto de la normalización estadística, se pone énfasis en estimar el desempeño real esperado en un escenario de desarrollo de fármacos. Los datos de dianas terapéuticas no se suelen conocer al nivel de proteína sino al de complejo o familia de proteínas. La mayoría de estudios no lo tiene en cuenta, llegando a estimaciones optimistas sobre el desempeño esperado. En esta tesis se propone un estudio exhaustivo que corrige el efecto de los complejos de proteínas, compara algoritmos de propagación con distintas métricas de rendimiento por su informatividad y explora el rol de la red biológica y de la enfermedad en cuestión. Se demuestra que la normalización estadística tiene poco efecto en el desempeño y que, en general, los métodos de propagación siguen siendo útiles en el desarrollo de fármacos después de corregir las estimaciones optimistas de su rendimiento.
ACKNOWLEDGEMENTS On the scientific principle behind this thesis: La teva tesi és de trileros — Alexandre Perera Lluna On how to properly manage a scientific project: la cencia no se ace sola ahi que acerla — @cientefico On how (not) to compromise one’s mental health: Your closest collaborator is you six months ago, but you don’t reply to emails. — Karl Broman This doctoral thesis contains many fingerprints, one for every person who helped, supported or beared me throughout the process. Actually, here are some numbers: this is the second time I write this section; my first attempt led to 5,757 characters, forming 978 words in 68 sentences that thanked exactly 107 people. As it felt a bit overwhelming to read, it was heavily summarised, but I must express my sincerest gratitude to all of them. Thank you. I feel fortunate for having been a part of the B2SLab in Barcelona under the guidance of Alex Perera. You are gifted both within and outside the scientific scope, coupled with a selfless, educational, caring and easy-going personality. You know, you could probably make a living out of your coaching skills! You have shaped my critical thinking into the scientist I am now. Thank you for your academic and personal support throughout these years. My deepest thanks to Francesc too, who can be seen as my older scientific brother, as we share our scientific father. Your support has been priceless and even now I struggle to find the right words to thank you properly. Despite your departure from the B2SLab (a sad moment for me), you kept in touch and eventually arranged an enriching internship in Takeda. I look forward to keep learning from you, as a professional and as a friend, and I wish you the very best as a scientist and as a dad. Another key person is Alfonso Buil, my supervisor in my recurrent scientific stays in the Sct. Hans institute, Denmark. You have been outstandingly hospitable and attentive. You taught me statistics, genetics, but especially to keep calm and carry on, owing to your smooth attitude. Off topic: one could easily add a joke about bones in the previous sentence. Anyway, thank you for making it possible. vii
viii One name that deserves a special place is Pere Caminal. Not only for being one of the fathers of bioengineering research in Spain, but for your unconditional kindness and personal help. You brought us together and still show genuine interest for each one of us. Thank you. Warm regards to the B2SLab, both the old and the young generations, for creating a cozy working environment. I wish the best to Pol S. We had an interesting time working together and our paths may encounter again in the future. A special mention to Josep, who took over some of my personal projects with enthusiasm. I had the opportunity to visit external research centers and enjoyed an enriching experience. I want to thank Oscar Yanes and his lab for hosting me in Reus and having the perseverance to guide me on the first (and bumpy) stages of my research career. My warmest regards to the people in Sct. Hans, Denmark. Thank you Xabi for being a good friend and integrating me in your social circles. Thank you Camellia for being a warm neighbour and for your exceptional indian cuisine. Thank you Wes for listening to boring and convoluted statistical stuff. I would also like to acknowledge the people involved in the GAIT project. I also keep good memories of the Takeda people. Special thanks to Ben for keeping in touch and leading a fruitful collaboration with GSK. To my University friends, highlighting Guillem B. for regularly checking on me and sharing our thoughts between beer and beer, Pablo, for our long calls and your caring attitude, and to Ferran, for the good old times. To Stack Overflow and Wikipedia, no need to explain. To Sandwichez, where I spent too many hours (and money). To the Penyafort team and to the Galileu community, led by Jas, for our fond memories. To Alvaro N., for being a good listener, preparing amazing food and keeping a relaxed atmosphere with his intelligent humour. To my closest people, especially to Xavier S. and Guillem M., who have accompanied me since high school. To my parents Manuel and Pilar, who did not hesitate for a moment about my capabilities. Thank you for looking after me and for your genuine interest in what I do, even if most of it looks weird and hard to understand. To my family, especially to my grandparents. To Angela, probably the one I annoyed the most, but with whom I shared the best moments too. I owe my mental sanity to you. You just knew how to take my mind out of vicious circles. You helped me keeping a work-life balance, putting everything in perspective and understanding the value of certain things. Thank you for being there for me. Wrapping up this chapter in my life feels somewhere between unreal and relief. I will miss some aspects of the PhD life. Time to move on. Thank you.
CONTENTS 1 introduction 1 1.1Omics sciences 1 1.1.1Introduction 1 1.1.2Genomics 1 1.1.3Transcriptomics 2 1.1.4Proteomics 3 1.1.5Metabolomics 3 1.1.6Other omics 3 1.2Data interpretation 4 2 state of the art 9 2.1Network representations in biology 9 2.1.1Database resources 10 2.1.2Specialised databases 11 2.1.3Network databases 14 2.1.4Comprehensive databases 17 2.2Network propagation algorithms 21 2.2.1Introduction to network propagation 21 2.2.2Introduction to graph theory 22 2.2.3Network propagation in computational biology 25 2.2.4Statistical properties of network propagation 29 2.3Applications of network propagation 33 2.3.1Metabolomics data enrichment 33 2.3.2Disease gene identification 34 2.4Open issues 35 2.4.1Heterogeneity and biases in network propagation 35 2.4.2Results interpretability in pathway analysis 36 2.4.3Performance overestimation in target gene prediction 36 2.4.4Free and open source software 36 3 goals 51 3.1Main objective 51 3.2Detailed objectives 51 3.2.1Conception of the statistical normalisation 51 3.2.2Application to metabolomics data enrichment 51 3.2.3Application to gene target discovery 51 3.3Expected contributions 52 4 statistical properties 53 4.1Introduction 53 4.2Approach 55 4.3Methods 55 4.3.1Unnormalised scores 55 4.3.2Normalised scores 56 4.3.3Metrics and baselines 56 4.3.4Bias quantification 57 4.3.5Performance explanatory models 57 ix
2introduction ...ATTGCGCCGAGATATTTAGCGGATGCAAATAAGATTTATTACC... Transcription factor Gene AUAUUUAG GAUGCAAA AAGAUUUA Messenger RNA Protein products A B C A B C D E F G H JK L K M Genomics Transcriptomics Proteomics Metabolomics ATP ADP Genotype Phenotype Figure 1: Overview of omics sciences. Genomics, trancriptomics, proteomics and metabolomics are the main omics sciences, whose measurements range from genotypic to phenotypic data. Figure adapted from ‘Figure 1’ in (Joyce and Palsson,2006). understanding of numerous complex traits through their findings (McCarthy et al.,2008). Currently, whole genome and whole exome (protein-coding regions) sequencing are improving our ability to discover genetic variants in human populations (Petersen et al.,2017). 1.1.3 Transcriptomics Transcriptomics measures the presence and abundance of RNA transcripts (Joyce and Palsson,2006). The total amount of messenger RNA in a cell or organism is called transcriptome. The assessment of differential gene expression was pioneering in the study of disease in the late 1990s and has generated a remarkable amount of biological knowledge. RNA-seq is a standard
1.1 omics sciences 3 technology to measure gene expression (Wang et al.,2009) and disposes of a rich array of tools to analyse its data. Despite being a valuable source of information, differences in transcript abundances do not necessarily imply the same changes at the protein level (Joyce and Palsson,2006), meaning that transcriptomics alone is unable to explain the whole cellular state. 1.1.4 Proteomics Proteomics focuses in identifing and quantifying the proteins within cells and tissues (Joyce and Palsson,2006). The proteome is the set of all the expressed proteins in a cell or organism. Proteins orchestrate the metabolism, albeit their activity is in turn affected by the metabolic state of the cells. The proteome is a dynamic reflection of the combination of genetic and environmental factors and is considered an excellent source of disease biomarkers (Horgan and Kenny,2011). Likewise, interaction events between proteins have been thoroughly studied and proven to be useful for applying networkbased algorithms (Cowen et al.,2017). 1.1.5 Metabolomics Metabolomics is the study of metabolites, the lightweight molecules that can be found within living organisms. The collection of metabolites found within cells or organisms is called the metabolome, also including those coming from the environment. Metabolite measurements are, in fact, quantitative phenotypes that give a snapshot of the functional readout of the cells (Joyce and Palsson,2006). This is particularly appealing since it displays the actual effect of genomic or transcriptomic events. Compared with transcriptomics and proteomics, metabolomics data poses further statistical challenges due to its technical limitations (Joyce and Palsson,2006), its physical and chemical complexity (Horgan and Kenny,2011) and the unknown extent of the human metabolome (Wishart et al.,2012). 1.1.6 Other omics Besides genomics, transcriptomics, proteomics and metabolomics, other omics sciences are emerging in terms of experimental techniques, public data and computational tools (Joyce and Palsson,2006). Metagenomics is the analysis of genetic material from environmental samples, typically involving microbial communities. Epigenomics measures epigenetic events on the genetic material of cells, such as histone modification, DNA methylation and chromatin accessibility assays. The study of microRNA data is sometimes referred to as miRNomics.Lipidomics denotes the study of lipids, whereas glycomics revolves around carbohydrates and glycans. This list is not exhaustive but illustrative on the richness of current data acquisition and generation capabilities.
4introduction 1.2 data interpretation Leveraging omics data has provided key biological insights about normal and disease states, novel drug targets, drug response, biomarkers and predictive models for diagnosis and prognosis (Horgan and Kenny,2011). The volume of high-throughput data generated in omics sciences has yielded high dimensional data that requires careful statistical treatment (Horgan and Kenny,2011) and is challenging to interpret and understand (Joyce and Palsson,2006). This issue first appeared with the advent of microarrays: a formal approach was needed to contextualise experimental results, usually extensive lists of differentially expressed genes. The solution was named functional analysis and relied on classical statistical tests to assess if any known molecular function appeared more than expected within the gene list. Grouping genes that functioned in the same biological processes and finding dysregulated processes reduced the complexity of the data while providing richer mechanistical insights (Khatri et al.,2012). This is still a simple yet powerful approach, usually referred as over representation analysis, allowing the test of virtually any known annotation type. It plays a prominent role in translating differentially abundant genes, proteins or other molecular entities stemming from high-throughput technologies into biological knowledge (Mitrea et al.,2013), as illustrated in figure 2. (B) Data interpretation (A) Experimental data G1 G2 G3 G4 G5 G6 G7 G8 G9 G10 WT1 KO1 WT2 WT3 KO2 KO3 Technical platform Preprocessed data Differential abundance F1: G1, G2, G3, G4, G10 F2: G4, G5, G6 F3: G1, G7, G8, G9, G10 G1** G2*** G3*** G4** G5 G6 G7 G8· G9 G10 WT vs KO Functional annotations -log(p-value) F1 F2 F3 G1 G2 G3 G4 Genes of interest Over represented functions Figure 2: Examplary omics data analysis. This small example illustrates a hypothetical transcriptomics workflow in a case-control experiment. (A) Measurement of gene expression and discovery of differentially expressed genes between wild type (WT) and knockout (KO) experimental groups. (B) Discovery of highly occurring functions among the differentially expressed genes.
1.2 data interpretation 5 Data interpretation approaches have been extended to work on quantitative data and leverage biological network databases (Khatri et al.,2012), but the interpretability of their outputs is still an area of active research. The integration of different omics, which provide complementary views of a common reality, is a promising approach to attain a holistic picture of the molecular processes underlying disease (Ge et al.,2003).
6introduction references Cowen, Lenore, Trey Ideker, Benjamin J Raphael, and Roded Sharan 2017 “Network propagation: a universal amplifier of genetic associations”, Nature Reviews Genetics,18,9, p. 551. Ge, Hui, Albertha JM Walhout, and Marc Vidal 2003 “Integrating ‘omic’ information: a bridge between genomics and systems biology”, TRENDS in Genetics,19,10, pp. 551-560. Gomez-Cabrero, David, Imad Abugessaisa, Dieter Maier, Andrew Teschendorff, Matthias Merkenschlager, Andreas Gisel, Esteban Ballestar, Erik Bongcam-Rudloff, Ana Conesa, and Jesper Tegnér 2014 Data integration in the era of omics: current and future challenges. Goodwin, Sara, John D McPherson, and W Richard McCombie 2016 “Coming of age: ten years of next-generation sequencing technologies”, Nature Reviews Genetics,17,6, p. 333. Horgan, Richard P and Louise C Kenny 2011 “‘Omic’ technologies: genomics, transcriptomics, proteomics and metabolomics”, The Obstetrician & Gynaecologist,13,3, pp. 189-195. Joyce, Andrew R and Bernhard Ø Palsson 2006 “The model organism as a system: integrating ’omics’ data sets”, Nature reviews Molecular cell biology,7,3, p. 198. Khatri, Purvesh, Marina Sirota, and Atul J Butte 2012 “Ten years of pathway analysis: current approaches and outstanding challenges”, PLoS computational biology,8,2, e1002375. McCarthy, Mark I, Gonçalo R Abecasis, Lon R Cardon, David B Goldstein, Julian Little, John PA Ioannidis, and Joel N Hirschhorn 2008 “Genome-wide association studies for complex traits: consensus, uncertainty and challenges”, Nature reviews genetics,9,5, p. 356. Mitrea, Cristina, Zeinab Taghavi, Behzad Bokanizad, Samer Hanoudi, Rebecca Tagett, Michele Donato, Calin Voichita, and Sorin Draghici 2013 “Methods and approaches in the topology-based analysis of biological pathways”, Frontiers in physiology,4, p. 278. Petersen, Britt-Sabina, Broder Fredrich, Marc P Hoeppner, David Ellinghaus, and Andre Franke 2017 “Opportunities and challenges of whole-genome and-exome sequencing”, BMC genetics,18,1, p. 14. Sanger, Frederick, Steven Nicklen, and Alan R Coulson 1977 “DNA sequencing with chain-terminating inhibitors”, Proceedings of the national academy of sciences,74,12, pp. 5463-5467.
References 7 Venter, J Craig, Mark D Adams, Eugene W Myers, Peter W Li, Richard J Mural, Granger G Sutton, Hamilton O Smith, Mark Yandell, Cheryl A Evans, Robert A Holt, et al. 2001 “The sequence of the human genome”, science,291,5507, pp. 1304- 1351. Voelkerding, Karl V, Shale A Dames, and Jacob D Durtschi 2009 “Next-generation sequencing: from basic research to diagnostics”, Clinical chemistry,55,4, pp. 641-658. Wang, Zhong, Mark Gerstein, and Michael Snyder 2009 “RNA-Seq: a revolutionary tool for transcriptomics”, Nature reviews genetics,10,1, p. 57. Wishart, David S, Timothy Jewison, An Chi Guo, Michael Wilson, Craig Knox, Yifeng Liu, Yannick Djoumbou, Rupasri Mandal, Farid Aziat, Edison Dong, et al. 2012 “HMDB 3.0–the human metabolome database in 2013”, Nucleic acids research,41, D1, pp. D801-D807.
2STATE OF T HE ART This chapter covers the concept of annotation databases, biological networks and propagation algorithms on them. Figure 3pictures how contextual data, in network format, enriches experimental data by enabling predictions on a variety of domains. It also points out the logical order of the sections, whose topics are biological networks (data origin, definition and construction), propagation algorithms (graph theory definitions, the guiltby-association principle, algorithm formulations) and two case studies (pathway analysis and disease gene prediction). Experimental data Network propagation Novel predictions 1. Network representation Network data 2. Network-based algorithms 3. Applications Figure 3: Overview of network propagation. Conceptual map of how network propagation takes advantage of network data to bring new insights from experimental data. The sections cover (1) the construction of biological networks from public data, (2) the definition and uses of network propagation algorithms and (3) two domains for the application of network propagation. 2.1 network representations in biology Network data is a central concept in computational biology, both as a way to represent current knowledge and a corpus that enables new predictions. A proper knowledge representation is essential to provide sensible predictions through network propagation algorithms. Figure 4displays these ideas, covered in this section. The logical order of the sections is as follows: section 2.1.1introduces large projects that contribute abundant data to the public domain, section 2.1.2covers specialised databases that annotate specific molecular levels, and both typically act as building blocks for network and 9
10 state of the art pathway databases (sections 2.1.3and 2.1.4). Network-based algorithms can be applied to network data and to pathway annotations, provided that the latter are represented as networks. Growing public data Network DBs Biological networks Specialised DBs Comprehensive DBs Novel predictions Figure 4: Overview of network representations. Conceptual map of public data is shaped into different types of database (DB), which in turn can be used to build biological networks. Network-based algorithms can be applied to biological networks to generate new knowledge. 2.1.1 Database resources National and international consortia regularly foster large scale studies, which aim to translate large sample sizes into meaningful knowledge. This section contains a list of varied initiatives, to highlight the outstanding value of large scale data recollection. One prominent example is the Encode project (ENCODE Project Consortium and others,2012), aimed at annotating all the regions of the human genome throught the intregration of thousands of datasets. Likewise, the 1000 Genomes project (1000 Genomes Project Consortium and others,2015) reconstructed 2,504 genomes from 26 populations in order to generate a high-quality reference panel of human genetic variation, annotating over 88 million variants. Another reference panel is provided by the UK10K project (UK10K consortium and others,2015). Owing to the recording of several phenotypes, they further assess the contribution of genomic variation to a set of traits and causal mutations for disease. In the topic of mental disorders, the iPSYCH cohort (C. B. Pedersen et al.,2018) aims at finding novel genetic and environmental factors of conditions like schizophrenia, autism, attention-deficit/hyperactivity disorder or bipolar disorder. iPSYCH disposes of a Danish population-wide registry, granting an important statistical advantage. The Genotype-Tissue Expression (GTEx) project (Lonsdale et al.,2013) is a systematic study on the effect of genetic variations on gene expression in
2.1 network representations in biology 11 human tissues. Post-mortem samples were collected and analysed using genomics (whole genome sequencing) and transcriptomics (RNA sequencing). Scanned images and clinical data are also available. The incipient Human Cell Atlas (Regev et al.,2017) will explore cell types exhaustively with single cell technologies. Regarding oncology, The Cancer Genome Atlas (Tomczak et al.,2015) is a landmark cancer genomics program that characterised primary cancer and matched normal samples from 11,000 patients and 33 cancer types1. Genomic, epigenomic, transcriptomic, proteomic and clinical data was leveraged to publish marker manuscripts for each cancer type2 and to elucidate key commonalities and differences across cancer types and tissues (Hoadley et al.,2014). The Connectivity Map, or CMap (Subramanian, Narayan, et al.,2017), aims at understanding cellular function by a large scale experiment on human cell lines with an array of perturbagens. Gene expression was quantified through L1000, a novel cost-effective profiling method that directly measures 978 landmark genes and accurately imputes 9,196 genes out of the 11,350 remaining transcripts. A total of 1,319,138 L1000 profiles are available and consolidated into 473,647 signatures that involve up to 77 cell lines. Along with data from large initiatives, the number of scientific articles grows steadily too. Open science policies and protocols are gaining traction and encouraging data deposition on public online repositories with common protocols, such as Gene Expression Omnibus, or GEO (Clough and Barrett, 2016) and MetaboLights (Kale et al.,2016). Likewise, platforms like GitHub3, Zenodo4and figshare5facilitate the storage of computer code and data. Combining the increasing data and knowledge availability, there has been a need of annotation resources that centralise, format and curate data from the public domain. Two examples are the GWAS catalog (MacArthur et al., 2016) for genetic associations and the Expression Atlas (Papatheodorou et al., 2017) for transcriptomics studies. These efforts ease large-scale analyses and meta-studies, which leverage large sample sizes to unravel novel biological insights. 2.1.2 Specialised databases A plethora of databases are publicly available, essentially at any molecular level ranging from genomic to phenotypic data. Table 1displays a selection of such resources, some of which were already mentioned in section 2.1.1. Those cover, in order: genes, proteins, metabolites, compounds and phenotypes. The Gene Ontology, or GO (G. O. Consortium,2016) terms are a common choice to annotate gene function in three aspects: molecular function, cellular component and biological process. GO terms conform a hierarchy with 1https://www.cancer.gov/about-nci/organization/ccg/research/structural-genomics/ tcga/history. Accessed on 31/12/2019. 2https://www.cancer.gov/about-nci/organization/ccg/research/structural-genomics/ tcga/publications. Accessed on 31/12/2019. 3https://github.com/. Accessed on 26/1/2020. 4https://zenodo.org/. Accessed on 26/1/2020. 5https://figshare.com/. Accessed on 26/1/2020.
18 state of the art correlated) do not necessarily imply physical interaction or genetic regulation. Table 3displays a selection of pathway databases, further elaborated in this section and arranged according to their principal subject: genetic events, metabolism, signalling, disease and integrative. Table 3: Selection of public pathway resources, sorted by their main focus: genetic, metabolic, general or integrative. Resource name Main subject Reference PANTHER Genetics (Mi et al.,2018) SMPDB Metabolic (Jewison et al.,2013) LIPID MAPS Metabolic, lipids (V. B. O’Donnell et al.,2019) SwissLipids Metabolic, lipids (Aimo et al.,2015) MetaCyc signalling, metabolic (Caspi et al.,2019) KEGG signalling, metabolic, disease (Kanehisa et al.,2016) Reactome signalling, metabolic, disease (Fabregat et al.,2017) WikiPathways signalling, metabolic, disease (Slenter et al.,2017) PathBank signalling, metabolic, disease (Wishart, Carin Li, et al.,2019) Pathway Commons Integrative: pathways, interactions (Rodchenkov et al.,2019) ConsensusPathDB Integrative: pathways, interactions (Herwig et al.,2016) BioModels Integrative: models (Malik-Sheriff et al.,2019) The PANTHER database (Protein ANalysis THrough Evolutionary Relationships) (Mi et al.,2018) provides evolutionary and functional annotations for genes in over 900 genomes. Besides nourishing from Gene Ontology (G. O. Consortium,2016), PANTHER contains a collection of 177 pathways with 3,092 pathway components, 53,548 associated sequences and capturing 6,000 references (PANTHER™ Pathway 3.6.3, released December 2019). Conversely, the SMPDB (Small Molecule Pathway Database) (Jewison et al.,2013) focuses in human pathways, for which small molecules play a central role. Its 2.0version includes the following pathway types: metabolic (92), disease (221), drug action (232), drug metabolism (53), physiological action (5) and small molecule signalling (15). In most of them, the cellular location, tissue or organ where reactions take place are available. Lipids encompass a fundamental part of the metabolism and contribute to its understanding, but the technical limitations hinder their characterisation (V. B. O’Donnell et al.,2019). The LIPID MAPS (Lipid Metabolites and Pathways Strategy) initiative (V. B. O’Donnell et al.,2019) categorised over 30,000 lipids from several organisms (Hartler,2015) and contributed with 10 pathways, including the metabolism of cholesterol, eicosanoids, glycerolipids, omega fatty acids and sphingolipids8. Alternatively, SwissLipids (Aimo et al.,2015) features an in silico library of 244,155 feasible lipid structures, more than 2,000 curated enzymatic reactions linking to over 800 proteins and involving glycerophospholipids, glycerolipids, sphingolipids, sterols, fatty acids, fatty alcohols and wax esters (Aimo et al.,2015). Certain pathway databases aim at a broader understanding of biological pathways, usually within multiple organisms, in the context of metabolism, signalling events, genetic regulation and disease states. They also rely on 8http://www.lipidmaps.org/resources/pathways/index.php. Accessed 18/12/2019.
2.1 network representations in biology 19 and link to multiple specialised databases for entities like sequences, proteins, enzymes, metabolites, lipids, drugs and diseases. MetaCyc (Caspi et al.,2019), the Kyoto Encyclopedia of Genes and Genomes (KEGG) (Kanehisa et al.,2016), Reactome (Fabregat et al.,2017) and WikiPathways (Slenter et al.,2017) are widespread public resources for this purpose. Commercial options include Ingenuity Pathway Analysis9and MetaCore10 but are out of the scope of this thesis. MetaCyc offers over 2,570 pathways based on experimental evidence from more than 54,000 publications, and involving 14,003 compounds and 15,691 reactions in version 21.1(Caspi et al.,2019). Although MetaCyc mainly covers small molecule metabolism, the amount of data on macromolecular mecanism is increasing. Specifically, 35 species are annotated with 20 or more pathways (294 in Homo sapiens). By August 2017, almost 11,000 organism-specific semi-automatic metabolic networks were described in pathway/genome databases. KEGG also offers curated organismal pathways, ranging from metabolic to signalling processes. KEGG contains the following database categories: systems information for pathway data, genomic information, chemical information for metabolites, KEGG LIGAND for reactions and enzymes, health information for diseases and KEGG MEDICUS for drugs (Kanehisa et al., 2016). As of October 2016, KEGG encompassed 496 manually drawn pathway reference maps. One example can be found in figure 7: the photosynthesis KEGG pathway, depicting its enzymatic reactions, metabolites and related pathways. Reactome is a knowledge representation that describes human signal transduction, transport, DNA replication, metabolism in a single, consistent data structure (Fabregat et al.,2017). Reactome version 62 covers 10,719 genes, 24,704 protein forms (including post-translational modifications and cellular localisations), 1,768 metabolites and 11,302 reactions drawn from 27,526 scientific articles. 2,012 human pathways are divided into 26 superpathways that represent broad biological domains. Disease annotations are present in the form of 906 disease-specific reactions, annotated from 1,334 mutated variants in 285 gene products. Wikipathways, on the other hand, is based on a crowdsourcing paradigm to annotate biological pathway data (Slenter et al.,2017). WikiPathways focuses on 25 reference species, annotating 2,614 pathways that were contributed by 634 individuals by September 2017. These involve 11,532 genes (7,982 related to metabolic processes) and 3,133 metabolites, with an increasing effort to improve the coverage of the latter. PathBank (Wishart, Carin Li, et al.,2019) is a recent effort on 10 model organisms to provide a pathway for every protein and a map for every metabolite. Over 110,000 pathways containing 78,488 compounds (including metabolites and drugs), 8,993 proteins and 176,535 reactions/interactions are available for metabolism, signalling, disease, drugs and physiology. Their metadata provide subcellular locations, cofactors and protein quaternary structures. PathBank aims at improving the coverage of lipid syn- 9https://www.qiagenbioinformatics.com/products/ingenuity-pathway-analysis/. Accessed 26/01/2020. 10 https://portal.genego.com/. Accessed 26/01/2020.
20 state of the art thesis, metabolite signalling, small molecule hormone signalling and small molecule drug action in KEGG and MetaCyc, while keeping the coverage of protein and cellular signalling from Reactome and Wikipathways. Integrative efforts have emerged to ease searching, downloading, querying, browsing and analysing a collection of varied pathway databases. Pathway Commons (Rodchenkov et al.,2019) aggregates 22 public databases (including PANTHER, KEGG, Reactome and WikiPathways) into 4,794 human pathways (18,490 genes and 11,437 metabolites) and 2.3milion interactions by February 2019. Likewise, ConsensusPathDB (Herwig et al.,2016) integrates 32 databases for human, 15 for mouse and 14 for yeast. In humans, 158,523 physical entities are annotated with 458,570 interactions and conform 4,593 pathway gene sets. ConsensusPathDB further allows tasks such as pathway analysis, heterogeneous network inference and module analysis, starting from genome-wide data or priority lists of genes, proteins or metabolites. BioModels (Malik-Sheriff et al.,2019) is based on a systems biology approach, storing about 2,000 mathematical models from the literature. Such models can predict the states of biological systems, ease the elaboration of novel hypotheses and improve our mechanistic understanding. Models on cell signalling, metabolic pathways and gene regulation are also contextualised with cross-references to standard data resources using machine-friendly controlled vocabularies. Figure 7: The photosynthesis KEGG pathway (identifier: map00195). Metabolites are represented by circles, reactions by arrows and their enzymes by superposed rectangles. Other neighbouring pathways are visible, like Carbon fixation in photosynthetic organisms. Note how also gene/protein names for ortholog groups are provided. Image downloaded from https://www.genome.jp by 9/1/2019. Despite is usefulness, leveraging biological pathway data suffers from limitations. Pathways are under continuous construction and considered to
2.2 network propagation algorithms 21 be highly incomplete (Ogris et al.,2016), limiting the statistical power of pathway-based approaches. Maintaining pathway databases is demanding, leading to an inability to keep up with the growing liteature or even to discontinuation (Rodchenkov et al.,2019). Another common critique arises from the fact that manual curation results in artificial borders between biological pathways. Consequently, notable differences exist between major pathway databases, in terms of focus, coverage, granularity and pathway definition (Domingo-Fernández et al.,2018), implying that the database choice has a considerable impact in any downstream analysis. In this context, the emergence of integrative resources provide a proxy to alleviate databaserelated biases in data mining. However, the degree of overlap, complementarity, pathway cross-talk and even disagreement needs careful investigation (Domingo-Fernandez et al.,2019). The PathMe platform (Domingo- Fernandez et al.,2019) is a pioneering effort to harmonise and understand the merging of the human data in KEGG, Reactome and WikiPathways under a common controlled vocabuary. Considering the wide spectrum of databases and organisms, there is still room to scale up pathway database harmonisation. 2.2 network propagation algorithms Once network data is available, network propagation allows the integration of experimental (or annotated) data with contextual knowledge. This section, conceptualised in figure 8, covers the mathematical definition of the networks, the algorithms applied on them (through specific examples in computational biology) and the examination of their statistical properties. 2.2.1 Introduction to network propagation High-throughput techniques are contributing in the prediction and identification of a vast collection of molecular interactions. This has led to a rich variety of biological network resources (section 2.1.3), such as proteinprotein interaction, gene regulatory, co-expression and metabolic networks (J. K. Huang et al.,2018). Such networks are usually defined as a set of nodes and a set of edges that connect pairs of nodes. For example, nodes can be proteins and edges can be experimentally proven interactions. Both nodes and edges can have additional attributes, like edge directionality and weight. On the other hand, the guilt by association principle states that interacting entities are more prone to share molecular functions. This basic concept has found ubiquitous application to problems like protein function prediction, disease gene prioritisation (figure 2B), module inference, cancer patients stratification, drug discovery and causal variants identification (Cowen et al.,2017). The paradigm behind network propagation is to infer the labels of molecular entities using the neighbouring connections from a biological network and a set of known, labelled entities. The simplest approach, neighbour vot-
22 state of the art Guilt by association Diffusion models Statistical properties Applications in computational biology Random walks Kernel methods Definitions Topological bias Figure 8: Overview of network propagation algorithms. Conceptual map of the section layout. Relying on basic network definitions, the guilt by association principle justifies the concept of network propagation, which accepts a variety of formulations based on random walks, heat (or another abstract entity) diffusion or graph kernels. The statistical description of the predictions of such methods raised concerns about the presence of a topological bias, i.e. related to the intrinsic properties of the networks. ing, infers the label of an entity from the known labels of its neighbours. A use case could be as follows: to infer whether a protein is related to the obesity phenotype, one counts the proportion of obesity-related interacting proteins. Proteins with the highest proportions are suggested as potential associations. More sophisticated network propagation or diffusion approaches allow the propagation to reach beyond direct neighbours and attain competitive performances in many computational biology applications (Cowen et al.,2017). Following the case study of obesity, one can conceive an abstract substance (heat, fluid, current) that flows from the obesity-related proteins to the rest of the network through the edges. The degree of association of each protein is then measured by the substance received at every protein. In this regard, the diversity of problem formulations and computational biology applications is notable, sometimes causing the re-discovery of equivalent methods under different names and domains. The term “network propagation” therefore stands for a general purpose, heterogeneous but unifying formalism for network analysis. 2.2.2 Introduction to graph theory This section covers basic notions on graph theory prior to the introduction of propagation methods in section 2.2.3. These definitions will be referenced throughout this thesis.
2.2 network propagation algorithms 23 Defining a graph The mathematical definition of a network or graph G, adapted from (Smola and Kondor,2003), consists of two sets: G= (V,E)(1) Vis the set of vertices, or nodes, typically numbered from 1to n. The number nof nodes in the graph is called the graph order.Eis the set of edges, consisting of pairs of nodes (i,j)that indicate that there is a connection from ito j, meaning that they are neighbours. Here, only finite graphs without multiple edges (i.e. ‘repeated’) or loops (edges of the kind (i,i)) are considered. G0= (V0,E0)is a subgraph of Gif G0is a graph with V0⊆Vand E0⊆E. If Gis weighted, each edge (i,j)comes with a weight Wij ∈R. The methods hereby discussed assume Wij ⩾0: the greater the weight, the easier it is to traverse the edge. If Gis unweighted,Wij =1if (i,j)∈Eand Wij =0 otherwise. If Gis undirected, the edges (i,j)and (j,i)are identical and usually denoted as unordered pairs {i,j}, and Wij =Wji. Conversely, in a directed graph, the edge (i,j)can only be traversed from ito jand does not imply the existence of the edge (j,i). Asimple graph is an unweighted, undirected graph without multiple edges or loops (Diestel,2000). The n×nreal matrix Wis called the adjacency matrix of G. The degree matrix of an undirected graph Gis the diagonal matrix Dwith Dii =Pn j=1Wij. Dii is called the degree of vertex i. Note that if Gis directed, one can either use the in-degree Dii =Pn j=1Wji or the out-degree Dii =Pn j=1Wij (Bang- Jensen and Gutin,2008). If Gis simple, Dii is the number of neighbours of the i-th node. Nodes with a high amount of neighbours are called hubs (Cowen et al.,2017), whereas nodes with no connections are isolated. Examples can be found in figure 9. Walks, paths and connectivity This section is based on the definitions in (Diestel,2000) for simple graphs. Awalk in a graph G= (V,E)is a sequence of alternating nodes vi∈Vand edges ej∈E, represented as v1,e1,v2,e2,...,ek−1,vk, starting and ending on nodes, with el= (vl,vl+1), i.e. every edge connects the nodes before and after it. A path is a walk without repeating nodes and its length is the number of edges in it. Ashortest path between two nodes viand vjis a path whose length is minimum among those that start at viand end at vj. If such a path exists, its length is called the shortest path distance dG(vi,vj)between viand vj, which are denoted as connected. Otherwise, if no path exists between viand vj, then dG(vi,vj) = ∞by convention and they are disconnected. A simple graph Gis connected if every possible pair of its nodes is connected. A connected component of a graph Gis a maximal connected subgraph of G.
24 state of the art Undirected, unweighted Directed Multiedge Loop 0.7 Weighted (A) Types of edges (B) Network definition B E A G C I D J H F (C) Network plot (D) Laplacian matrix A B C D E F G H I J J I H G F E D C B A (E) Diffusion kernel - regul. Laplacian A B C D E F G H I J J I H G F E D C B A (F) Random walk matrix P Figure 9: Examples on network definitions. Each sub-figure represents one stage in the network definition and analysis. (A) A network can be built using different types of edges. (B) Definition of the vertex and edge sets. (C) Plot of the network. (D) Unnormalised graph Laplacian matrix. The diagonal contains the node degree. Node F would be the equivalent of a biological hub, given its high degree. (E) Kernel matrix for label propagation, derived from the Laplacian matrix. Darker colours reflect a higher node similarity. Equation 14 details its formal definition, with σ2=1 and using Linstead of ˜ L.(F) Random walk matrix Pfrom equation 4. The darker the colour, the higher the probability of transitioning from the row-indexed to the column-indexed vertex. Note how Pis asymmetric. The graph Laplacian matrix The starting point of many network-based propagation methods is the graph Laplacian matrix, which comes in two flavours: the unnormalised graph Laplacian Land its normalised version ˜ L. These are defined for simple graphs but naturally work on weighted graphs (Smola and Kondor,2003):
2.2 network propagation algorithms 25 L:= D−W(2) ˜ L:= D−1 2LD−1 2=In−D−1 2WD−1 2, (3) where Inis the n×nidentity matrix, being nthe order of the graph. Isolated nodes typically require a special treatment in ˜ Lbecause their degree is 0. Figure 9D contains a small example on how to compute L. Because Gis undirected, Land ˜ Lare symmetric and therefore diagonalisable. The spectral properties of Land ˜ Lhave been extensively studied and are tightly connected to the topological properties of G. Random walks A notable branch of graph theory is the study of random walks, which are markov chains on graphs (Lovász,1993). In each step of a random walk, a fictional walker sits on a node iand randomly chooses an edge to resume her random walk, or decides to start another walk. The stationary probability distributions of such processes exist and have been extensively studied. The matrices in equations 4and 5are normalised versions of the adjacency matrix, commonly used to compute random walk-based scores on undirected graphs (Cowen et al.,2017). Again, isolated nodes require special treatment. An example of the former can be found in figure 9F. P:= WD−1(4) ˜ P:= D−1 2WD−1 2(5) 2.2.3 Network propagation in computational biology The starting point of network-based algorithms has its roots in the Guilt By Association (GBA) principle (Oliver,2000). In a succint formulation, it states that molecular entities that interact are prone to share biological properties. An exemplary instance can be found in (Lavi et al.,2012): genes that appear co-expressed tend to be closer in a network of interactions. The straightforward approach for GBA is neighbour voting, where the label of a given node is predicted by letting its neighbours vote with their own labels (Ballouz et al.,2016). Neighbour voting has been improved into the so-called label propagation and network diffusion approaches, which generally allow further propagation into higher-order neighbourhoods (Cowen et al.,2017). Resorting to propagation beyond neighbour nodes is endorsed by the network parsimony principle, that supports that the underlying perturbations propagate through the shortest paths within the complex molecular networks (Massucci et al.,2016). However, purely shortest paths-based approaches suffer from the “small world” property of biological networks: most nodes can be reached from every other node in a small number of steps due to the presence of hubs, i.e. highly connected nodes (Cowen et al.,2017). The use of label propagation techniques has been extensively
26 state of the art reviewed for finding social communities (Zoidi et al.,2015) and genetic associations (Cowen et al.,2017). (A) (B) Figure 10: Network diffusion example. (A) Three seed nodes are labelled as positive before applying the network diffusion. (B) After the propagation, all the nodes earn a diffusion score. Its magnitude is represented by a heat colour scale: red scores are the highest whereas white are the lowest. The closer to the seed nodes, the higher the score. The following sections provide several angles to tackle with the network diffusion and the random walks paradigms. A recurrent topic is how different approaches and physical models lead to equivalent formulations. Physical models An intuitive way to define diffusion-based approaches is through physical models, illustrated in figure 10. Equation 6contains the equation of a fluid propagation model: ∂fs(t) ∂t = −Lγfs(t) + bsu(t)(6) fs(t)is the column vector containing the amount of fluid in every node at time t.Lγ=L+γI, being Lthe graph Laplacian from equation 2,Ithe identity matrix and γ∈Ra parameter to control the rate of fluid leaking in every node. bsis the vector indicating the rate at which the fluid is pumped on the source nodes, and u(t)is the unit step function. The fluid densities fsin the stationary state (t→∞) are given by equation 7. fs=L−1 γbs(7) HotNet (Vandin et al.,2011) uses the model in equation 6, placing sources on one gene at a time and regarding fsin equation 7as the influence of that gene to all the genes in the network. Influence values are used to build an influence graph, in order to find subnetworks with mutations in a statistically significant number of patients. In HotNet, the source node diffuses 1 positive unit and the rest of nodes diffuse 0units of flow; therefore, bsis a binary vector. HotNet2is a second iteration with the same purpose, based on insulated diffusion processes, which can also be formulated in terms of random walks with restarts (Mark DM Leiserson et al.,2015) as in equation 17.
2.2 network propagation algorithms 27 Likewise, TieDIE (Paull et al.,2013) defines two diffusion processes to find ‘linker genes’ between two sets of genes: the source set (mutated genes) and the target set (transcription factors). The authors try three scoring schemes: the HotNet formulation (equation 7), PageRank (equation 17) and the pathway impact score from SPIA (Adi Laurentiu Tarca et al.,2008). eQED is a modified version of an electrical model to accomodate for directionality (Suthram et al.,2008), based on the random walk matrix in equation 4(Cowen et al.,2017). eQED prioritises causal genes, among those close to a genetic marker, for downstream gene expression changes. GeneMANIA (Mostafavi et al.,2008) predicts gene functions through diffusion processes on multiple networks that represent complementary data sources. The networks are combined using weights that maximise a kind of kernel-target alignment (Cristianini et al.,2002). The diffusion process is solved as in equation 8, equivalent to equation 7with y=bs,f=fs,γ=1. f= (I+L)−1y(8) The input yis defined differently; nodes are divided into positives (genes with the property of interest), negatives (genes with other properties) and unlabelled (nodes to be prioritised). The positives diffuse 1positive unit, like HotNet. However, negative nodes diffuse −1units, whereas unlabelled nodes diffuse a bias term k=n+−n− nthat accounts for the balance between the number of positive n+and negative n−instances over the number of nodes n. On the other hand, the diffusion problem in equation 6can be posed as the convex optimisation instance in equation 9with y=bsand f=fs(Tsuda et al.,2005); the parameter cis a tradeoff between loss and smoothness. min f(f−y)T(f−y) + cfTLf (9) Its solution is again similar to the classical diffusion problem (equation 7): f= (I+cL)−1y(10) The authors in (Tsuda et al.,2005) reformulated it to: min f,γ(f−y)T(f−y) + cγ fTLf ⩽γ(11) The convex problem was further generalised to accomodate knetworks (equation 12), with an application to protein function prediction. Analogously to GeneMANIA, the negatives are forced to diffuse −1units in this approach, whereas unlabelled nodes diffuse 0units. min f,γ(f−y)T(f−y) + cγ fTLkf⩽γ;k=1,...,m(12) Similar formulations have been used to derive supervised and unsupervised classification algorithms that favour smooth predictors on the network. This idea has found its application in microarray data, in order to leverage a priori network data (Rapaport et al.,2007).
34 state of the art the distinct roles and properties of the metabolites within biological processes. MetPA (Xia and Wishart,2010), part of MetaboAnalyst, computes centrality measures on a metabolite level and reports a topological impact of the metabolites in the user-provided list. mummichog predicts functional metabolic activity through module analysis within metabolic networks followed by pathway analysis (S. Li et al.,2013). Other tools revolve around curating and visualising network data. For instance, Metscape 2(Karnovsky et al.,2011) buils and displays the so-called CREG networks, connecting compounds, reactions, enzymes and genes. TP methods are bounded by the network data limitations, such as biases and incompleteness (Bayerlová et al.,2015). Current challenges include accounting for pathway cross-talk (Donato et al.,2013), the consideration of organism-specific data and the interpretability of the results (Booth et al.,2013). In a real scenario, however, TP analysis does not necessarily outpeform simpler tests in gene expression data (Bayerlová et al.,2015). The differential traits of certain enrichment methods are not an automatic guarantee of their optimality, which should be assessed by their ability to recover truly affected pathways (Mitrea et al.,2013). This is further hindered by the lack of standard datasets for the evaluation of such methods (Mitrea et al.,2013) and the statistical challenges of metabolomics, such as the unknown size of the metabolome and the sparsity of annotations compared to other omics sciences (Chagoyen and Pazos,2012). In addition, current evidence supports the inexistence of a universally optimal pathway enrichment technique (Adi L Tarca et al.,2013), adding an extra layer of complexity when defining the direction of future efforts. 2.3.2 Disease gene identification The identification of novel therapeutic targets is an area of active research. A wide spectrum of network-based approaches have been developed for this purpose or similar problems. Efforts include neighbour voting (Ballouz et al.,2016), semi-supervised learning (Valentini, Armano, et al.,2016), propagation and random walks (Vanunu et al.,2010), artificial neural networks (Muslu et al.,2019), supervised learning on diffusion-based features (H. Cho et al.,2016) or on diffusion-based distances (M. Cao et al.,2013). Methods make use of a variety of networks, ranging from a single interactome (Vanunu et al.,2010) to supervised weighted combinations of networks from various data sources (Mostafavi et al.,2008;Tsuda et al.,2005;Valentini, Paccanaro, et al.,2014). The assessment of the real benefit of integrating multiple sources is endorsed by some authors (Valentini, Paccanaro, et al., 2014), whereas others have found only a marginal improvement, if any, over a plain averaging of the networks (Mordelet and Vert,2011;Tsuda et al., 2005). The impact of the network coverage has also been examined. In line with the robustness to noise in diffusion-based methods, which are able to downweigh spurious predictions (Cowen et al.,2017), it has been reported that the usage of a larger network outweighs the higher proportion of noisy and low-confidence edges (J. K. Huang et al.,2018).
2.4 open issues 35 Drug target data suffers from the true negative issue: reliable data about truly unsuccessful targets is extremely rare, so the negative class becomes fuzzier (Ferrero et al.,2017). Another limitation derives from the fact that targeting is usually known at the protein complex level (or even for a protein family) instead of the protein sub-unit level (Bento et al.,2014). The label is therefore shared by all the genes that code the proteins within that protein complex, implying that the data is structured. Usual cross validation techniques yield misleading performance estimates on structured data if uncorrected (D. R. Roberts et al.,2017). Similar problems with cross validation on structured data have been identified and corrected in other fields, namely ligand-target binding models (Lopez-del Rio et al.,2018) and ecology (D. R. Roberts et al.,2017). 2.4 open issues Several limitations and challenges in the application of network-based and pathway analysis approaches in computational biology were mentioned throughout this chapter. The following sections highlight specific issues of special interest within the scope of this thesis. 2.4.1 Heterogeneity and biases in network propagation One of the first steps when applying network propagation to a new problem is the decision of how to propagate the data on the network. The choice of the graph kernel, the treatment of positive, negative and unlabelled nodes and the need of a statistical normalisation are open questions. An implementation that allows a systematic benchmark of these options is still missing. Besides, both the networks –for instance, protein-protein interaction networks (Edwards et al.,2002)– and the data that is propagated on them suffer from incompleteness and spurious associations. Despite their robutsness, diffusion-based approaches are still affected to a certain extent that has not been thoroughly characterised. On the other hand, network topology has been proven to affect diffusion scores (Erten, Bebek, et al.,2011). A plethora of network data resources is publicly available, with differences in data sources, coverage, topology and confidence (J. K. Huang et al.,2018). The network choice greatly affects downstream analysis (J. K. Huang et al.,2018), including diffusion-based approaches, which were used in their original study. In addition, the coexistence of well-studied and barely known genes or proteins, respectively turning into high and low-degree nodes, poses a challenge. In (Erten, Bebek, et al.,2011), propagation methods are shown to better predict highly connected proteins, but known disease genes are also biased towards highly connected genes. This circularity hampers the discovery of novel disease genes among the less studied ones. Few publications have explored how sensitive diffusion scores are. The authors in (Bersanelli et al.,2016) quantify an empiric p-value by permuting the input labels, in order to account for nodes with systematically low
36 state of the art or high scores. A similar concept has been applied to gene expression tstatistics (Cun and Fröhlich,2013) and to disease gene prioritisation (Erten, Bebek, et al.,2011), whereas (Biran et al.,2019) rewire network connections instead of permuting. There is enough evidence supporting the existence of topology-related biases and a noticeable impact upon their removal. The node degree seems to have a clear biasing role, but is unable to explain the whole casuistic of score behaviour (Hill et al.,2019). This further encourages its characterisation and quantification, including factors like the non-observability of certain nodes (i.e. due to experimental limitations). 2.4.2 Results interpretability in pathway analysis Even though pathway analysis was conceived to improve the biological interpretation of experimental data, understanding a list of affected pathways still remains as an outstanding challenge due to pathway overlap and cross-talk effects (Donato et al.,2013). Active research lines include the addition of richer organism-specific contextual representations (Booth et al., 2013), the modelling of pathway cross-talk (Donato et al.,2013) and the creation of aggregated pathway databases that better reflect current knowledge (Domingo-Fernandez et al.,2019). 2.4.3 Performance overestimation in target gene prediction Drug target data is often known for protein complexes instead of their individual protein sub-units (Bento et al.,2014). The presence of data structure can artificially inflate the peformance estimates from classical crossvalidation (Lopez-del Rio et al.,2018;D. R. Roberts et al.,2017). In addition, the performance metrics require a careful consideration. Classical metrics like the Area Under the Receiver Operating Characteristic can be misleading in early-retrieval (Saito and Rehmsmeier,2015), like the practical scenario in which only few targets can be tested. A comprehensive study controlling both factors is needed to obtain a realistic snapshot of the expected benefit, if any, of applying network propagation methods for drug discovery. 2.4.4 Free and open source software Competitive algorithms can be found in the literature for most of the areas in computational biology. However, the availability of their software, source code and the interaction with the user is variable across their spectrum. It is essential to provide the source code and the data that generates the conclusions of any manuscript to achieve reproducible science (Peng,2011). The lack of the data or source code that support the findings hinders their replication and a wider method adoption. Certain algorithms are available through a web server, preventing the user from customising its settings, modifying the algoritm or its application to other salient problems in computational biology. For instance, (Kamburov et al.,2011) offers a user-friendly web server for pathway enrichment that, on the other hand, does not contemplate changing the pathway libraries. En-
2.4 open issues 37 richNet (Glaab et al.,2012) is available on a web server that also offers an API using RESTful calls, enabling the user some programmatic options and network customisation. The web server MetaboAnalyst (Chong et al.,2018) has deployed a companion R package to provide batch analysis, reproducibility and transparency. Another common practice is to provide the raw data and the scripts that were used in the publication, like in MashUp (H. Cho et al.,2016). This policy is commendable, albeit still limited by the lack of maintenance over time and the non-standard distribution of the software, sometimes depending on a private or institutional server. Public repositories, either general purpose like CRAN11 (R Core Team, 2018) or specialised like Bioconductor12 (Huber et al.,2015), are a robust solution to endorse good coding practices, sofware maintenance and support, reproducibility, data availability and standard distribution channels for the R computing language. This is the case of RANKS (Valentini, Armano, et al.,2016), available in CRAN, and EGAD (Ballouz et al.,2016), published in Bioconductor. There is an analogous initiative for the python programming language community, called Biopython (Cock et al.,2009). An effort is needed not only in software publication in public repositories, but also in proper maintenance and long term support. The netClass R package (Cun and Fröhlich,2014) serves as an example: it was published in CRAN in 2013, but archived by July 29th, 2017 due to uncorrected check problems. This becomes an obstacle for its adoption and for reproducible science, since the package might need bug fixes to work on latest R versions. 11 https://cran.r-project.org. Accessed on 31/12/2019. 12 https://www.bioconductor.org. Accessed on 31/12/2019.
38 state of the art references 1000 Genomes Project Consortium and others 2015 “A global reference for human genetic variation”, Nature,526,7571, p. 68. Aggio, Raphael BM, Katya Ruggiero, and Silas Granato Villas-Bôas 2010 “Pathway Activity Profiling (PAPi): from the metabolite profile to the metabolic pathway activity”, Bioinformatics,26,23, pp. 2969- 2976. Aguirre-Plans, Joaquim, Janet Piñero, Jörg Menche, Ferran Sanz, Laura Furlong, Harald Schmidt, Baldo Oliva, and Emre Guney 2018 “Proximal pathway enrichment analysis for targeting comorbid diseases via network endopharmacology”, Pharmaceuticals,11,3, p. 61. Aimo, Lucila, Robin Liechti, Nevila Hyka-Nouspikel, Anne Niknejad, Anne Gleizes, Lou Götz, Dmitry Kuznetsov, Fabrice PA David, F Gisou van der Goot, Howard Riezman, et al. 2015 “The SwissLipids knowledgebase for lipid biology”, Bioinformatics, 31,17, pp. 2860-2866. Alanis-Lobato, Gregorio, Miguel A Andrade-Navarro, and Martin H Schaefer 2016 “HIPPIE v2.0: enhancing meaningfulness and reliability of protein– protein interaction networks”, Nucleic acids research, gkw985. Amberger, Joanna S, Carol A Bocchini, Alan F Scott, and Ada Hamosh 2018 “OMIM.org: leveraging knowledge across phenotype–gene relationships”, Nucleic acids research,47, D1, pp. D1038-D1043. Bader, Gary D, Michael P Cary, and Chris Sander 2006 “Pathguide: a pathway resource list”, Nucleic acids research,34, suppl 1, pp. D504-D506. Ballouz, Sara, Melanie Weber, Paul Pavlidis, and Jesse Gillis 2016 “EGAD: ultra-fast functional analysis of gene networks”, Bioinformatics,33,4, pp. 612-614. Bang-Jensen, Jørgen and Gregory Z Gutin 2008 Digraphs: theory, algorithms and applications, Springer Science & Business Media. Bánky, Dániel, Gábor Iván, and Vince Grolmusz 2013 “Equal opportunity for low-degree network nodes: a PageRankbased method for protein target identification in metabolic graphs”, PLoS One,8,1, e54204. Bayerlová, Michaela, Klaus Jung, Frank Kramer, Florian Klemm, Annalen Bleckmann, and Tim Beißbarth 2015 “Comparative study on gene set and pathway topology-based enrichment methods”, BMC bioinformatics,16,1, p. 334.
References 39 Beißbarth, Tim and Terence P Speed 2004 “GOstat: find statistically overrepresented Gene Ontologies within a group of genes”, Bioinformatics,20,9, pp. 1464-1465. Bento, A Patrícia, Anna Gaulton, Anne Hersey, Louisa J Bellis, Jon Chambers, Mark Davies, Felix A Krüger, Yvonne Light, Lora Mak, Shaun McGlinchey, et al. 2014 “The ChEMBL bioactivity database: an update”, Nucleic acids research,42, D1, pp. D1083-D1090. Bersanelli, Matteo, Ettore Mosca, Daniel Remondini, Gastone Castellani, and Luciano Milanesi 2016 “Network diffusion-based analysis of high-throughput data for the detection of differentially enriched modules”, Scientific reports,6, p. 34841. Biran, Hadas, Martin Kupiec, and Roded Sharan 2019 “Comparative analysis of normalization methods for network propagation”, Frontiers in genetics,10, p. 4. Booth, Sean C, Aalim M Weljie, and Raymond J Turner 2013 “Computational tools for the secondary analysis of metabolomics experiments”, Computational and structural biotechnology journal,4,5, e201301003. Cao, Mengfei, Hao Zhang, Jisoo Park, Noah M Daniels, Mark E Crovella, Lenore J Cowen, and Benjamin Hescott 2013 “Going the distance for protein function prediction: a new distance metric for protein interaction networks”, PloS one,8,10, e76339. Carter, Hannah, Matan Hofree, and Trey Ideker 2013 “Genotype to phenotype via network analysis”, Current opinion in genetics & development,23,6, pp. 611-621. Caspi, Ron, Richard Billington, Ingrid M Keseler, Anamika Kothari, Markus Krummenacker, Peter E Midford, Wai Kit Ong, Suzanne Paley, Pallavi Subhraveti, and Peter D Karp 2019 “The MetaCyc database of metabolic pathways and enzymes-a 2019 update”, Nucleic acids research. Chagoyen, Monica and Florencio Pazos 2012 “Tools for the functional interpretation of metabolomic experiments”, Briefings in bioinformatics,14,6, pp. 737-744. Chatr-Aryamontri, Andrew, Rose Oughtred, Lorrie Boucher, Jennifer Rust, Christie Chang, Nadine K Kolas, Lara O’Donnell, Sara Oster, Chandra Theesfeld, Adnane Sellam, et al. 2017 “The BioGRID interaction database: 2017 update”, Nucleic acids research,45, D1, pp. D369-D379. Cho, Hyunghoon, Bonnie Berger, and Jian Peng 2016 “Compact integration of multi-network topology for functional analysis of genes”, Cell systems,3,6, pp. 540-548.
40 state of the art Chong, Jasmine, Othman Soufan, Carin Li, Iurie Caraus, Shuzhao Li, Guillaume Bourque, David S Wishart, and Jianguo Xia 2018 “MetaboAnalyst 4.0: towards more transparent and integrative metabolomics analysis”, Nucleic acids research. Clough, Emily and Tanya Barrett 2016 “The gene expression omnibus database”, in Statistical Genomics, Springer, pp. 93-110. Cock, Peter JA, Tiago Antao, Jeffrey T Chang, Brad A Chapman, Cymon J Cox, Andrew Dalke, Iddo Friedberg, Thomas Hamelryck, Frank Kauff, Bartek Wilczynski, et al. 2009 “Biopython: freely available Python tools for computational molecular biology and bioinformatics”, Bioinformatics,25,11, pp. 1422- 1423. Consortium, Gene Ontology 2016 “Expansion of the Gene Ontology knowledgebase and resources”, Nucleic acids research,45, D1, pp. D331-D338. Consortium, UniProt 2018 “UniProt: a worldwide hub of protein knowledge”, Nucleic acids research,47, D1, pp. D506-D515. Cowen, Lenore, Trey Ideker, Benjamin J Raphael, and Roded Sharan 2017 “Network propagation: a universal amplifier of genetic associations”, Nature Reviews Genetics,18,9, p. 551. Cristianini, Nello, John Shawe-Taylor, Andre Elisseeff, and Jaz S Kandola 2002 “On kernel-target alignment”, in Advances in neural information processing systems, pp. 367-373. Cun, Yupeng and Holger Fröhlich 2013 “Network and data integration for biomarker signature discovery via network smoothed t-statistics”, PloS one,8,9, e73074. 2014 “Netclass: an r-package for network based, integrative biomarker signature discovery”, Bioinformatics,30,9, pp. 1325-1326. Diestel, Reinhard 2000 Graph Theory, Second edition, Graduate Texts in Mathematics, Springer, vol. 173. Domingo-Fernández, Daniel, Charles Tapley Hoyt, Carlos Bobis-Álvarez, Josep Marín-Llaó, and Martin Hofmann-Apitius 2018 “ComPath: an ecosystem for exploring, analyzing, and curating mappings across pathway databases”, NPJ systems biology and applications,5,1, p. 3. Domingo-Fernandez, Daniel, Sarah Mubeen, Josep Marin-Llao, Charles Hoyt, and Martin Hofmann-Apitius 2019 “PathMe: Merging and exploring mechanistic pathway knowledge”, bioRxiv, p. 451625.
References 41 Donato, Michele, Zhonghui Xu, Alin Tomoiaga, James G Granneman, Robert G MacKenzie, Riyue Bao, Nandor Gabor Than, Peter H Westfall, Roberto Romero, and Sorin Draghici 2013 “Analysis and correction of crosstalk effects in pathway analysis”, Genome research. Edwards, Aled M, Bart Kus, Ronald Jansen, Dov Greenbaum, Jack Greenblatt, and Mark Gerstein 2002 “Bridging structural biology and genomics: assessing protein interaction data with known complexes”, TRENDS in Genetics,18,10, pp. 529-536. ENCODE Project Consortium and others 2012 “An integrated encyclopedia of DNA elements in the human genome”, Nature,489,7414, p. 57. Erten, Sinan, Gurkan Bebek, Rob M Ewing, and Mehmet Koyutürk 2011 “DADA: degree-aware algorithms for network-based disease gene prioritization”, BioData mining,4,1, p. 19. Erten, Sinan and Mehmet Koyutürk 2010 “Role of centrality in network-based prioritization of disease genes”, in European Conference on Evolutionary Computation, Machine Learning and Data Mining in Bioinformatics, Springer, pp. 13-25. Everitt, Brian S 1992 The analysis of contingency tables, Chapman and Hall/CRC, chap. 2x2 Contingency tables. Fabregat, Antonio, Steven Jupe, Lisa Matthews, Konstantinos Sidiropoulos, Marc Gillespie, Phani Garapati, Robin Haw, Bijay Jassal, Florian Korninger, Bruce May, et al. 2017 “The reactome pathway knowledgebase”, Nucleic acids research,46, D1, pp. D649-D655. Ferrero, Enrico, Ian Dunham, and Philippe Sanseau 2017 “In silico prediction of novel therapeutic targets using gene–disease association data”, Journal of translational medicine,15,1, p. 182. El-Gebali, Sara, Jaina Mistry, Alex Bateman, Sean R Eddy, Aurélien Luciani, Simon C Potter, Matloob Qureshi, Lorna J Richardson, Gustavo A Salazar, Alfredo Smart, et al. 2018 “The Pfam protein families database in 2019”, Nucleic acids research, 47, D1, pp. D427-D432. Glaab, Enrico, Anaïs Baudot, Natalio Krasnogor, Reinhard Schneider, and Alfonso Valencia 2012 “EnrichNet: network-based gene set enrichment analysis”, Bioinformatics,28,18, pp. i451-i457.
42 state of the art Greene, Casey S, Arjun Krishnan, Aaron K Wong, Emanuela Ricciotti, Rene A Zelaya, Daniel S Himmelstein, Ran Zhang, Boris M Hartmann, Elena Zaslavsky, Stuart C Sealfon, et al. 2015 “Understanding multicellular function and disease with human tissuespecific networks”, Nature genetics,47,6, p. 569. Han, Heonjong, Jae-Won Cho, Sangyoung Lee, Ayoung Yun, Hyojin Kim, Dasom Bae, Sunmo Yang, Chan Yeong Kim, Muyoung Lee, Eunbeen Kim, et al. 2017 “TRRUST v2: an expanded reference database of human and mouse transcriptional regulatory interactions”, Nucleic acids research,46, D1, pp. D380-D386. Hartler, Jürgen 2015 “LIPID MAPS: Tools and Databases”, in Encyclopedia of Lipidomics, ed. by Markus R. Wenk, Springer Netherlands, Dordrecht, pp. 1-4, isbn:978-94-007-7864-1. Herwig, Ralf, Christopher Hardt, Matthias Lienhard, and Atanas Kamburov 2016 “Analyzing and interpreting genome data at the network level with ConsensusPathDB”, Nature protocols,11,10, p. 1889. Hill, Abby, Scott Gleim, Florian Kiefer, Frederic Sigoillot, Joseph Loureiro, Jeremy Jenkins, and Melody K Morris 2019 “Benchmarking network algorithms for contextualizing genes of interest”, PLOS Computational Biology,15,12, e1007403. Hoadley, Katherine A, Christina Yau, Denise M Wolf, Andrew D Cherniack, David Tamborero, Sam Ng, Max DM Leiserson, Beifang Niu, Michael D McLellan, Vladislav Uzunangelov, et al. 2014 “Multiplatform analysis of 12 cancer types reveals molecular classification within and across tissues of origin”, Cell,158,4, pp. 929- 944. Huang, Justin K, Daniel E Carlin, Michael Ku Yu, Wei Zhang, Jason F Kreisberg, Pablo Tamayo, and Trey Ideker 2018 “Systematic Evaluation of Molecular Networks for Discovery of Disease Genes”, Cell systems,6,4, pp. 484-495. Huang, Zhou, Jiangcheng Shi, Yuanxu Gao, Chunmei Cui, Shan Zhang, Jianwei Li, Yuan Zhou, and Qinghua Cui 2018 “HMDD v3.0: a database for experimentally supported human microRNA– disease associations”, Nucleic acids research,47, D1, pp. D1013-D1017. Huber, Wolfgang, Vincent J Carey, Robert Gentleman, Simon Anders, Marc Carlson, Benilton S Carvalho, Hector Corrada Bravo, Sean Davis, Laurent Gatto, Thomas Girke, et al. 2015 “Orchestrating high-throughput genomic analysis with Bioconductor”, Nature methods,12,2, p. 115.
References 43 Hwang, Sohyun, Chan Yeong Kim, Sunmo Yang, Eiru Kim, Traver Hart, Edward M Marcotte, and Insuk Lee 2018 “HumanNet v2: human gene networks for disease research”, Nucleic acids research,47, D1, pp. D573-D580. Jeske, Lisa, Sandra Placzek, Ida Schomburg, Antje Chang, and Dietmar Schomburg 2018 “BRENDA in 2019: a European ELIXIR core data resource”, Nucleic acids research,47, D1, pp. D542-D549. Jewison, Timothy, Yilu Su, Fatemeh Miri Disfany, Yongjie Liang, Craig Knox, Adam Maciejewski, Jenna Poelzer, Jessica Huynh, You Zhou, David Arndt, et al. 2013 “SMPDB 2.0: big improvements to the Small Molecule Pathway Database”, Nucleic acids research,42, D1, pp. D478-D484. Kale, Namrata S, Kenneth Haug, Pablo Conesa, Kalaivani Jayseelan, Pablo Moreno, Philippe Rocca-Serra, Venkata Chandrasekhar Nainala, Rachel A Spicer, Mark Williams, Xuefei Li, et al. 2016 “MetaboLights: An Open-Access Database Repository for Metabolomics Data”, Current protocols in bioinformatics,53,1, pp. 14-13. Kamburov, Atanas, Rachel Cavill, Timothy MD Ebbels, Ralf Herwig, and Hector C Keun 2011 “Integrated pathway-level analysis of transcriptomics and metabolomics data with IMPaLA”, Bioinformatics,27,20, pp. 2917-2918. Kanehisa, Minoru, Miho Furumichi, Mao Tanabe, Yoko Sato, and Kanae Morishima 2016 “KEGG: new perspectives on genomes, pathways, diseases and drugs”, Nucleic acids research,45, D1, pp. D353-D361. Karnovsky, Alla, Terry Weymouth, Tim Hull, V Glenn Tarcea, Giovanni Scardoni, Carlo Laudanna, Maureen A Sartor, Kathleen A Stringer, HV Jagadish, Charles Burant, et al. 2011 “Metscape 2bioinformatics tool for the analysis and visualization of metabolomics and gene expression data”, Bioinformatics,28,3, pp. 373-380. Khatri, Purvesh, Marina Sirota, and Atul J Butte 2012 “Ten years of pathway analysis: current approaches and outstanding challenges”, PLoS computational biology,8,2, e1002375. Köhler, Sebastian, Leigh Carmody, Nicole Vasilevsky, Julius O B Jacobsen, Daniel Danis, Jean-Philippe Gourdine, Michael Gargano, Nomi L Harris, Nicolas Matentzoglu, Julie A McMurry, et al. 2018 “Expansion of the Human Phenotype Ontology (HPO) knowledge base and resources”, Nucleic acids research,47, D1, pp. D1018-D1027.
50 state of the art Wishart, David S, Yannick Djoumbou Feunang, An C Guo, Elvis J Lo, Ana Marcu, Jason R Grant, Tanvir Sajed, Daniel Johnson, Carin Li, Zinat Sayeeda, et al. 2017 “DrugBank 5.0: a major update to the DrugBank database for 2018”, Nucleic acids research,46, D1, pp. D1074-D1082. Wishart, David S, Yannick Djoumbou Feunang, Ana Marcu, An Chi Guo, Kevin Liang, Rosa Vázquez-Fresno, Tanvir Sajed, Daniel Johnson, Carin Li, Naama Karu, et al. 2017 “HMDB 4.0: the human metabolome database for 2018”, Nucleic acids research,46, D1, pp. D608-D617. Wishart, David S, Carin Li, Ana Marcu, Hasan Badran, Allison Pon, Zachary Budinski, Jonas Patron, Debra Lipton, Xuan Cao, Eponine Oler, et al. 2019 “PathBank: a comprehensive pathway database for model organisms”, Nucleic acids research. Xia, Jianguo and David S Wishart 2010 “MetPA: a web-based metabolomics tool for pathway analysis and visualization”, Bioinformatics,26,18, pp. 2342-2344. Zerbino, Daniel R, Premanand Achuthan, Wasiu Akanni, M Ridwan Amode, Daniel Barrell, Jyothish Bhai, Konstantinos Billis, Carla Cummins, Astrid Gall, Carlos García Girón, et al. 2017 “Ensembl 2018”, Nucleic acids research,46, D1, pp. D754-D761. Zhang, Wei, Jeremy Chien, Jeongsik Yong, and Rui Kuang 2017 “Network-based machine learning and graph theory algorithms for precision oncology”, NPJ precision oncology,1,1, pp. 1-15. Zoidi, Olga, Eftychia Fotiadou, Nikos Nikolaidis, and Ioannis Pitas 2015 “Graph-based label propagation in digital media: A review”, ACM Computing Surveys (CSUR),47,3, p. 48.
3GOALS 3.1 main objective Diffusion scores are used in every discipline of computational biology that involves biological networks. On the other hand, concerns have arisen about the existence of a bias within the scores, potentially reaching numerous areas of active research. The main objective of this thesis is to develop, characterise and implement a statistical normalisation of diffusion scores. The potential benefits (or the absence thereof) will be assessed for salient problems in computational biology: pathway enrichment for metabolomics data and novel gene target discovery. 3.2 detailed objectives The main objective of this thesis can be achieved through three conceptual steps. First, a generic formulation of the normalisation used to address the bias. Then, its application to two computational biology domains: metabolomics data enrichment and prediction of sensible disease gene targets. 3.2.1 Conception of the statistical normalisation •Characterise and understand the bias in diffusion scores. •Define statistical models to normalise diffusion scores, focusing on providing a deterministic formulation. •Give a general guideline about when and how should diffusion scores be normalised. 3.2.2 Application to metabolomics data enrichment •Build a contextual representation linking metabolites to pathways as a knowledge graph. •Define a diffusion-based enrichment method, examine the need of a statistical normalisation. •Validate the method on in-house and public datasets. 3.2.3 Application to gene target discovery •Define a validation framework suitable for structured data and a performance metric oriented to drug development. 51
52 goals •Benchmark diffusion-based methods on a protein interaction network, including unnormalised and normalised scores. •Quantify the impact of the network choice and the disease under study. 3.3 expected contributions The main contribution will revolve around quantifying the presence of bias within the diffusion algorithms and providing ways to address it. On the other hand, the knowledge graph for metabolomics data enrichment has its interest per se, as it delivers a new paradigm for data interpretation. Every detailed objective is expected to lead to one or more publications in indexed scientific journals. To encourage open and reproducible science, all the algorithms and models will be released as free, open source tools. They will be encapsulated in R packages with extensive documentation and published in the Bioconductor repository.
4STATISTICAL PROPERTIES the effect of statistical normalisation on diffusion scores in computational biology Network diffusion and label propagation are fundamental tools in computational biology, with applications like gene-disease association, protein function prediction and module discovery. More recently, several publications have introduced a permutation analysis after the propagation process, due to concerns that network topology can bias diffusion scores. This opens the question of the statistical properties and the presence of bias of such diffusion processes in each of its applications. In this work, we characterised some common null models behind the permutation analysis and the statistical properties of the diffusion scores. We benchmarked seven diffusion scores on three case studies: synthetic signals on a yeast interactome, simulated differential gene expression on a protein-protein interaction network and prospective gene set prediction on another interaction network. For clarity, all the datasets were based on binary labels, but we also present theoretical results for quantitative labels. Diffusion scores starting from binary labels were affected by the label codification, and exhibited a problem-dependent topological bias that could be removed by the statistical normalisation. Parametric and non-parametric normalisation addressed both points by being codification-independent and by equalising the bias. We identified and quantified two sources of bias -mean value and variancethat yielded performance differences when normalising the scores. We provided closed formulae for both and showed how the null covariance is related to the spectral properties of the graph. Despite none of the proposed scores systematically outperformed the others, normalisation was preferred when the sought positive labels were not aligned with the bias. We conclude that the decision on bias removal should be problem and data-driven, i.e. based on a quantitative analysis of the bias and its relation to the positive entities. The code is publicly available at https://github.com/b2slab/diffuBench 4.1 introduction The guilt by association principle states that two proteins that interact with one another are prone to participate in the same, or related, cellular functions (Oliver,2000). This cornerstone fact has motivated the exploration This chapter is a reproduction of the following preprint, with minor section title changes: Picart-Armada, Sergio, Wesley K. Thompson, Alfonso Buil, and Alexandre Perera-Lluna. “The effect of statistical normalisation on network propagation scores”. BioRxiv (2020). 53
54 statistical properties of network algorithms on interaction networks for protein function prediction (Sharan et al.,2007). Network analysis has further proven its usefulness in other computational biology problems, such as prioritising candidate disease genes (Barabási et al.,2011), finding modular structures (Mitra et al., 2013) and modelling organisms (Aderem,2005). Network propagation is a fundamental formalism to leverage network data in computational biology. Its theoretical basis revolves around graph spectral theory, graph kernels and random walks (Smola and Kondor,2003). The central concept is that nodes carry abstract labels that, following the guilt by association principle, are propagated to the neighbouring nodes (Zoidi et al.,2015). Unlabelled nodes can therefore be inferred a label based on the available data of their neighbours. Label propagation can be defined in several ways, such as the heat diffusion, the electrical model or random walks with restarts (RWR), some of which lead to equivalent formulations (Cowen et al.,2017). One of the most common diffusion formulations relies on the regularised Laplacian graph kernel (Smola and Kondor,2003) - examples are provided throughout this paragraph. HotNet (Vandin et al.,2010) is a tool for finding modules with a statistically high number of mutated genes in cancer, after propagating the labels of mutated genes. The authors in (Bersanelli et al.,2016) have found relevant modules from gene expression and mutation data, based on a diffusion process followed by an automatic subgraph mining. GeneMANIA (Mostafavi et al.,2008) is a web server that predicts gene function by optimising a combination of knowledge networks and running a diffusion process on the resulting network. TieDIE (Paull et al.,2013) defines two diffusion processes in order to connect two sets of genes, applied to link perturbation in the genome with changes in the transcriptome. More generally, the predictive power of label propagation using graph kernels has been benchmarked in gene-disease association (Guala and Sonnhammer, 2017;Lee et al.,2011;Valentini et al.,2014). Some studies have pointed out biases in diffusion scores and explored the effect of their removal. The authors of DADA (Erten et al.,2011) have found that prioritisation using RWR favours highly connected genes and suggest several normalisation strategies. One of them computes a z-score that adjusts for the mean value and standard deviation estimated from propagation scores from random degree-preserving inputs. Another possibility is to normalise diffusion scores into empirical p-values, as used in the diffusion of t-statistics derived from gene expression (Cun and Fröhlich,2013). The aim was to quantify robust biomarkers, whose diffusion score is unlikely to arise from a permuted input. In the discovery of enriched modules (Bersanelli et al.,2016), the effect of the topology has been mitigated by combining diffusion scores with their empirical p-values. Similarly, exact z-scores and empirical p-values have been used for pathway analysis of metabolomics data (Sergio Picart-Armada, Fernández-Albert, et al.,2017). A recent study (Biran et al.,2019) has normalised RWR into an empirical p-value, obtained from edge rewiring. Specifically, random degree-preserving networks have been built to re-run the propagation and draw values from the null distributions of scores. Another recent manuscript (Hill et al.,2019) highlights biases in certain network propagation algorithms, related to the node degree.
4.2 approach 55 Overall, a variety of measures to address the bias have emerged, but a systematic quantification and evaluation of the biases is missing. The normalisation can potentially backfire, for instance by missing highly connected nodes that are associated with the property under study (Erten et al.,2011). The goal of this manuscript is to provide a quantitative way to assess the presence of the bias and its alignment with the node labels, in order to understand the impact and adequateness of the normalisation. 4.2 approach Here, we address basic statistical properties of the normalisation of singlenetwork diffusion scores to remove topology-related biases. We define and quantify two sources of bias. Both are derived from a statistical standpoint, based on the exact means and variances of the null distributions of the diffusion scores under input permutation. Differences in mean values between nodes should be the first indicator of systematic advantages: nodes with highest means will often be prioritised over those with lowest means. In their absence, differences in variances should be examined instead, as nodes with highest spread can be more likely to reach extreme scores. We compare classical and normalised propagation, as implemented in diffuStats (Sergio Picart-Armada, Thompson, et al.,2017), in data with and without bias. The main results are derived for the commonly used regularised Laplacian kernel, although most of them apply to other graph kernels and, to a lesser extent, to random walks with restarts. Special emphasis is placed on identifying scenarios under which normalisation is beneficial or detrimental and on understanding the underlying reasons why. 4.3 methods We include seven diffusion scores that are part of the diffuStats package (Sergio Picart-Armada, Thompson, et al.,2017): fraw,fml,fgm,fbers,fmc, fzand fberp. These scores are variations of the original diffusion model with a regularised unnormalised Laplacian kernel (Smola and Kondor,2003). Labelled nodes are referred to as positives if they have the property of interest, and negatives otherwise. 4.3.1 Unnormalised scores The starting point is the fraw score, which requires a graph kernel K (Smola and Kondor,2003) and input vector yraw and is computed as: fraw =Kyraw (27) This work focuses on the unnormalised, regularised Laplacian kernel for K, for being a widespread choice in the computational biology literature (electrical model, heat or fluid propagation). The values in yraw reflect
56 statistical properties the weights of each type of node: 1for positives and 0for negative and unlabelled entities. fml and fgm differ from fraw by setting a weight of −1on negative nodes. fgm also weighs unlabelled nodes with a bias term adapted from GeneMA- NIA (not to be confused with the diffusion bias). On the other hand, fbers measures the relative change between fraw and yraw, with a moderating parameter : fbers(i) = fraw(i) yraw(i) + (28) 4.3.2 Normalised scores Normalised scores attempt to equalise nodes that systematically show low or high scores, regardless of the input and due to the specific topology of the network. The lynchpin of normalisation is the null distribution of the diffusion scores under a random permutation πof the labelled nodes. The null scores arise from applying fraw to a randomised input Xy=π(yraw) and comparing, for the i-th node, fraw(i)to its null distribution Xf(i), where Xf=KXy. An empirical p-value can be computed throught Monte Carlo trials for the i-th node on Ntrials: p(i) = ri+1 N+1, (29) where riis the number of randomised trials having an equal or higher diffusion score in node i. In order to assign high scores to relevant nodes, the score is defined as fmc(i) = 1−p(i). We also include a parametric alternative to fmc by computing z-scores for each node i: fz(i) = fraw(i) − E(Xf(i)) pVar(Xf(i)) (30) The expected value and variance of the null distributions are analytically determined (see Supplement 1). Thus, fzhas a computational advantage over Monte Carlo trials. Finally, a hybrid combining an unnormalised and a normalised score is provided, inspired by how (Bersanelli et al.,2016) moderated the effect of hubs: fraw:fberp(i)=−log10(p(i))fraw(i). 4.3.3 Metrics and baselines Two baseline methods were used. First pagerank, regarded as an inputnaïve centrality measure (default damping factor of 0.85), to measure the predictive power of a basic network property. Second, a random predictor, to set an absolute baseline. Performances were quantified with two metrics: the area under the Receiver Operating Characteristic curve (AUROC) and the area under the Precision-Recall curve (AUPRC), as implemented in the precrec package (Saito and Rehmsmeier,2017). For clarity, the ranking (ordering) of the nodes for any given score and instance was normalised to lie
4.4 materials 57 in [0,1]by dividing it by the number of ranked nodes, so that top suggestions corresponded to ranks close to 0. 4.3.4 Bias quantification The reference expected value of the i-th node bK µ(i)(eq. 33) was defined as proportional to the expected value of its null distribution Xf(i)(eq. 31). Reference expected values that vary across nodes can indicate systematic differences in the diffusion scores of such nodes. In the absence of differences in the reference expected value, variancerelated bias was analysed instead. The reference variance of the i-th gene bK σ2(i)(eq. 34) was defined as, up to an additive constant, the base 10 logarithm of the variance of Xf(i), straightforward to obtain from the covariance matrix (eq. 32). The rationale is that the scores of nodes with varying dispersion measures should not be compared directly. 4.3.5 Performance explanatory models Explanatory models have found use in the formal description of differences in performance as a function of design factors (Lopez-del Rio et al., 2019;S Picart-Armada et al.,2019). Following (S Picart-Armada et al.,2019), the trends in AUROC and AUPRC were described through logistic-like quasibinomial models with a logit link function, as a generalisation of logistic models to prevent over and under-dispersion issues. Table 4presents the main model for each case study. The categorical regressors were: method,metric (AUROC or AUPRC), biased (refers to the signal, true or false), strat (labelled, unlabelled or overall), array (ALL or Lym), and the parameters k,rand p_max for the second case study. path_var_ref was quantitative, equal to the reference pathway variance bpK σ2(eq. 35). The responses were either AUROC, AUPRC, or both mixed, the latter denoted by Performance. 4.4 materials The evaluation of the diffusion scores was performed on three datasets of different nature, as described in Table 4: (1) synthetic signals on a yeast interactome, (2) pathway-based synthetic signals on a human network and (3) real signals on another human network. Table 4: Case studies for characterising biases and benchmarking diffusion scores. Interactions in explanatory models are denoted by a colon. Case Network Positive nodes Signal Bias type Purpose Explanatory model for hypothesis testing (1) Yeast Synthetic Synthetic, bias-based Mean value Proof of concept Performance ∼method +method :biased +metric (2) HPRD KEGG pathways Pathway sub-sampling Mean value Background influence in bias AUPRC ∼method +method :strat +array+k+r+pmax (3) BioGRID KEGG pathways Prospective pathway prediction Variance Bias in a common scenario AUROC ∼method +method :path_var_ref
58 statistical properties 4.4.1 Networks Yeast network A small yeast network was used to demonstrate the casuistic of diffusion scores properties. Medium and high confidence interactions from several sources were provided by the original study (Von Mering et al.,2002), as found in the igraphdata R package (Csardi,2015). It contains 2,617 proteins and 11,855 unweighted edges, but we worked only with its largest connected component (2,375 proteins, 11,693 edges). HPRD network The diffuse large B-cell lymphoma study, available in the R package DLBCL (Dittrich and Beisser,2010), contains a differential expression dataset accompanied by a human interactome network extracted from the Human Protein Reference Database, HPRD (Mishra et al.,2006). The original network encompasses 9,385 proteins with 36,504 interactions, whose largest connected component (8,989 nodes, 34,325 interactions) was extracted to compute the diffusion scores. We derived two gene backgrounds based on expression arrays. The first background (Lym) was taken from the expression data from 2,557 genes (2,482 in the network) in the lymphoma study (Rosenwald et al.,2002). The second background (ALL) was based on the acute lymphocytic leukemia array (Chiaretti et al.,2004), available in the ALL R package (Li,2009), encompassing 6,133 genes (5,921 in the network). BioGRID network The Biological General Repository for Interaction Datasets (BioGRID) (Chatraryamontri et al.,2017) is a public database with curated genetic and protein interaction from Homo sapiens and other organisms. BioGRID was retrieved in January 2017, but only keeping interactions dating from 2010 or older. The interactions were weighted according to (Cao et al.,2014), under the assumptions that more publications about an interaction boost its confidence and that low-throughput technologies are more reliable that high-throughput ones. The network encompassed 11,394 nodes and 67,573 edges and was connected. 4.4.2 Datasets Synthetic bias-based dataset 100 biased and 100 unbiased instances of positive, negative and unlabelled nodes were generated in dataset (1) from table 4, by sampling positive nodes with probabilities proportional to biased and unbiased scores. By construction, the frequencies of the positives drawn for biased signals were positively correlated with the reference expected value, whereas those of the unbiased signals were uncorrelated with it.
4.4 materials 59 Nodes were partitioned into three equally sized pools, from which positive nodes were drawn: (a) labelled nodes that were fed to the diffusion methods, (b) target nodes, the ones to be ranked and whose ground truth was known, and (c) filler nodes that were neither target nor labelled. For each instance, a fixed fraction of labelled nodes xewere uniformly sampled as positives, the rest of labelled nodes were deemed negatives and the target and filler nodes were left unlabelled. This input served two purposes: generate the ground truth in target nodes, and be the input for all the diffusion scores. To generate the ground truth in target nodes of biased signals, the raw diffusion scores were computed from the input above. A fixed fraction of target nodes xswas sampled with probabilities proportional to their raw scores, i.e. p(i)∝fraw(i), to become positives. The remaining target nodes would remain negatives, completing the ground truth. The regularised unnormalised Laplacian kernel is endorsed by physical models that ensure fraw(i)> 0 provided that inputs have one or more positives and the graph is connected. Analogously, unbiased signals were generated by sampling a fraction of target nodes xs, but with probabilities roughly proportional to the unbiased diffusion scores mc:p(i)∝fmc(i) + 1 N+1. By definition, the frequency of appearance of the target nodes was independent of the bias, and the small offset ensured p(i)> 0. In both cases, after sampling the ground truth, the same input was used again for all the diffusion scores, in order to rank the target nodes and compute the corresponding AUROC and AUPRC. Pathway sub-sampling dataset Synthetic gene expression statistics were generated, based on pathways in the Kyoto Encyclopedia of Genes and Genomes (KEGG) (Kanehisa et al.,2017), and on two array-based gene backgrounds described within the HPRD network. Genes outside the background were hidden (unlabelled), and genes inside were given p-values for differential expression. Each signal derived from krandom KEGG pathways. The pathways were assumed to be affected as a whole, but only a sampled portion of rgenes showed differential expression patterns. The p-values of the differential expressed genes were uniformly sampled from [0,pmax], whereas the rest of genes were uniform in [0,1], following a previously study (Rajagopalan and Agarwal,2004). For both expression arrays, genes with an FDR < 10% within their background were used as positives, the remaining background genes as negatives and the hidden nodes were deemed unlabelled. Notice that, by definition, this procedure generated false positives and false negatives among the input genes. The target genes were those belonging to the kaffected pathways, including those with no apparent differential expression and those among the unlabelled nodes. Methods were compared using the AUROC and AUPRC, computed separately on labelled, unlabelled genes, and overall, on a grid of parameters: k∈{1,3,5},r∈{0.3,0.5,0.7}and pmax ∈{10−2,10−3,10−4}. For each combination of parameters, N=50 instances were simulated.
66 statistical properties 4.6 conclusion In this study, we ratified that diffusion scores are biased due to the graph topology. We introduced two direct quantifications of the bias, in terms of the expected value and variance of the null distribution of the diffusion scores under input permutation. We analysed the benefits and pitfalls of using unbiased, statistically normalised scores and discussed several choices of the label weights when defining the diffusion process. We proved equivalences between scores under certain conditions, helping simplify the setup of the diffusion, and discovered that normalised alternatives are invariant under label weights changes. We found an explicit link between principal directions of the null covariance and the spectral features of the network. We applied the diffusion-based prioritisation on three scenarios: two with a mean value-related bias and one with a variance-related bias. Class imbalance and node topology had an impact in unnormalised scores, whereas normalised scores were more robust to both phenomena given their weightindependent definition. The parametric normalisation requires no permutations compared to Monte Carlo trials and performed equally or better, providing a convenient way to normalise. While mean value bias was straightforward to characterise, variance bias was less intuitive albeit of noticeable impact. In general terms, the statistical normalisation is advised if the positives are not aligned with the bias, and discouraged otherwise. The statistical background, i.e. which nodes are permuted, is a key piece that should be clearly stated in every application. Bias assessment should be carried through its direct quantification instead of indirect indicators, which can be misleading. We conclude that the statistical normalisation can be beneficial or detrimental, and the decision should follow from the dependence between the node bias and the hypothetical or desired properties of the new positives. Topology-related bias can manifest in different ways (mean value- or variancerelated bias) and each instance should be properly characterised. acknowledgements SP thanks Guillem Belda-Ferrín for reviewing the mathematical proofs. SP thanks Imanol Morata Martínez and Camellia Sarkar for fruitful discussions and useful suggestions. Conflict of Interest: none declared funding This work was supported by the Spanish Ministry of Economy and Competitiveness (MINECO) [TEC2014-60337-R and DPI2017-89827-R to A.P.] and the National Institutes of Health (NIH) [R01GM104400 to W.T.]. AP and SP thank CIBERDEM and CIBER-BBN for funding, both initiatives of the Spanish ISCIII. SP thanks the AGAUR FI-scholarship programme.
References 67 references Aderem, Alan 2005 “Systems biology: its practice and challenges”, Cell,121,4, pp. 511- 513. Barabási, Albert-László, Natali Gulbahce, and Joseph Loscalzo 2011 “Network medicine: a network-based approach to human disease”, Nature reviews. Genetics,12,1, p. 56. Bersanelli, Matteo, Ettore Mosca, Daniel Remondini, Gastone Castellani, and Luciano Milanesi 2016 “Network diffusion-based analysis of high-throughput data for the detection of differentially enriched modules.” Scientific Reports,6, August, p. 34841. Biran, Hadas, Martin Kupiec, and Roded Sharan 2019 “Comparative analysis of normalization methods for network propagation”, Frontiers in genetics,10, p. 4. Cao, Mengfei, Christopher M Pietras, Xian Feng, Kathryn J Doroschak, Thomas Schaffner, Jisoo Park, Hao Zhang, Lenore J Cowen, and Benjamin J Hescott 2014 “New directions for diffusion-based network prediction of protein function: incorporating pathways with confidence”, Bioinformatics, 30,12, pp. i219-i227. Chatr-aryamontri, Andrew, Rose Oughtred, Lorrie Boucher, Jennifer Rust, Christie Chang, Nadine K Kolas, Lara O’Donnell, Sara Oster, Chandra Theesfeld, Adnane Sellam, et al. 2017 “The BioGRID interaction database: 2017 update”, Nucleic acids research,45, D1, pp. D369-D379. Chiaretti, Sabina, Xiaochun Li, Robert Gentleman, Antonella Vitale, Marco Vignetti, Franco Mandelli, Jerome Ritz, and Robin Foa 2004 “Gene expression profile of adult T-cell acute lymphocytic leukemia identifies distinct subsets of patients with different response to therapy and survival”, Blood,103,7, pp. 2771-2778. Cowen, Lenore, Trey Ideker, Benjamin J Raphael, and Roded Sharan 2017 “Network propagation: a universal amplifier of genetic associations”, Nature Reviews Genetics,18,9, pp. 551-562. Csardi, Gabor 2015 igraphdata: A Collection of Network Data Sets for the ’igraph’ Package, R package version 1.0.1,https://CRAN.R-project.org/package= igraphdata. Cun, Yupeng and Holger Fröhlich 2013 “Network and Data Integration for Biomarker Signature Discovery via Network Smoothed T-Statistics”, PLoS One,8,9.
68 statistical properties Dittrich, Marcus and Daniela Beisser 2010 DLBCL: Diffuse large B-cell lymphoma expression data, R package version 1.16.0,http://bionet.bioapps.biozentrum.uni-wuerzburg. de/. Erten, Sinan, Gurkan Bebek, Rob M Ewing, and Mehmet Koyutürk 2011 “DADA: degree-aware algorithms for network-based disease gene prioritization”, BioData mining,4,1, p. 19. Guala, Dimitri and Erik LL Sonnhammer 2017 “A large-scale benchmark of gene prioritization methods”, Scientific Reports,7. Hill, Abby, Scott Gleim, Florian Kiefer, Frederic Sigoillot, Joseph Loureiro, Jeremy Jenkins, and Melody K Morris 2019 “Benchmarking network algorithms for contextualizing genes of interest.” PLoS Computational Biology,15,12. Kanehisa, Minoru, Miho Furumichi, Mao Tanabe, Yoko Sato, and Kanae Morishima 2017 “KEGG: new perspectives on genomes, pathways, diseases and drugs”, Nucleic acids research,45, D1, pp. D353-D361. Lee, Insuk, U Martin Blom, Peggy I Wang, Jung Eun Shim, and Edward M Marcotte 2011 “Prioritizing candidate disease genes by network-based boosting of genome-wide association data”, Genome Research,21,7, pp. 1109- 1121. Li, Xiaochun 2009 ALL: A data package, R package version 1.20.0. Lopez-del Rio, Angela, Alfons Nonell-Canals, David Vidal, and Alexandre Perera-Lluna 2019 “Evaluation of cross-validation strategies in sequence-based binding prediction using Deep Learning”, Journal of chemical information and modeling,59,4, pp. 1645-1657. Mishra, Gopa R, M Suresh, K Kumaran, N Kannabiran, Shubha Suresh, P Bala, K Shivakumar, N Anuradha, Raghunath Reddy, T Madhan Raghavan, et al. 2006 “Human protein reference database—2006 update”, Nucleic acids research,34, suppl 1, pp. D411-D414. Mitra, Koyel, Anne-Ruxandra Carvunis, Sanath Kumar Ramesh, and Trey Ideker 2013 “Integrative approaches for finding modular structure in biological networks”, Nat. Rev. Genet.,14,10, pp. 719-732.
References 69 Mostafavi, Sara, Debajyoti Ray, David Warde-Farley, Chris Grouios, and Quaid Morris 2008 “GeneMANIA: a real-time multiple association network integration algorithm for predicting gene function.” Genome Biology,9Suppl 1, S4. Oliver, Stephen 2000 “Guilt-by-association goes global”, Nature,403, February, pp. 601- 603. Paull, Evan O., Daniel E. Carlin, Mario Niepel, Peter K. Sorger, David Haussler, and Joshua M. Stuart 2013 “Discovering causal pathways linking genomic events to transcriptional states using Tied Diffusion Through Interacting Events (TieDIE)”, Bioinformatics,29,21, pp. 2757-2764. Picart-Armada, S, SJ Barrett, DR Willé, A Perera-Lluna, A Gutteridge, and BH Dessailly 2019 “Benchmarking network propagation methods for disease gene identification”, PLoS Comput Biol,15,9, e1007276. Picart-Armada, Sergio, Francesc Fernández-Albert, Maria Vinaixa, Miguel A Rodríguez, Suvi Aivio, Travis H Stracker, Oscar Yanes, and Alexandre Perera-Lluna 2017 “Null diffusion-based enrichment for metabolomics data”, PloS one, 12,12, e0189012. Picart-Armada, Sergio, Wesley K Thompson, Alfonso Buil, and Alexandre Perera-Lluna 2017 “diffuStats: an R package to compute diffusion-based scores on biological networks”, Bioinformatics,34,3, pp. 533-534. Rajagopalan, Dilip and Pankaj Agarwal 2004 “Inferring pathways from gene lists using a literature-derived network of biological relationships”, Bioinformatics,21,6, pp. 788-793. Rosenwald, Andreas, George Wright, Wing C Chan, Joseph M Connors, Elias Campo, Richard I Fisher, Randy D Gascoyne, H Konrad Muller- Hermelink, Erlend B Smeland, Jena M Giltnane, et al. 2002 “The use of molecular profiling to predict survival after chemotherapy for diffuse large-B-cell lymphoma”, New England Journal of Medicine, 346,25, pp. 1937-1947. Saito, Takaya and Marc Rehmsmeier 2017 “Precrec: fast and accurate precision–recall and ROC curve calculations in R”, Bioinformatics,33,1, pp. 145-147. Sharan, Roded, Igor Ulitsky, and Ron Shamir 2007 “Network-based prediction of protein function”, Molecular systems biology,3,1, p. 88.
70 statistical properties Smola, Alexander J and Risi Kondor 2003 “Kernels and regularization on graphs”, in Learning theory and kernel machines, Springer, pp. 144-158. Valentini, Giorgio, Alberto Paccanaro, Horacio Caniza, Alfonso E. Romero, and Matteo Re 2014 “An extensive analysis of disease-gene associations using network integration and fast kernel-based gene prioritization methods”, Artificial Intelligence in Medicine,61,2, pp. 63-78. Vandin, Fabio, Eli Upfal, and Benjamin J. Raphael 2010 “Algorithms for detecting significantly mutated pathways in cancer”, Lect. Notes Comput. Sci.,6044 LNBI, 3, pp. 506-521. Von Mering, Christian, Roland Krause, Berend Snel, Michael Cornell, Stephen G Oliver, Stanley Fields, and Peer Bork 2002 “Comparative assessment of large-scale data sets of protein–protein interactions”, Nature,417,6887, p. 399. Zoidi, Olga, Eftychia Fotiadou, Nikos Nikolaidis, and Ioannis Pitas 2015 “Graph-Based Label Propagation in Digital Media: A Review”, ACM Computing Surveys,47,3,48:1-48:35.
5THE R PACKAGE DIFFUSTATS diffustats: an r package to compute diffusion-based scores on biological networks Label propagation and diffusion over biological networks are a common mathematical formalism in computational biology for giving context to molecular entities and prioritising novel candidates in the area of study. There are several choices in conceiving the diffusion process (involving the graph kernel, the score definitions and the presence of a posterior statistical normalisation) which have an impact on the results. This manuscript describes diffuStats, an R package that provides a collection of graph kernels and diffusion scores, as well as a parallel permutation analysis for the normalised scores, that eases the computation of the scores and their benchmarking for an optimal choice. The R package diffuStats is publicly available in Bioconductor, https://bioconductor.org, under the GPL-3license. 5.1 introduction Network analysis can help finding therapeutic targets and understanding biology in networks obtained from protein-protein interactions, gene regulation and metabolic reactions. In this context, label propagation and diffusion algorithms (Zoidi et al.,2015) address a general problem of molecular entity ranking according to a seed node list. Examples include finding significantly mutated subnetworks in cancer (Vandin et al.,2010), predicting gene function (Mostafavi et al.,2008), prioritising genome-wide association hits (Lee et al.,2011) and classifying proteins (Tsuda et al.,2005). In general, the mentioned methods involve diffusion processes with adhoc parameter and network settings, making comparisons fundamentally difficult. The RANKS R package (Valentini et al.,2016) is an effort to collect a range of diffusion kernels and scores, but the effects of label codification and a recently proposed statistical normalisation (Bersanelli et al.,2016) have not been explored. To that end, we introduce the diffuStats R package gathering diffusion kernels, input codifications and statistical normalisations to benchmark single-network diffusion settings. This chapter is a postprint of the following journal article: Picart-Armada, Sergio, Wesley K. Thompson, Alfonso Buil, and Alexandre Perera-Lluna. “diffuStats: an R package to compute diffusion-based scores on biological networks”. Bioinformatics 34, no. 3(2018): 533-534. 71
72 the r package diffustats Table 5: Implemented diffusion scores. fraw,fml and fgm differ on the weights of the positive, negative and unlabelled nodes. fbersquantifies the change of the fraw scores relative to the input scores. fberp,fmc and fzare statistically normalised by permuting the labelled examples, but not the unlabelled*. fmc derives from an empirical p-value, whereas fberpcombines fmc and fraw.fzis a parametric alternative to fmc requiring no stochastic permutations. Quantitative inputs are allowed in all the scores except fml and fgm. Score y+y−yuNormalised Stochastic Quantitative Reference raw1 0 0 No No Yes (Vandin et al.,2010) ml 1-1 0 No No No (Tsuda et al.,2005) gm 1-1kNo No No (Mostafavi et al.,2008) bers1 0 0 No No Yes (Bersanelli et al.,2016) berp1 0 0* Yes Yes Yes (Bersanelli et al.,2016) mc 1 0 0* Yes Yes Yes (Bersanelli et al.,2016) z1 0 0* Yes No Yes (Harchaoui et al.,2013) 5.2 methods The diffuStats R package offers scoring schemes for diffusing a label vector on a network, determined by (a) the graph kernel, (b) the translation of labels into a numeric vector yto be smoothed, and (c) the statistical normalisation. In general, diffusion scores fare based on modifications of the quantity f=K·y, where yare the input labels, fthe diffusion scores and Kis a graph kernel. Regarding (a), most of the cited applications use the regularised Laplacian kernel, but our package also offers the diffusion kernel, the p-step random walk kernel, the inverse cosine kernel (Smola and Kondor,2003) and the commute time kernel (Yen et al.,2007). In practice, they differ on the reach and the behaviour of the spreading inside the network - further detail for its choice can be found in the documentation. The decision can be data-driven, based on prior studies on the subject or on desirable properties of the kernel. For (b) and (c), the implemented scores are variations of the propagation of a binary vector whose ones are the positive labels y+and whose zeroes are the negative y−and unlabelled yuentities (Table 5). The statistical normalisation (c) compares the diffusion scores with the distribution of scores that arise from a permuted input, in order to spot nodes whose score is systematically high or low regardless of the input. However, not all normalised scores require stochastic simulations. The diffuStats R package contains proper documentation and unit testing to facilitate its development. Its algorithms are written in R except the permutations, which use C++; further details on the implementation can be found in the supplementary materials. Manipulating networks with more than 10,000 nodes might require extra RAM memory and computational power due to the kernel matrix size.
5.3 results 73 0.5 0.6 0.7 0.8 0.9 raw ml gm ber_s ber_p mc z Method Area under the curve Benchmark of protein categories Figure 15: Performance comparison for diffusion scores in predicting 12 functions on half of the proteins using the area under the Receiver Operating Characteristic curve. 5.3 results The example data is a yeast interactome with 12 annotated protein functions (Von Mering et al.,2002). The functionalities of our package are demonstrated by (i) obtaining a prioritised list of annotations given a set of labels, and (ii) benchmarking all the available diffusion scores in a dataset containing validation data. For both analyses, half of the proteins in the interactome will be treated as unlabelled and will receive a score from the propagation of the other half using the default regularised Laplacian kernel. Regarding (i), the function “diffuse” allows to compute the desired diffusion scores with a starting set of scores (labels) and a network: scores_diff <- diffuse(graph = yeast, scores = scores, method = "raw") When assessing the performance of different diffusion scores in a given dataset (ii), the desired metrics involving the diffusion scores and the validation can be computed on a grid of parameters: performance <- perf(graph = yeast, scores = scores, validation = validation, grid_param = grid_param) The results are returned as a table that can be directly plotted (Fig. 15). In this case study, statistically normalised scores fmc,fzand fberpseem preferable than their unnormalised counterparts comparing the area under the ROC curves. For instance, fzoutperforms fraw,fml,fgm and fbers(FDR <25%, Wilcoxon test), thus highlighting the usefulness of a prior screening in score performance and its potential impact. In summary, the R package diffuStats gathers diffusion kernels and scores with statistical normalisation that are object of active research in bioinformatics, like functional prediction or module identification. In addition, it
74 the r package diffustats facilitates benchmarking diffusion scoring methods to find the optimal configuration for the application of interest. funding: This work was supported by the Spanish Ministry of Economy and Competitiveness (MINECO) [TEC2014-60337-R to A.P.] and the National Institutes of Health (NIH) [R01GM104400 to W.T.]. AP. and S.P. thank CIBERDEM and CIBER-BBN for funding, both initiatives of the Spanish ISCIII. SP. thanks the AGAUR FI-scholarship programme. Conflict of Interest: none declared
References 75 references Bersanelli, Matteo, Ettore Mosca, Daniel Remondini, Gastone Castellani, and Luciano Milanesi 2016 “Network diffusion-based analysis of high-throughput data for the detection of differentially enriched modules”, Sci. Rep.,6. Harchaoui, Zaid, Francis Bach, Olivier Cappe, and Eric Moulines 2013 “Kernel-based methods for hypothesis testing: A unified view”, IEEE Signal Process Mag,30,4, pp. 87-97. Lee, Insuk, U Martin Blom, Peggy I Wang, Jung Eun Shim, and Edward M Marcotte 2011 “Prioritizing candidate disease genes by network-based boosting of genome-wide association data”, Genome Res.,21,7, pp. 1109-1121. Mostafavi, Sara, Debajyoti Ray, David Warde-Farley, Chris Grouios, and Quaid Morris 2008 “GeneMANIA: a real-time multiple association network integration algorithm for predicting gene function”, Genome Biol.,9,1, S4. Smola, Aj Alexander J and Risi Kondor 2003 “Kernels and Regularization on Graphs”, Mach. Learn.,2777, pp. 1- 15. Tsuda, Koji, HyunJung J. Shin, and Bernhard Schölkopf 2005 “Fast protein classification with multiple networks”, Bioinformatics, 21, SUPPL. 2, pp. 59-65. Valentini, Giorgio, Giuliano Armano, Marco Frasca, Jianyi Lin, Marco Mesiti, and Matteo Re 2016 “RANKS: A flexible tool for node label ranking and classification in biological networks”, Bioinformatics,32,18, pp. 2872-2874. Vandin, Fabio, Eli Upfal, and Benjamin J. Raphael 2010 “Algorithms for detecting significantly mutated pathways in cancer”, Lect. Notes Comput. Sci.,6044 LNBI, 3, pp. 506-521. Von Mering, Christian, Roland Krause, Berend Snel, Michael Cornell, Stephen G Oliver, Stanley Fields, and Peer Bork 2002 “Comparative assessment of large-scale data sets of protein–protein interactions”, Nature,417,6887, pp. 399-403. Yen, Luh, Francois Fouss, Christine Decaestecker, Pascal Francq, and Marco Saerens 2007 “Graph nodes clustering based on the commute-time kernel”, Advances in Knowledge Discovery and Data Mining, pp. 1037-1045. Zoidi, Olga, Eftychia Fotiadou, Nikos Nikolaidis, and Ioannis Pitas 2015 “Graph-Based Label Propagation in Digital Media: A Review”, ACM Comput. Surv.,47,3,48:1-48:35.
82 metabolomics enrichment The arrangement of nodes for the PageRank calculation is identical to the one for diffusion (Fig 17b), being edges directed towards the upper levels. Random walks start only at the affected metabolites and explore all the reachable nodes. Further details are available in S3Appendix. 6.2.3 Null models The ranking of the network nodes is not achieved through raw scores, due to potential biases related to topological features. This is also the case in classical over-representation analysis, as it can be rephrased as a particular case of heat diffusion (Fig 18) where the observed statistic is the node temperature and its null distribution is the hypergeometric distribution. In view of this, our approach also includes a permutation analysis in the input, leading to a null distribution of scores for each node. Node scores are normalised using their null distributions and ranked, allowing a subgraph (Fig 16) to be extracted. Further details can be found in S4Appendix. Compounds Pathway Pathway ARest of nodes Figure 18: Toy example of an over-representation analysis of a hypothetical "pathway A" containing 3metabolites out of a total of 10. The list to be enriched contains 4metabolites, showing 2hits in the pathway. The corresponding (Fisher’s exact test) over-representation can be understood as a diffusion process on the depicted network followed by a null model. The temperature of pathway A is always coincident with the number of hits in the pathway, implying that its null distribution is the hypergeometric distribution, to which a one-tailed temperature comparison is made. The null model will be introduced in the heat diffusion scenario (the PageRank case is analogous). Let nin be the number of compounds in the input. Then, exactly nin different KEGG compounds are chosen at random following dependent Bernoulli distributions, so that Xi=1if iis chosen and Xi=0otherwise. Normalisation can be performed using (i) the theoretical mean and variance of the scores, which can be obtained from Eq. 36, using the fact that, for the null model, Gis a random vector Xwith known mean and covariance matrix: E(Tnull) = RHD ·E(X)(37) Σ(Tnull) = RHD ·Σ(X)·RT HD (38) The normalised score (z-score) of node iis defined in terms of the expected value µi=E(Tnull)iand standard deviation σi=pΣ(Tnull)i,i zi=Ti−µi σi (39)
6.2 materials and methods 83 Then, nodes with the top kscores are kept and reported. Alternatively, scores can be normalised through (ii) Monte Carlo simulations with nperm permutations, which provide an estimate of the probability pithat the null distribution attains a score greater than or equal to the observed one. Estimation of piinvolves the empirical cumulative distribution function with a small correction (North et al.,2002), ribeing the number of permutations in which the null score of node iis greater or equal than Ti: pi=ri+1 nperm +1(40) A consensus solution is derived from nvote independent sets of Monte Carlo trials, each trial reporting the top knodes. The consensus solution may therefore report a node count not exactly equal to k. 6.2.4 NMR validation The reported subgraphs contain entities other than pathways and compounds that can be useful for the researchers. Among these, the highlighted reactions have been partially validated by quantifying their distance to an independent second set of affected metabolites. In order to analyse the reactions in the scope of a metabolic network, distances are computed on the unweighted, maximal connected subgraph containing all the compounds and reactions from the KEGG graph, referred to as the reaction-compound graph. The validation metric is the resistance distance, previously used in the chemical literature (Bapat,2004). Under these settings, the reported reactions are compared to all the reactions that involve the input metabolites (their nearest neighbours) in terms of their resistance distance to the second set of metabolites. 6.2.5 Evaluation with synthetic signals In order to deploy an analysis of true and false positive pathway identifications, we opted to statistically characterize the pathway prioritisation induced by the diffusion scores. Artificial pathway signals have been generated to (a) find biases in the absence of a signal that might cause false positives, and to (b) quantify the ability to recover true positive pathways. The proposed methods are not directly compared to IMPaLA and MetaboAnalyst due to the lack of a batch analysis mode, but instead to their underlying distribution using Fisher’s exact test. Our Monte Carlo approaches have not been aggregated into consensus solutions. The performance metric is the pathway rank in the list ordered by a method, where 1 npis the best rank and 1is the worst one, npbeing the number of pathways in the KEGG graph. Ranks in Fisher’s exact test are computed using the raw p-values, so that top ranked pathways correspond to lowest p-values. To compute the p-values, a metabolite is considered to belong to a pathway if it can be reached via the pathway in our directed KEGG graph (Fig 17). In (a), noisy signals are generated and the ranks of all the pathways are calculated within signals. Then, the mean rank of a specific pathway iis
84 metabolomics enrichment computed across all the signals. This measure can reveal pathways that tend to have an extreme rank irrespective of the input. In (b), a target pathway generates the signal and its rank is used as the metric of interest. Methods able to recover the signal will show low ranks in general terms. 6.2.6 Description of the experimental data Our method has been tested using data from a case-control experiment aimed at determining the function of an uncharacterised mitochondrial protein by silencing the gene using short hairpin RNAs (shRNA). Metabolites abundances were determined from five replicates of cell cultures expressing either control or experimental shRNA. Metabolite measurements were performed by Metabolon platform (www. metabolon.com) using GC/MS (Thermo-Finnigan Trace DSQ single quadrupole) and LC/MS (Waters ACQUITY UPLC and a Thermo-Finnigan LTQFT). The proprietary Metabolon analysis reported 168 quantified metabolites annotated in the KEGG database. In addition, we have used NMR following the labelling of the same cells with [U-13C] glucose (DeBerardinis et al.,2007) to trace carbon atoms, in order to further validate the conclusions of our new method. The reported reactions are evaluated in terms of their resistance distance to the affected metabolites found by NMR. 6.2.7 Description of the synthetic data All the signals generate a list with fixed length nin =35 for each one of the nppathway nodes in the KEGG graph. Three sampling types have been defined – differences arise in the specification of how much more probable compounds in the target pathway are. The first signal is a uniform sampling of nin compounds that imitates noise: the probability of drawing a compound jwithin pathway i,pi,j, is ki=1times more likely to be drawn than compounds outside the pathway, and thus does not depend on the pathway. In the second signal, compounds belonging to pathway iare ki=10 times more likely to be drawn. Therefore, there are two different probability values: inside pathway and outside pathway. This sampling is affine to the assumptions in Fisher’s exact test from ORA. As for the third signal, pi,jis proportional to the quantity RHDij, which is greater in compounds close to the pathway. This takes into account the whole KEGG graph, thus being influenced by indirect connections and compound specificity.
6.3 results 85 6.3 results 6.3.1 Input for the algorithms After the curation step, our knowledge base graph contains 10,183 nodes and 31,539 edges. The nodes are stratified in 288 pathways, 178 modules, 1,149 enzymes, 4,699 reactions and 3,869 compounds. The degree distribution of its vertices follow a scale-free network model, where P(k)∼k−γ, with γ=2.084 ∈[2,3], see S1Appendix. On the other hand, MS led to 168 quantified metabolites from KEGG. Two identifiers that each appeared twice have been dropped, as well as a KEGG drug, excluded from the KEGG compound category. The remaining 163 metabolites have been tested between both conditions, leading to 38 significant metabolites (two-tailed non-parametric Wilcoxon, FDR < 0.05), of which 33 have been mapped to our KEGG graph. The 33 MS-derived compounds served as input for each of the proposed enrichment algorithms. Heat Diffusion (HD) and PageRank (PR) are followed by norm (z-score normalisation) or sim (Monte Carlo permutations). Normalised scores have been computed through the null models with nin = 33, followed with subgraph selection with a desired number of nodes k= 250. For simulated methods, a consensus subgraph using nvote =9runs of nperm =10,000 permutations each has been derived by majority vote on each node. Regardless of the specific details, high diffusion scores are an indicator of overall closeness to the MS-derived metabolites and potential relevance in the condition being studied. This intuition, known as guilt-by-association, can be phrased in the context of heat diffusion: high temperatures are found close to the heat sources. Therefore, warm nodes are candidates for further study as they are easily reached through database annotations from the input metabolites. 6.3.2 Null model impact The impact of using the null model in HD and an overview of the random temperatures behaviour is described in Fig 19. The null model is closely related to the graph structure and node topology, quantified through the vertex degree. In Fig 19a, the mean temperatures show different trends for the five levels in the graph; in particular, there is an increase in the mean pathway temperature as the pathway becomes larger. This implies that, regardless of the input, larger pathways will generally show warmer temperatures and the results will be biased towards them. Likewise, the standard deviations of the null temperatures show level-specific changes (Fig 19b), with the compounds being the most affected entities – the higher the degree of the compound, the lower its standard deviation. The usage of z-scores instead of raw temperatures has consequences in the highlighted nodes. Reporting the nodes with the top 250 raw temperatures does not reveal any pathway (Fig 19c), whereas five pathways lay among the top 250 z-scores (Fig 19d). Likewise, if only pathway nodes are considered, their ranking using raw temperatures is closely related to the ranking using
86 metabolomics enrichment (c) (d) Legend ●Pathway ●Module ●Enzyme ●Reaction ●Compound (input: ) (a) (b) Top 250 raw temperatures Top 250 normalised temperatures Figure 19: Expected value (a) and standard deviation (b) of the null temperatures, stratified by level – jitter applied for visual purposes and 0.95 confidence intervals computed by the default GAM models in ggplot2R library (Wickham,2009). Clear biases arise due to the node degree, a topological property of the nodes: the larger the pathway, the higher its mean value, and the more connected a compound is, the smaller its variance. If pathways are ranked by raw temperatures, a large pathway will have an undesired, consistent advantage over small ones and will be reported too often. The usage of z-scores (d) instead of raw temperatures (c) to select the top 250 nodes addresses these biases and highlights pathway and module nodes that were eclipsed by other compounds and reactions with higher mean null temperatures. the mean temperatures from the null model (Fig 20a), which is a property of the graph but not of the experimental data; using z-scores instead corrects this bias (Fig 20b). If the top 20 pathways are selected through their raw temperature, some of them are even below their mean null temperature (Fig 20c), whereas keeping the top 20 z-scores removes the bias towards larger pathways and suggests otherwise overlooked pathways (Fig 20d).
6.3 results 87 Figure 20: Ranking the 288 KEGG pathways – lower is best– using raw temperatures (a) biases the ranks towards pathways with higher mean null temperature, which in turn tend to be large pathways. Using the z-scores instead (b) breaks this clear dependence and avoids reporting pathways just because of their size. The top 20 pathways through raw temperatures (c), depicted as black dots, include pathways that are even below their mean value, while the top 20 z-scores (d) suggest smaller pathways that were penalised by the aforementioned bias. 6.3.3 Subgraph extraction Four subgraphs have been extracted using the MS-derived compounds. The desired number of nodes kfor each approach, together with the actual number of reported nodes and the number of KEGG pathways, are shown in Table 6. A connected component (CC) of an undirected graph is a maximal connected subgraph so that any two nodes in the subgraph are connected by a path. For the directed graphs, the weak CC definition is used, in which directed edges are considered as undirected when computing the CC. The number of nodes belonging to each solution subgraph, along with its largest CC and the number of CCs, are also reported. Additional details regarding the largest CC and number of CCs for other values of kcan be found in S5 Appendix. Defining the overlap coefficient between two solutions G1and G2as
88 metabolomics enrichment Table 6: Summary of the outputs Name k Pathways Nodes #CC Largest CC HD norm 250 hsa00250, hsa00270, hsa00480, hsa05230, hsa05231 250 8 206 HD sim 250 hsa00250, hsa00270, hsa00330, hsa00480, hsa05230, hsa05231 261 8 221 PR norm 250 hsa00250, hsa00270, hsa00480, hsa05231 250 9 187 PR sim 250 hsa00250, hsa00270, hsa00480, hsa05231 279 10 152 Summary of the outputs, using diffusion (HD) as well as PageRank (PR), and normalising the scores with Monte Carlo simulations (sim) or z-scores (norm). Monte Carlo simulations have been run 10,000 times per solution, and 9solutions have been computed to build a consensus solution. Note that the desired number of nodes kis slightly different to the number of nodes actually reported in the Monte Carlo simulations. The last two columns contain the number of connected components (CC) and the number of nodes in the largest CC. overlap(G1,G2) = |G1∩G2| min(|G1|,|G2|), (41) solutions tend to overlap despite their differences (Table 7). Regarding the stratification of the subgraphs in terms of KEGG categories, they follow a trend similar to the KEGG graph (S5Appendix). Table 7: Solutions overlap HD norm HD sim PR norm PR sim HD norm 1.00 0.82 0.88 0.82 HD sim 0.82 1.00 0.77 0.83 PR norm 0.88 0.77 1.00 0.84 PR sim 0.82 0.83 0.84 1.00 Overlap coefficient statistics for HD and PR. The overlapping nature of solutions is a sign of consistency among approaches. 6.3.4 Pathway analysis Our methods are compared to IMPaLA and MetaboAnalyst to verify the concordance in terms of metabolic pathways. All the approaches have been compared using the example data from IMPaLA (S2Table) and MetaboAnalyst (S3Table), and they show consistent and compatible reports. The results for our dataset are summarised in Table 8and described in S1Table, together with further details about the reports of the alternative tools. The metabolic pathways Alanine, aspartate and glutamate metabolism (hsa00250), Cysteine and methionine metabolism (hsa00270) and especially the Glutathione metabolism (hsa00480) recur in all of the approaches. Some of our solutions are more specific, suggesting the module Glutathione Biosynthesis (M00118) as well. Our null model takes pathway overlap and crosstalk into account and allows a visualisation of the pathway structure through the null diffusion correlation matrix (S4Appendix). The subgraph resulting from applying HD sim (Fig 21) inherits the scalefree structure from the whole graph and enrols the three recurrently reported pathways in the same connected component: hsa00250, hsa00270 and hsa00480. The biological perturbation stemming from the MS-derived com-
6.3 results 89 Table 8: Reported pathways KEGG id Pathway name HD norm HD sim PR norm PR sim MA FCS MA ORA IMPaLA ORA hsa00250 Alanine, aspartate and glutamate metabolism + + + + + + - hsa00270 Cysteine and methionine metabolism + + + + + + + hsa00480 Glutathione metabolism + + + + + + + hsa05230 (hsa00970) Central carbon metabolism in cancer + + - - * - + hsa05231 (hsa00564) Choline metabolism in cancer + + + + * - - hsa00260 (M00020) Glycine, serine and threonine metabolism * * - - + - - hsa00330 (M00133) Arginine and proline metabolism * + - - + - + hsa00510 (M00073) N-Glycan biosynthesis - - * * - - - Pathways reported by our methods. ’+’ means a hit for the term reported in the KEGG id column, ’*’ stands for a hit of the closely related term in parenthesis in the same column and ’-’ states no hit. Our 4solutions are compared to MetaboAnalyst (MA), using ORA and FCS, and IMPaLA using ORA. Pathways hsa00250, hsa00270 and hsa00480 are repeatedly reported by all the methodologies. Pathways hsa05230 and hsa05231 are reported by some of our methods, while alternative approaches find some close (*) and exact (+) matches. In some cases, instead of reporting a whole pathway, only specific modules within it are reported as relevant; this is the case of M00133 and M00073. Furthermore, module M00073 does not contain any compounds, being out of the scope of MetaboAnalyst and IMPaLA, but is reported by one of our methods due to the presence of other indirect relationships through enzymes in the graph. pounds can be tracked in terms of reactions, enzymes and modules, up to the relevant pathways. On the other hand, results on the recovery of synthetic signals can be found in Fig 22. In (a) absence of signal, HD ranks pathways with a mean rank close to 0.5, and only a few are biased to the top or the bottom of the list. Mean ranks in Fisher’s exact test and PR are also centered around 0.5, but have more dispersion. In (b) the presence of a target pathway, three sampling schemes have been explored. In (1) the signal is actually noise and the target pathway is a decoy. The rank of the target pathway for HD and PR is uniformly spread in [0,1], whereas Fisher’s exact test shows some asymmetry in the rank distribution. In (2), the sampling probability depends on the presence or absence of the metabolite in the pathway. Fisher’s exact test outperforms HD and PR as the median rank of the target pathway is closer to 0, as expected by its optimality. However, in (3), the sampling probability is network-based and HD outperforms PR, which in turn outperforms Fisher’s exact test. Differences between sim (Monte Carlo trials) and norm (parametric approach) are subtle. 6.3.5 NMR analysis NMR carbon tracking revealed 13 isotopically enriched metabolites from 13C-glucose, showing differential fractional enrichment between case-control, of which 5had already been found through MS; some of these metabolites can be seen in Fig 23 in the context of the Glutathione metabolism. Our solutions are assessed in terms of the resistance distance from the reported reactions to the remaining 8metabolites. The smaller the overall distance of a solution, the more related its nodes are to the 8metabolites proven affected by NMR. The resistance distances have been computed on the reaction-compound graph, which is the largest CC of the subgraph that contains all the reactions and compounds in the KEGG graph. The reactions suggested in our subgraphs show lower resistance distances to the 8NMR-derived metabolites than the totality of reactions in the reactioncompound graph (Table 9). Furthermore, they are also lower than the resis-
90 metabolomics enrichment Top 250 z-scores for diffusion m e N-Acetyl-L-glutamate P lic u fide G a L-2-Aminoadipate -Acetylaspartylglu e te 4-Aminob anoate Cys-Gly a-L-Glu am ne st L-Cystine Alanine L-Proline Glyc L mocysteine FAD Sphingomyelin 2-Hydroxyglutarate GABA (gamm u rat Polyamine biosynthesis, arginine => agmatine => utrescine = spermidine Cysteine biosynthesis, homocysteine + serine => cysteine Glutathione biosyn ate utathi Serine biosynthesis, glycerate-3P => serine Module Enzyme Reaction Compound Input compound Pathway Legend Icosadienoic acid Pantothenate (11Z,14Z)-Icosadienoyl-CoA N-Acetylneuraminate D-Glucose D-Glucose 1-phosphate UPD-Glucuronate 3-(4-Hydroxyphenyl)lactate N-Formylmethionine 5'-Methylthioadenosine sn-Glycero phate sn-Glycerol-3-phosphocholine myo-In tral carbon metabolism in cancer Cysteine and onine metabo n m Choline metabolism in cancer A an e, and g uta ism N-Ace rtate Figure 21: Subgraph reported through HD norm, the names of reactions and enzymes have been omitted for clarity. Compounds are green, reactions are blue, enzymes are orange, modules are purple and pathways are red. The compounds in the input are highlighted as green squares to ease the tracing of the biological perturbation up to the pathways. The presence of reactions and enzymes that link pathways in this subgraph might suggest relevant entities by which affected pathways crosstalk. All the reported pathways and modules lie in a large CC, as well as a newly proposed metabolite (L-Glutamate). tance distances from the neighbouring reactions of the MS-derived metabolites to the 8NMR metabolites (FDR < 0.01). 6.4 discussion Our approach for enriching summary metabolomics data, Fig 16, is based on diffusion processes over a graph drawn from several KEGG categories (Fig 17). KEGG is the database of choice due to its level of curation and structure, which eases the graph representation. Specifically, the definition of KEGG categories naturally allows a hierarchical arrangement of levels. Our graph design is enhanced by the compound-reaction-enzyme-gene networks built by MetScape (S1Appendix), and the inclusion of modules and
6.4 discussion 91 Figure 22: Synthetic signals evaluation using the pathway rank as a metric to assess orderings. Lowest ranks correspond to best ranked pathways. The proposed methodology is compared to ORA, represented by Fisher’s exact test. (a) 288 noisy signals have been generated, and every pathway has been ranked in each of the 288 runs. Data points for a given methodology are the mean rank of each pathway, giving 288 data points per box. (b) 288 signals with a target pathway have been generated, in three scenarios: pure noise, proportion-based sampling and network-based sampling. Each box contains the rank of the target pathway, leading to 288 data points per box. Table 9: Distance to NMR metabolites Method Graph order C00299 C00122 C00116 C00105 C00020 C00581 C00300 C00025 Reaction-compound graph 4539[8008]0.56(0.62)0.56(0.62)0.57(0.62)0.54(0.62)0.47(0.62)0.93(0.62)0.82(0.62)0.47(0.62) First neighbours 414[447]0.42(0.12)0.43(0.12)0.44(0.12)0.40(0.12)0.33(0.12)0.79(0.12)0.68(0.12)0.33(0.12) HD norm 147[250]0.39(0.10)0.39(0.10)0.40(0.10)0.37(0.10)0.30(0.10)0.76(0.10)0.65(0.10)0.30(0.10) HD sim 148[261]0.39(0.09)0.39(0.09)0.40(0.10)0.37(0.09)0.30(0.09)0.76(0.09)0.65(0.09)0.30(0.09) PR norm 143[250]0.39(0.10)0.39(0.10)0.40(0.10)0.37(0.10)0.30(0.10)0.75(0.10)0.65(0.10)0.30(0.10) PR sim 172[279]0.40(0.12)0.41(0.12)0.42(0.12)0.38(0.12)0.31(0.12)0.77(0.12)0.66(0.12)0.31(0.12) Mean resistance distance between the reactions reported in our solutions and each compound reported using NMR, with their standard deviations in parentheses. For each subgraph of KEGG graph, the number of reactions and the total number of nodes (in square brackets) are displayed. The reactioncompound subgraph contains the largest connected component having all the reactions and compounds in the KEGG graph. The first neighbours subgraph contains the MS-derived metabolites and all the reactions in which they participate. Resistance distances are computed on the reaction-compound graph. For every NMR-derived metabolite, there is a significant difference in resistance distances between the reactions proposed in our solutions and the reactions involving any of the MS-derived metabolite (one-sided Wilcoxon test, FDR < 0.01 for the 32 possible comparisons: 8NMR metabolites, tests of 4solutions against the first neighbours reactions). This implies that the reported reactions are closer to the NMR- derived compounds than the bulk of neighbouring reactions. pathways in our arrangement allows a comprehensive picture of the affected biology. The graph contains all the KEGG compounds and the subset of affected metabolites forced to diffuse inside it (Fig 17). The closer a node is to the affected compounds, the higher its score becomes. Likewise, the top scoring candidates naturally involve higher flow and become relevant in the flow discharge from the graph. Because our KEGG graph is conceived and curated in a bottom-up manner, diffusion is expected to follow that trend too: the perturbation in the lowest level will diffuse to the upper levels to exit the graph. Ideally, a relevant subgraph found through this diffusion (Fig 21) would inherit the stratification of the KEGG graph, thus allowing the extrapolation of knowledge in terms of compounds to the rest of categories. This allows holistic picturing of pathways of interest, such as Glutathione
98 metabolomics enrichment references Aivio, Suvi Marjaana and Travis H. Stracker 2014 The Role of EXD2in the maintenance of mithocondrial homeostasis, Doctoral Thesis, Universitat Pompeu Fabra, Departament de Ciències Experimentals i de la Salut. Alonso, Arnald, Sara Marsal, and Antonio Julià 2015 “Analytical methods in untargeted metabolomics: state of the art in 2015”, Front. Bioeng. Biotechnol.,3,23,issn:2296-4185. Bapat, RB 2004 “Resistance matrix of a weighted graph”, MATCH-COMMUN MATH CO,50, pp. 73-82. Bonals, Lluís Albert 2005 Transferència de calor: apunts de classe. Chagoyen, Monica and Florencio Pazos 2011 “MBRole: enrichment analysis of metabolomic data”, Bioinformatics, 27,5, pp. 730-731. 2013 “Tools for the functional interpretation of metabolomic experiments”, Brief. Bioinform.,14,6, pp. 737-744. Croft, David, Antonio Fabregat Mundo, Robin Haw, Marija Milacic, Joel Weiser, Guanming Wu, Michael Caudy, Phani Garapati, Marc Gillespie, Maulik R. Kamdar, Bijay Jassal, Steven Jupe, Lisa Matthews, Bruce May, Stanislav Palatnik, Karen Rothfels, Veronica Shamovsky, Heeyeon Song, Mark Williams, Ewan Birney, Henning Hermjakob, Lincoln Stein, and Peter D’Eustachio 2014 “The Reactome pathway knowledgebase”, Nucleic Acids Res.,42, Database issue, pp. D472-D477. Csardi, Gabor and Tamas Nepusz 2006 “The igraph software package for complex network research”, InterJournal, Complex Systems, p. 1695,http://igraph.org. DeBerardinis, Ralph J, Anthony Mancuso, Evgueni Daikhin, Ilana Nissim, Marc Yudkoff, Suzanne Wehrli, and Craig B Thompson 2007 “Beyond aerobic glycolysis: transformed cells can engage in glutamine metabolism that exceeds the requirement for protein and nucleotide synthesis”, Proc. Natl. Acad. Sci. U.S.A.,104,49, pp. 19345- 19350. Donato, Michele, Zhonghui Xu, Alin Tomoiaga, James G Granneman, Robert G MacKenzie, Riyue Bao, Nandor Gabor Than, Peter H Westfall, Roberto Romero, and Sorin Draghici 2013 “Analysis and correction of crosstalk effects in pathway analysis”, Genome Res.,23,11, pp. 1885-1893.
References 99 Draghici, Sorin, Purvesh Khatri, Adi Laurentiu Tarca, Kashyap Amin, Arina Done, Calin Voichita, Constantin Georgescu, and Roberto Romero 2007 “A systems biology approach for pathway level analysis”, Genome Res.,17,10, pp. 1537-1545. Faust, Karoline, Pierre Dupont, Jérôme Callut, and Jacques van Helden 2010 “Pathway discovery in metabolic networks by subgraph extraction”, Bioinformatics,26,9, pp. 1211-1218,issn:13674803. Fernández-Albert, Francesc, Rafael Llorach, Cristina Andrés-Lacueva, and Alexandre Perera 2014 “An R package to analyse LC/MS metabolomic data: MAIT (Metabolite Automatic Identification Toolkit)”, Bioinformatics,30,13, pp. 1937- 1939. Haynes, Winston A, Roger Higdon, Larissa Stanberry, Dwayne Collins, and Eugene Kolker 2013 “Differential expression analysis for pathways”, PLOS Comput. Biol., 9,3, e1002967. Huang, Da Wei, Brad T Sherman, and Richard A Lempicki 2009 “Bioinformatics enrichment tools: paths toward the comprehensive functional analysis of large gene lists”, Nucleic Acids Res.,37,1, pp. 1-13. Ideker, Trey, Owen Ozier, Benno Schwikowski, and Andrew F Siegel 2002 “Discovering regulatory and signalling circuits in molecular interaction networks”, Bioinformatics,18, suppl 1, S233-S240. Kamburov, Atanas, Rachel Cavill, Timothy M. D. Ebbels, Ralf Herwig, and Hector C. Keun 2011 “Integrated pathway-level analysis of transcriptomics and metabolomics data with IMPaLA”, Bioinformatics,27,20, pp. 2917-2918. Kanehisa, Minoru, Michihiro Araki, Susumu Goto, Masahiro Hattori, Mika Hirakawa, Masumi Itoh, Toshiaki Katayama, Shuichi Kawashima, Shujiro Okuda, Toshiaki Tokimatsu, and Yoshihiro Yamanishi 2008 “KEGG for linking genomes to life and the environment”, Nucleic Acids Res.,36, Database issue, pp. D480-D484. Kankainen, Matti, Peddinti Gopalacharyulu, Liisa Holm, and Matej Orešiˇc 2011 “MPEA – metabolite pathway enrichment analysis”, Bioinformatics, 27,13, pp. 1878-1879. Karnovsky, Alla, Terry E. Weymouth, Tim Hull, V. Glenn Tarcea, Giovanni Scardoni, Carlo Laudanna, Maureen A. Sartor, Kathleen A. Stringer, H. V. Jagadish, Charles F. Burant, Brian D. Athey, and Gilbert S. Omenn 2012 “Metscape 2bioinformatics tool for the analysis and visualization of metabolomics and gene expression data”, Bioinformatics,28,3, pp. 373-380.
100 metabolomics enrichment Kelder, Thomas, Martijn P van Iersel, Kristina Hanspers, Martina Kutmon, Bruce R Conklin, Chris T Evelo, and Alexander R Pico 2012 “WikiPathways: building research communities on biological pathways”, Nucleic Acids Res.,40, Database issue, pp. D1301-D1307. Kessler, Nikolas, Heiko Neuweger, Anja Bonte, Georg Langenkämper, Karsten Niehaus, Tim W. Nattkemper, and Alexander Goesmann 2013 “MeltDB 2.0-advances of the metabolomics software system”, Bioinformatics,29,19, pp. 2452-2459. Khatri, Purvesh, Marina Sirota, and Atul J. Butte 2012 “Ten Years of Pathway Analysis: Current Approaches and Outstanding Challenges”, PLOS Comput. Biol.,8,2. Li, Xianbin, Liangzhong Shen, Xuequn Shang, and Wenbin Liu 2015 “Subpathway analysis based on signaling-pathway impact analysis of signaling pathway”, PLOS ONE,10,7, e0132813. Mitra, Koyel, Anne-Ruxandra Carvunis, Sanath Kumar Ramesh, and Trey Ideker 2013 “Integrative approaches for finding modular structure in biological networks”, Nat. Rev. Genet.,14,10, pp. 719-732. Nicholson, Jeremy K., John Connelly, John C. Lindon, and Elaine Holmes 2002 “Metabonomics: a platform for studying drug toxicity and gene function”, Nat. Rev. Drug Discov.,1,2, pp. 153-161. North, Bernard V, David Curtis, and Pak C Sham 2002 “A note on the calculation of empirical P values from Monte Carlo procedures”, Am. J. Hum. Genet.,71,2, p. 439. Page, Lawrence, Sergey Brin, Rajeev Motwani, and Terry Winograd 1999 The PageRank citation ranking: bringing order to the Web, tech. rep., Stanford InfoLab. Paull, Evan O, Daniel E Carlin, Mario Niepel, Peter K Sorger, David Haussler, and Joshua M Stuart 2013 “Discovering causal pathways linking genomic events to transcriptional states using Tied Diffusion Through Interacting Events (TieDIE)”, Bioinformatics,29,21, pp. 2757-2764. R Core Team 2015 R: A Language and Environment for Statistical Computing, R Foundation for Statistical Computing, Vienna, Austria, https :// www.R - project.org/. Rahnenführer, Jörg, Francisco S Domingues, Jochen Maydt, and Thomas Lengauer 2004 “Calculating the statistical significance of changes in pathway activity from gene expression data”, Stat. Appl. Genet. Mol.,3,1. Reddy, Junuthula Narasimha and David K Gartling 2010 The finite element method in heat transfer and fluid dynamics.
References 101 Smoot, Michael E., Keiichiro Ono, Johannes Ruscheinski, Peng-Liang Wang, and Trey Ideker 2011 “Cytoscape 2.8: new features for data integration and network visualization”, Bioinformatics,27,3, pp. 431-432. Subramanian, Aravind, Pablo Tamayo, Vamsi K. Mootha, Sayan Mukherjee, Benjamin L. Ebert, Michael A. Gillette, Amanda Paulovich, Scott L. Pomeroy, Todd R. Golub, Eric S, Lander, and Jill P. Mesirov 2005 “Gene set enrichment analysis: A knowledge-based approach for interpreting genome-wide expression profiles”, Proc. Natl. Acad. Sci. U.S.A.,102,43, pp. 15545-15550. Tarca, Adi Laurentiu, Sorin Draghici, Gaurav Bhatti, and Roberto Romero 2012 “Down-weighting overlapping genes improves gene set analysis”, BMC Bioinform.,13,1, p. 136. Vandin, Fabio, Eli Upfal, and Benjamin J Raphael 2011 “Algorithms for detecting significantly mutated pathways in cancer”, J. Comput. Biol.,18,3, pp. 507-522. Vinaixa, Maria, Emma L Schymanski, Steffen Neumann, Miriam Navarro, Reza M Salek, and Oscar Yanes 2015 “Mass spectral databases for LC/MS and GC/MS-based metabolomics: state of the field and future prospects”, TrAC-Trend. Anal. Chem.,78, pp. 23-25. Weckwerth, Wolfram 2003 “Metabolomics in Systems Biology”, Annu. Rev. Plant Biol.,54,1, pp. 669-689. Wickham, Hadley 2009 ggplot2: Elegant Graphics for Data Analysis, Springer-Verlag New York, isbn:978-0-387-98140-6,http://ggplot2.org. Wishart, David S., Timothy Jewison, Anchi Guo, Michael Wilson, Craig Knox, Yifeng Liu, Yannick Djoumbou, Rupasri Mandal, Farid Aziat, Edison Dong, et al. 2013 “HMDB 3.0- The Human Metabolome Database in 2013”, Nucleic Acids Res.,41, Database issue, pp. D801-D807. Xia, Jianguo, Igor V Sinelnikov, Beomsoo Han, and David S Wishart 2015 “MetaboAnalyst 3.0– making metabolomics more meaningful”, Nucleic Acids Res.,43, Web Server issue, W251-W257. Xia, Jianguo and David S. Wishart 2010 “MSEA: a web-based tool to identify biologically meaningful patterns in quantitative metabolomic data”, Nucleic Acids Res.,38, Web Server issue, W71-W77.
7THE R PACKAGE FELLA fella: an r package to enrich metabolomics data Pathway enrichment techniques are useful for understanding experimental metabolomics data. Their purpose is to give context to the affected metabolites in terms of the prior knowledge contained in metabolic pathways. However, the interpretation of a prioritized pathway list is still challenging, as pathways show overlap and cross talk effects. We introduce FELLA, an R package to perform a network-based enrichment of a list of affected metabolites. FELLA builds a hierarchical representation of an organism biochemistry from the Kyoto Encyclopedia of Genes and Genomes (KEGG), containing pathways, modules, enzymes, reactions and metabolites. In addition to providing a list of pathways, FELLA reports intermediate entities (modules, enzymes, reactions) that link the input metabolites to them. This sheds light on pathway cross talk and potential enzymes or metabolites as targets for the condition under study. FELLA has been applied to six public datasets –three from Homo sapiens, two from Danio rerio and one from Mus musculus– and has reproduced findings from the original studies and from independent literature. The R package FELLA offers an innovative enrichment concept starting from a list of metabolites, based on a knowledge graph representation of the KEGG database that focuses on interpretability. Besides reporting a list of pathways, FELLA suggests intermediate entities that are of interest per se. Its usefulness has been shown at several molecular levels on six public datasets, including human and animal models. The user can run the enrichment analysis through a simple interactive graphical interface or programmatically. FELLA is publicly available in Bioconductor under the GPL-3 license. 7.1 background Metabolomics is the science that measures lightweight molecules in living organisms and stands as a valuable source of biomarkers and biological knowledge (Madsen et al.,2010). The preprocessing of such data can be achieved through pipelines like MeltDB (Kessler et al.,2013) or MAIT (Fernández-Albert et al.,2014). Once metabolite abundances are available, pathway analysis tools ease data interpretation (Khatri et al.,2012) by framing the affected metabolites in terms of contextual knowledge. Databases This chapter is a postprint of the following journal article: Picart-Armada, Sergio, Francesc Fernández-Albert, Maria Vinaixa, Oscar Yanes, and Alexandre Perera-Lluna. “FELLA: an R package to enrich metabolomics data”. BMC bioinformatics 19, no. 1(2018): 538. 103
104 the r package fella like the Kyoto Encyclopedia of Genes and Genomes (KEGG) (Kanehisa et al., 2011) are sources of curated pathway data. The classification of enrichment techniques used here follows the review in (Khatri et al.,2012). Over representation analysis (ORA) approaches are based on testing the proportion of a list of affected metabolites inside a pathway. ORA is available in tools like the web server MetaboAnalyst (Xia et al.,2015) and the R package clusterProfiler (Yu et al.,2012). Functional class scoring (FCS) approaches use quantitative data instead and seek subtle but coordinated changes in the metabolites belonging to a pathway. MSEA in MetaboAnalyst and IMPaLA (Kamburov et al.,2011) contain implementations of FCS for metabolomics. Pathway topology-based (PT) approaches further include topological measures of the metabolites in the statistic, accounting for their inequivalence within the pathway. PT analyses can be performed using MetaboAnalyst. Here, we introduce the R package FELLA, available in Bioconductor (Huber et al.,2015), for metabolomics data interpretation that combines pathway enrichment with network analysis. The list of affected metabolites and the reported pathways are connected through intermediate entities -reactions, enzymes, modulesin a heterogeneous network layout. This suggests how the perturbation spreads at the pathway level and how pathways cross talk, enhancing the interpretability of the output. 7.2 implementation FELLA is an R package that performs metabolomics data enrichment starting from (I) a network derived from KEGG and (II) a list of KEGG compounds (Fig. 25). A sub-network relevant to the input is extracted from (I) using network propagation algorithms that start from the labels in (II), providing a data enrichment that goes beyond a pathway list. The purpose of FELLA is to elaborate a biological explanation that justifies how the input metabolites can reach the reported pathways, as well as perspective on pathway cross talk. Two user guides illustrate the principles and the usage of FELLA: a quickstart (additional file 3) and an in-depth vignette with implementations details and three real examples (additional file 1). Two additional vignettes (additional files 4and 6) serve as case studies for non-human organisms. 7.2.1 Methodology The cornerstone of FELLA is its knowledge graph representation of the biochemistry in KEGG at several molecular levels. The network is hierarchical and connects KEGG compounds (metabolites) to KEGG pathways through intermediate entities, namely reactions, enzymes and KEGG modules, see figure 24. Such connections (edges) are obtained directly from KEGG annotations. The presence of intermediate levels allows inference at their level, meaning that relevant reactions, enzymes and KEGG modules can be suggested just by starting from a list of affected metabolites. This feature is
7.2 implementation 105 evaluated in several case studies, by linking the suggested enzymatic families and reactions to literature and to original findings within the studies. Compounds Reactions Enzymes Modules Pathways (a) (b) Figure 24: Node arrangement for the knowledge model used by FELLA. Entities are organised in a hierarchical manner, from bottom to top: KEGG compounds or metabolites, reactions, enzymes, KEGG modules and pathways. Binary labels at the level of metabolites are propagated to the rest of the network and a relevant, small sub-network is automatically reported. Nodes are ranked using the network propagation algorithms (a) heat diffusion and (b) PageRank. The affected metabolites are highlighted with a black ring. For heat diffusion (a), affected metabolites are forced to generate unitary flow. Every pathway is highlighted with a blue ring, representing its connection to a cool boundary node. In equilibrium, the highest temperature pathways (and nodes) will have the greatest heat flow, suggesting a relevant role in the experiment. For PageRank (b), affected metabolites are the start of random walks. PageR- ank scores, represented by the intensity of the blue colour, will attain higher values in the frequently reached random walk nodes. Figure extracted from (Picart-Armada et al.,2017). In order to report a sub-network, nodes are ranked according to a scoring function –based on network propagation– and only the top scoring nodes are returned. Two algorithms are supported for propagating the labels from the affected metabolites: a classical heat diffusion approach (Vandin et al.,2011) and the PageRank web ranking algorithm (Page et al.,1999). Further details on the network propagation settings can be found in (Picart-Armada et al., 2017) and in additional file 1. The main difference between both algorithms is that heat diffusion is undirected whereas PageRank is directed upwards. In practice, contrary to PageRank, heat diffusion will frequently report new metabolites because heat is allowed to propagate back to compounds from the upper levels (Picart-Armada et al.,2017). This behaviour can ease the discovery of intermediate metabolites that lay close to the input metabolites and tend to connect them. An example of its usefulness can be found in the gilt-head bream study. As exposed in (Picart-Armada et al.,2017), ranking nodes according to their raw diffusion scores suffers from a strong bias, related to the node level and topological features. This is addressed by normalising the diffusion score of every node using its background distribution under input permutations. Permutations can be simulated through Monte Carlo trials to obtain an empirical p-value, labelled as p-score. Alternatively, a parametric z-score can be obtained without requiring Monte Carlo trials. The p-score is obtained by transforming the z-score to lie in the [0,1]interval through the cumulative distribution function of a standard normal distribution. Under both statistical approximations, nodes with the lowest p-scores are reported
106 the r package fella as the suggested sub-network. Note that p-scores are used as a ranker rather than for testing hypotheses. An optional filter allows the removal of small connected components from the reported sub-network. When building the database, a number of random sub-networks are sampled to characterise how infrequent a connected component of order at least ris when knodes are uniformly sampled. The assumption behind this filter is that meaningful inputs encompass metabolites relatively close to each other within the knowledge graph, prone to be reported in large connected components involving most of them. 7.2.2 Classes FELLA relies on two classes: FELLA.DATA for the internal knowledge representation, based on the igraph R package (Csardi and Nepusz,2006), and FELLA.USER for the user analysis, see figure 25. These classes contain subclasses, invisible to the user and described in the additional file 1. The functions to manipulate both classes are described below, following the three blocks from figure 25. buildGraphFromKEGGREST() buildDataFromGraph() loadKEGGdata() defineCompounds() runHypergeom() runDiffusion() runPagerank() generateResultsTable() generateEnzymesTable()generateResultsGraph() addGOtoGraph() plotGraph() Organism e.g. Homo sapiens Input compounds e.g. C1, C3, C8, C9, C24 C1 C3 C8 FELLA.USER C1 C3 C8 FELLA.USER 0.5 0.1 1.2 . . . FELLA.DATA III III Figure 25: Design of the R package FELLA.(I) creation of a graph object from an organism code and its database, (II) ID mapping and propagation algorithms (diffusion, PageRank) to score all the nodes, (III) node prioritisation and results exporting.
7.2 implementation 107 Block I: local database The function buildGraphFromKEGGREST() retrieves the tabular KEGG data for the desired organism and builds the knowledge graph as described in (Picart-Armada et al.,2017). Then, a database can be built from the graph and stored in a local folder using buildDataFromGraph(). Databases are needed for the enrichment and should be loaded through the function loadKEGGdata(). Block II: enrichment analysis Once the database is loaded, i.e. the FELLA.DATA object is in memory, defineCompounds() maps the list of input metabolites, in the form of KEGG identifiers, to the internal representation, providing a FELLA.USER object. Then, the propagation algorithms in (Picart-Armada et al.,2017) are run to score the graph nodes. runDiffusion() uses the undirected heat diffusion model (Vandin et al.,2011) whereas runPagerank() runs the directed PageRank algorithm (Page et al.,1999). Both approaches are automatically followed by the statistical normalisation, either as a parametric z-score (approx = "normality") or as a simulated permutation analysis (approx = "simulation"), see table 10. The wrapper enrich() performs the metabolite mapping and the desired propagation algorithm (argument method) and statistical normalisation with a single call. Table 10: Scoring methods offered in FELLA, chosen by the enrich function arguments method and approx. Each row corresponds to a method mentioned in the original publication (Picart-Armada et al.,2017). The method hypergeom is Fisher’s exact test, included for reference. Method diffusion scores the nodes using the heat diffusion model to score the nodes. Method pagerank uses the PageRank algorithm on an upwardsdirected version of the network. Both scores undergo a statistical normalisation to remove structural biases, controlled through the approx argument. The user can choose the fast, parametric z-scores (normality) or the slower, non-parametric permutation analysis (simulation). N/A: non-applicable. Method Approx Notation in (Picart-Armada et al.,2017) Comment hypergeom N/A hypergeometric test Included for reference diffusion normality HD norm Heat diffusion scores followed by z-scores diffusion simulation HD sim Heat diffusion scores followed by permutations pagerank normality PR norm PageRank scores followed by z-scores pagerank simulation PR sim PageRank scores followed by permutations Block III: exporting results Finally, the best scoring KEGG entries can be visualised through plot(), exported as a sub-network with generateResultsGraph(), or in tabular format with generateResultsTable(). A dedicated table with the reported enzymes and its associated genes can be obtained with generateEnzymesTable(). Alternatively, exportResults() allows writing such objects directly to files.
114 the r package fella declarations ethics approval and consent to participate Not applicable. consent to publish Not applicable. availability of data and materials All data generated or analysed during this study are included in this published article (additional files 2,5and 7). competing interests The authors declare that they have no competing interests. funding This work was supported by the Spanish Ministry of Economy and Competitiveness (MINECO) [BFU 2014-57466-P to OY, TEC 2014-60337-R and DPI 2017-89827-R to AP]. OY, AP and SP thank for funding CIBERDEM and CIBER-BBN, both initiatives of Instituto de Investigación Carlos III (ISCIII). SP thanks the AGAUR FI-scholarship programme. The funding body had no role in the design of the study and collection, analysis, and interpretation of data and in writing the manuscript. authors’ contributions SP, FF, MV, OY and AP conceived the software. SP implemented the software and analysed the data. SP wrote the original manuscript. FF, MV, OY and AP critically revised the original manuscript. OY and AP supervised the project. All authors read and approved the final manuscript. acknowledgements We would like to thank Haizea Ziarrusta and our collaboration with the Department of Analytical Chemistry, University of the Basque Country (UP- V/EHU), Leioa, for using, discussing and helping improve our software.
7.4 conclusions 115 We would also like to thank the anonymous reviewers for their valuable comments. additional files Additional file 1 — FELLA.pdf User guide within the R package FELLA with background, implementation details and three real examples on its usage. Additional file 2 — datasets.zip Descriptive files on the three human datasets: a summary of the inputs (descriptive_input.csv), input and reported subgraph in each dataset (dataset_input.csv, dataset_subgraph.csv and dataset_subgraph.pdf), hits discussed in the results section (descriptive_hits.csv). Also contains the database object (fella_data.RData) and metadata about the database (info_fella_data.txt), the KEGG version (info_kegg.txt) and the R session (info_session.txt). Additional file 3 — quickstart.html User guide within FELLA showing fast and concise toy examples of its application. Additional file 4 — zebrafish.pdf Case study with FELLA: two datasets on the effect of oxybenzone exposition on gilt-head bream. Additional file 5 — zebrafish.zip R workspace from the gilt-head bream datasets. Additional file 6 — musmusculus.pdf Case study with FELLA: a multi-omic mouse model of non-alcoholic fatty liver disease. Additional file 7 — musmusculus.zip R workspace from the mouse model study.
116 the r package fella references Beisel, William R 1975 “Metabolic response to infection”, Annu. Rev. Med.,26,1, pp. 9-20. Chang, Winston, Joe Cheng, JJ Allaire, Yihui Xie, and Jonathan McPherson 2018 shiny: Web Application Framework for R, R package version 1.1.0,htt ps://CRAN.R-project.org/package=shiny. Chen, Liyan, Jing Li, Tiannan Guo, Sujoy Ghosh, Siew Kwan Koh, Dechao Tian, Liang Zhang, Deyong Jia, Roger W Beuerman, Ruedi Aebersold, et al. 2015 “Global metabonomic and proteomic analysis of human conjunctival epithelial cells (IOBA-NHC) in response to hyperosmotic stress”, J. Proteome Res.,14,9, pp. 3982-3995. Consortium, Gene Ontology et al. 2015 “Gene ontology consortium: going forward”, Nucleic Acids Res.,43, D1, pp. D1049-D1056. Csardi, Gabor and Tamas Nepusz 2006 “The igraph software package for complex network research”, InterJournal, Complex Systems, p. 1695. Decuypere, Saskia, Jessica Maltha, Stijn Deborggraeve, Nicholas JW Rattray, Guiraud Issa, Kaboré Bérenger, Palpouguini Lompo, Marc C Tahita, Thusitha Ruspasinghe, Malcolm McConville, et al. 2016 “Towards Improving Point-of-Care Diagnosis of Non-malaria Febrile Illness: A Metabolomics Approach”, PLoS Negl Trop Dis,10,3, e0004480. Fernández-Albert, Francesc, Rafael Llorach, Cristina Andrés-Lacueva, and Alexandre Perera 2014 “An R package to analyse LC/MS metabolomic data: MAIT (Metabolite Automatic Identification Toolkit)”, Bioinformatics,30,13, pp. 1937- 1939. Fitch, Coy D, Guang-zuan Cai, and James D Shoemaker 2000 “A role for linoleic acid in erythrocytes infected with Plasmodium berghei”, Biochim. Biophys. Acta-Mol. Basis Dis.,1535,1, pp. 45-49. Godoy, Patricio, Agata Widera, Wolfgang Schmidt-Heck, Gisela Campos, Christoph Meyer, Cristina Cadenas, Raymond Reif, Regina Stöber, Seddik Hammad, Larissa Pütter, et al. 2016 “Gene network activity in cultivated primary hepatocytes is highly similar to diseased mammalian liver tissue”, Arch. Toxicol.,90,10, pp. 2513-2529.
References 117 Gogiashvili, Mikheil, Karolina Edlund, Kathrin Gianmoena, Rosemarie Marchan, Alexander Brik, Jan T Andersson, Jörg Lambert, Katrin Madjar, Birte Hellwig, Jörg Rahnenführer, et al. 2017 “Metabolic profiling of ob/ob mouse fatty liver using HR-MAS 1H-NMR combined with gene expression analysis reveals alterations in betaine metabolism and the transsulfuration pathway”, Anal. Bioanal. Chem.,409,6, pp. 1591-1606. Huber, W., V. J. Carey, R. Gentleman, S. Anders, M. Carlson, B. S. Carvalho, H. C. Bravo, S. Davis, L. Gatto, T. Girke, R. Gottardo, F. Hahne, K. D. Hansen, R. A. Irizarry, M. Lawrence, M. I. Love, J. MacDonald, V. Obenchain, A. K. Ole’s, H. Pagès, A. Reyes, P. Shannon, G. K. Smyth, D. Tenenbaum, L. Waldron, and M. Morgan 2015 “Orchestrating high-throughput genomic analysis with Bioconductor”, Nature Methods,12,2, pp. 115-121. Kaelin, William G and Craig B Thompson 2010 “Q&A: Cancer: clues from cell metabolism.” Nature,465,7298, pp. 562- 564. Kamburov, Atanas, Rachel Cavill, Timothy MD Ebbels, Ralf Herwig, and Hector C Keun 2011 “Integrated pathway-level analysis of transcriptomics and metabolomics data with IMPaLA”, Bioinformatics,27,20, pp. 2917-2918. Kanehisa, Minoru, Susumu Goto, Yoko Sato, Miho Furumichi, and Mao Tanabe 2011 “KEGG for integration and interpretation of large-scale molecular data sets”, Nucleic Acids Res.,40, D1, pp. D109-D114. Kessler, Nikolas, Heiko Neuweger, Anja Bonte, Georg Langenkämper, Karsten Niehaus, Tim W Nattkemper, and Alexander Goesmann 2013 “MeltDB 2.0–advances of the metabolomics software system”, Bioinformatics,29,19, pp. 2452-2459. Khatri, Purvesh, Marina Sirota, and Atul J. Butte 2012 “Ten Years of Pathway Analysis: Current Approaches and Outstanding Challenges”, PLoS Comput. Biol.,8,2. Kirkegaard, Thomas and Marja Jäättelä 2009 “Lysosomal involvement in cell death and cancer”, Biochim. Biophys. Acta-Mol. Cell Res.,1793,4, pp. 746-754. Lehtonen, Heli J., Ignacio Blanco, Jose M. Piulats, Riitta Herva, Virpi Launonen, and Lauri A. Aaltonen 2007 “Conventional renal cancer in a patient with fumarate hydratase mutation”, Hum. Pathol.,38,5, pp. 793-796. Maceyka, Michael and Sarah Spiegel 2014 “Sphingolipid metabolites in inflammatory disease”, Nature,510, 7503, p. 58.
118 the r package fella Madsen, Rasmus, Torbjörn Lundstedt, and Johan Trygg 2010 “Chemometrics in metabolomics – a review in human disease diagnosis”, Anal. Chim. Acta,659,1, pp. 23-33. Ni, Ying, Kevin M. Zbuk, Tammy Sadler, Attila Patocs, Glenn Lobo, Emily Edelman, Petra Platzer, Mohammed S. Orloff, Kristin A. Waite, and Charis Eng 2008 “Germline Mutations and Variants in the Succinate Dehydrogenase Genes in Cowden and Cowden-like Syndromes”, Am. J. Hum. Genet., 83,2, pp. 261-268. Page, Lawrence, Sergey Brin, Rajeev Motwani, and Terry Winograd 1999 The PageRank citation ranking: Bringing order to the web. Tech. rep., Stanford InfoLab. Picart-Armada, Sergio, Francesc Fernández-Albert, Maria Vinaixa, Miguel Angel Rodriguez, Suvi Aivio, Travis H. Stracker, Oscar Yanes, and Alexandre Perera-Lluna 2017 “Null diffusion-based enrichment for metabolomics data”, PloS one, 12,12, e0189012. Pithukpakorn, M, M-H Wei, O Toure, P J Steinbach, G M Glenn, B Zbar, W M Linehan, and J R Toro 2006 “Fumarate hydratase enzyme activity in lymphoblastoid cells and fibroblasts of individuals in families with hereditary leiomyomatosis and renal cell cancer.” J. Med. Genet.,43,9, pp. 755-62. Pollard, Patrick, Noel Wortham, and Ian Tomlinson 2003 “The TCA cycle and tumorigenesis: the examples of fumarate hydratase and succinate dehydrogenase”, Ann. Med.,35,8, pp. 634- 639. Seo, Young-Jin, Stephen Alexander, and Bumsuk Hahm 2011 “Does cytokine signaling link sphingolipid metabolism to host defense and immunity against virus infections?”, Cytokine Growth Factor Rev.,22,1, pp. 55-61. Singh, Keshav K, Mohamed M Desouki, Renty B Franklin, and Leslie C Costello 2006 “Mitochondrial aconitase and citrate metabolism in malignant and nonmalignant human prostate tissues.” Mol. cancer,5, p. 14. Vandin, Fabio, Eli Upfal, and Benjamin J Raphael 2011 “Algorithms for detecting significantly mutated pathways in cancer”, J. Comput. Biol.,18,3, pp. 507-522. Xia, Jianguo, Igor V Sinelnikov, Beomsoo Han, and David S Wishart 2015 “MetaboAnalyst 3.0– making metabolomics more meaningful”, Nucleic Acids Res.,43, Web Server issue, W251-W257.
References 119 Yu, Guangchuang, Fei Li, Yide Qin, Xiaochen Bo, Yibo Wu, and Shengqi Wang 2010 “GOSemSim: an R package for measuring semantic similarity among GO terms and gene products”, Bioinformatics,26,7, pp. 976-978. Yu, Guangchuang, Li-Gen Wang, Yanyan Han, and Qing-Yu He 2012 “clusterProfiler: an R package for comparing biological themes among gene clusters”, OMICS,16,5, pp. 284-287. 2014 “Distinct metabolic responses of an ovarian cancer stem cell line”, BMC Syst Biol,8,1, p. 134. Ziarrusta, Haizea, Leire Mijangos, Sergio Picart-Armada, Mireia Irazola, Alexandre Perera-Lluna, Aresatz Usobiaga, Ailette Prieto, Nestor Etxebarria, Maitane Olivares, and Olatz Zuloaga 2018 “Non-targeted metabolomics reveals alterations in liver and plasma of gilt-head bream exposed to oxybenzone”, Chemosphere,211, pp. 624- 631.
8DISEASE GENE IDENTIFICATION benchmarking network propagation methods for disease gene identification In-silico identification of potential target genes for disease is an essential aspect of drug target discovery. Recent studies suggest that successful targets can be found through by leveraging genetic, genomic and protein interaction information. Here, we systematically tested the ability of 12 varied algorithms, based on network propagation, to identify genes that have been targeted by any drug, on gene-disease data from 22 common non-cancerous diseases in OpenTargets. We considered two biological networks, six performance metrics and compared two types of input gene-disease association scores. The impact of the design factors in performance was quantified through additive explanatory models. Standard cross-validation led to over-optimistic performance estimates due to the presence of protein complexes. In order to obtain realistic estimates, we introduced two novel protein complex-aware crossvalidation schemes. When seeding biological networks with known drug targets, machine learning and diffusion-based methods found around 2-4 true targets within the top 20 suggestions. Seeding the networks with genes associated to disease by genetics decreased performance below 1true hit on average. The use of a larger network, although noisier, improved overall performance. We conclude that diffusion-based prioritisers and machine learning applied to diffusion-based features are suited for drug discovery in practice and improve over simpler neighbour-voting methods. We also demonstrate the large impact of choosing an adequate validation strategy and the definition of seed disease genes. 8.1 author summary The use of biological network data has proven its effectiveness in many areas from computational biology. Networks consist of nodes, usually genes or proteins, and edges that connect pairs of nodes, representing information such as physical interactions, regulatory roles or co-occurrence. In order to find new candidate nodes for a given biological property, the so-called network propagation algorithms start from the set of known nodes with that This chapter is a postprint of the following journal article: Picart-Armada, Sergio, Steven J. Barrett, David R. Willé, Alexandre Perera-Lluna, Alex Gutteridge, and Benoit H. Dessailly. “Benchmarking network propagation methods for disease gene identification”. PLoS computational biology 15, no. 9(2019): e1007276. 121
122 disease gene identification property and leverage the connections from the biological network to make predictions. Here, we assess the performance of several network propagation algorithms to find sensible gene targets for 22 common non-cancerous diseases, i.e. those that have been found promising enough to start the clinical trials with any compound. We focus on obtaining performance metrics that reflect a practical scenario in drug development where only a small set of genes can be essayed. We found that the presence of protein complexes biased the performance estimates, leading to over-optimistic conclusions, and introduced two novel strategies to address it. Our results support that network propagation is still a viable approach to find drug targets, but that special care needs to be put on the validation strategy. Algorithms benefitted from the use of a larger -although noisiernetwork and of direct evidence data, rather than indirect genetic associations to disease. 8.2 introduction The pharmaceutical industry faces considerable challenges in the efficiency of commercial drug research and development (Scannell et al.,2012) and in particular in improving its ability to identify future successful drug targets. It has been suggested that using genetic association information is one of the best ways to identify such drug targets (Nelson et al.,2015). In recent years, a large number of highly powered GWAS studies have been published for numerous common traits (for example, (Schizophrenia Working Group of the Psychiatric Genomics Consortium,2014;Verstockt et al.,2018)) and have yielded many candidate genes. Further potential targets can be identified by adding contextual data to the genetic associations, such as genes involved in similar biological processes (Boyle et al.,2017;Jia and Zhao, 2013). Biological networks and biological pathways can be used as a source of contextual data. Biological networks are widely used in bioinformatics and can be constructed from multiple data sources, ranging from macromolecular interaction data collected from the literature (Orchard et al.,2012) to correlation of expression in transcriptomics or proteomics samples of interest (Langfelder and Horvath,2008). A large number of interaction network resources have been made available over the years, many of which are now in the public domain, combining thousands of interactions in a single location (Razick et al.,2008;Szklarczyk, Morris, et al.,2016). They are based on three different fundamental types of data: (1) data-driven networks such as those built by WGCNA (Langfelder and Horvath,2008) for co-expression; (2) interactions extracted from the literature using a human curation process as exemplified by IntAct (Kerrien et al.,2011) or BioGRID (Chatr-Aryamontri et al., 2017); and (3) interactions extracted from the literature using text mining approaches (Al-Aamri et al.,2017). On the other hand, a plethora of network analysis algorithms are available for extracting useful information from such large biological networks in a variety of contexts. Algorithms range in complexity from simple firstneighbour approaches, where the direct neighbours of a gene of interest are assumed to be implicated in similar processes (Piovesan et al.,2015), to ma-
8.3 results 123 chine learning (ML) algorithms designed to learn from the features of the network to make more useful biological predictions (Re et al.,2012). One broad family of network analysis algorithms are the so-called Network Propagation approaches (Cowen et al.,2017), used in contexts such as protein function prediction (Sharan et al.,2007), disease gene identification (Cowen et al.,2017) and cancer gene mutation identification (Leiserson et al.,2014). In this paper, we perform a systematic review of the usefulness of network analysis methods for the purpose of identification of disease genes. As further explained in Methods, we define our test set of disease genes as genes for which the relationship with a disease was sufficiently clear to justify the start of a drug development programme. Claims that such methods are helpful in that context have been made on numerous occasions but a comprehensive validation study is lacking. One major challenge in doing such a study is to define a list of true disease genes for this purpose. To address this, the Open Targets collaboration between pharmaceutical companies and public institutions collects information on known drug targets to help identify new ones (Koscielny et al.,2016). A dedicated internet platform provides a free-to-use accessible resource summarising known data on gene-disease relationships from a number of data sources, like known released drugs and genetic associations from GWAS (Koscielny et al.,2016). The purpose of this work is to quantify the performance of network propagation methods to prioritise novel drug targets, using various networks and validation schemes, and aiming at a faithful reflection of a realistic drug development scenario. We are not predicting gene targets for specific drugs, but rather sensible genes to target for a specific disease. Data on actual compounds targeting a gene is ignored: as long as the gene has been targeted by one or more compounds reaching the clinical trials, it is considered a sensible drug target. We select a number of network propagation approaches that are representative of several classes of algorithms, and test their ability to recover known target genes for several non-cancerous diseases by crossvalidation. We benchmark multiple definitions of disease genes as input for the prioritisers, computational methods, biological networks, validation schemes and performance metrics. We account for all possible combinations of such factors and derive guidelines for future disease target identification studies. The code and data that support our conclusions can be found in https: //github.com/b2slab/genedise. 8.3 results 8.3.1 Benchmark framework Our general approach, summarised in Fig 28, consisted in using a biological network and a list of genes with prior disease-association scores as input to a network propagation approach. We tested some variations of classical network propagation -ppr,raw,gm,mc and zwhich differ on the directedness of the propagation, the input weights and the presence of a statistical normalisation of the scores. Semi-supervised methods included,
130 disease gene identification Figure 32: Ranking of all the methods. Ranking according to the predictions of the main explanatory models (left) and the reduced explanatory models within the STRING network and block cross-validation (right), in both cases on the drugs input and averaging over diseases. The main models serve as a global description of the metrics, whereas the reduced models are specific to the scenario of most interest. A column-wise z-score on the predicted mean is depicted, in order to illustrate the magnitude of the difference. Note how the top 20 hits and the AUPRC metrics lead to similar conclusions, as opposed to AUROC. ilar results; for instance, raw and z. Besides, methods EGAD (arguably one of the simplest) and COSNet (arguably one of the most complex) seemed to result in similar predictions. Fully supervised and semi-supervised approaches largely group in the top right hand quadrant of the STRING plot away from diffusion methods, possibly showing better learning capability with the larger network. Figure 33: Multi-view MDS plot displaying the preserved Spearman’s footrule distances between methods. The differential ranking of their top 100 novel predictions using known drug target inputs are taken into account across all 22 diseases. Results are shown separately for the 2networks considered in this study. Seed genes are excluded from the distance calculations.
8.3 results 131 When comparing overall performances shown in Fig 32 with the prediction differences from the MDS plot (Fig 33), the best methods owed their performance to different reasons as they do not occur within the same region of the plot (e.g. rf and raw). MDS plots on the eight possible combinations of network, input type and inclusion of seed genes are displayed in Figures O and P in S1Appendix. Focusing only on the STRING network and the block validation scheme, we fitted six additive explanatory models, called the reduced models, to model the six metrics for the drugs data input as a function of the method and the disease (see Table G in S1Appendix). Methods were prioritised according to their main effects (Fig 32). The reduced models better described this particular scenario, as they were not forced to fit the trends in all networks and validation schemes in an additive way. Considering the top 20 hits, rf and svm were the optimal choices, followed by wsld and knn. Comparing diseases The top 20 hits model in Fig 29 shows that allergy (the figure’s baseline reference), ulcerative colitis and rheumatoid arthritis (group I) are the diseases for which prediction of target genes was worst, whereas cardiac arrhythmia, Parkinson’s disease, stroke and multiple sclerosis (group II) are those for which it was best. As shown in Fig 34, group I diseases had fewer known target genes and lower modularity compared to group II diseases. Prediction methods worked better when more known target genes were available as input in the network, with two possible underlying reasons: the greater data availability to train the methods, and the natural bias of top 20 hits towards datasets with more positives. Likewise, a stronger modularity within target genes justifies the guilt-by-association principle and led to better performances. In turn, the number of genes and the modularity were positively correlated, see Figure N in S1Appendix. 8.3.3 Performance using genetic associations as input Using genetically associated genes as input to a prediction approach to find known drug targets mimicked a realistic scenario where novel genetic associations are screened as potential targets. However, inferring known drug targets through the indirect genetic evidence posed problems to prediction strategies, especially those based on machine learning. Learning is done using one class of genes in order to predict genes that belong to another class, and the learning space suffers from intrinsic uncertainties in the genetic associations to disease. Both classes are inherently different: certain genes can be difficult to target, and a gene does not require to have been formally associated genetically to a disease to become a valid target. Consequently, we observed a major performance drop on all the prioritisation methods: using any network and cross-validation scheme, the predicted top 20 hits were practically bounded by 1. This was more pronounced on supervised machine learning-focused strategies, as rf and svm lost their edge on diffusion-based strategies. The fact that the genetic associations of the
132 disease gene identification allergy Alzheimers disease arthritis asthma bipolar disorder cardiac arrhythmia COPD coronary heart disease drug dependence hypertension obesity Parkinson's disease psoriasis rheumatoid arthritis schizophrenia stroke lupus multiple sclerosis type I diabetes type II diabetes ulcerative colitis unipolar depression 0 5 10 15 20 0 5 10 15 20 25 Ranking by overall performance (drugs input, top 20 hits). Lower is better. Ranking by modularity (STRING network). 100 150 Genes Lowest rankings (left and down) correspond to largest magnitudes Disease performance in terms of modularity and gene list size Group II Group I More modular Less modular Figure 34: Disease performance in terms of input size and modularity. Disease performance ranked by the number of known target genes and their modularity (obtained using the igraph package, see Figure F in S1Ap- pendix). Modularity is a measure of the tendency of known target genes to form modules or clusters in the network. Diseases have been ranked using their explanatory model coefficient from the top 20 hits metric with known drug targets as input (x axis) and their modularity (y axis). As discussed in the text, best predicted diseases tend to have longer gene lists and be highly modular. validation fold were hidden further hindered the predictions and can be a cause of our pessimistic performance estimates. Comparing cross-validation schemes For reference, we also ran all three cross-validation schemes on the genetic data to quantify and account for complex-related bias. The models confirm that, contrary to the drugs-related input, the differences between the results for the different cross-validation schemes were rather modest. For example, method raw with the STRING network attains 0.59-0.64,0.50-0.54 and 0.37-0.40 hits in the top 20 under the classical, block and representative cross-validation strategies. The slightly larger negative effect on top 20 hits observed with the representative scheme is expected because the number of positives that act as validation decreased and this metric is biased by the class imbalance. The agreement between method ranking using AUPRC and top 20 hits was less consistent, possibly due to the performance drop, whilst AUROC yielded a noticeably different ranking again. Further data can be found in Tables O and P in S1Appendix.
8.4 discussion 133 Comparing networks The change in performance for using the OmniPath network instead of the filtered STRING network was also limited. For AUROC the effect was negative, whereas for the top 20 hits metric the performance improved. Method raw changed from 0.50-0.54 top 20 hits in STRING to 0.61-0.66 in OmniPath under the block validation strategy. Comparing methods To be consistent with the drugs section, we take as reference the block cross-validation strategy and the STRING network. The baseline approach pr that effectively makes use of the network topology alone proved difficult to improve upon, with 0.43-0.47 expected true hits in the top 20. Methods raw and rf respectively achieved 0.50-0.54 and 0.23- 0.26 – although significant, the difference in practice would be minimal. The best performing method was mc with 0.65-0.7hits. All the performance estimates can be found in Table P in S1Appendix. To give an idea of the effort that would be required in a realistic setting to find novel targets, the number of correct hits in the top 100 hits was 3.29-3.45 with the best performing method (in this case, ppr), against 2.25-2.38 of pr. Two main conclusions can be drawn from these results. First, the network topology baseline retained some predictive power upon which most diffusion-based methods, as well as machine-learning approaches COSNet and bagsvm, only managed to add minor improvements, if any. Second, drug targets could still be found by combining network analysis and genes with genetic associations to disease, but with a substantially lower performance and with a marginal gain compared to a baseline approach that would only use the network topology to find targets (e.g. by screening the most connected genes in the network). It is worth noting that gene-disease genetic association scores themselves have drawbacks and that better prediction accuracy could result as genetic association data improves. 8.4 discussion We performed an extensive analysis of the ability of several approaches based on network propagation to identify novel non-cancerous disease target genes. We explored the effect of various choices in factors including the biological network, the definition of disease genes acting as seeds, and the statistical framework being used to evaluate methods performance. We show that carefully choosing an appropriate cross-validation framework and suitable performance metric has an important effect in evaluating the utility of these methods. Our main conclusion is that network propagation seems effective for drug target discovery, reflecting the fact that drug targets tend to cluster within the network. This may be due to the fact that the scientific community has so
134 disease gene identification far been focusing on testing the same proven mechanisms, which can induce some ascertainment bias In a strict cross-validation setting, we found that even the most basic guiltby-association method was useful, with ∼2correct hits in its top 20 predictions, compared to ∼0.1when using a random ranking. The best diffusion based algorithm improved that figure to ∼3.75, and the best overall performing method was a random forest classifier on network-based features (∼4.4 hits). Leading approaches can be notably different in terms of their top predictions, suggesting potential complementarity. We found a better performance when using a network with more coverage at the expense of more false positive interactions. In a more conservative network, random forest performance dropped to ∼3.1hits. Comparing performance on different diseases shows that the more known target genes, and the more clustered these are in the network, the better the performance of network propagation approaches for finding novel targets for it. We also explored the prediction of known drug target genes by seeding the network with an indirect data stream, in particular, genetic association data. Here, the best performing methods were diffusion-based and presented a statistically significant, but marginal, improvement over approaches that only look at network centrality. We conclude that network propagation methods can help identify novel targets for disease, but that the choice of the input network and the seed scores on the genes needs careful consideration. Based on our approach and endorsed benchmarks, we recommend the use of methods employing representations of diffusion-based information (the MashUp network-based features and the diffusion kernels), namely random forest, the support vector machine variants, and raw diffusion algorithms for optimal results. 8.5 materials and methods 8.5.1 Selection of methods for investigation Network propagation algorithms were selected for validation based on the following criteria: 1. Published in a peer-reviewed journal, with evidence of improved performance in gene disease prediction relative to contenders. 2. Implemented as a well documented open source package, that is efficient, robust and usable within a batch testing framework. 3. Directly applicable for gene disease identification from a single gene or protein interaction network, without requiring fundamental changes to the approach or additional annotation information. 4. Capable of outputting a ranked list of individual genes (as opposed to gene modules, for example). In addition, we selected methods that were representative of a diverse panel of algorithms, including diffusion variants, supervised learning on
8.5 materials and methods 135 features derived from network propagation, and a number of baseline approaches (see Table 12). 8.5.2 Testing framework, algorithms and parameterisation All tests and batch runs were set-up and conducted using the R statistical programming language (R Core Team,2016). When no R package was available, the methodology was re-implemented, building upon existing R packages whenever possible. Standard R machine learning libraries were used to train the support vector machine and random forest classifiers. Only the MashUp algorithm (Cho et al.,2016) required feature generation outside of the R environment, using the Matlab code from their publication. Further details on the methods implementation can be found in S1Appendix, section “Method details”. EGAD (Ballouz et al.,2017), a pure neighbour-voting approach, was used here as a baseline comparator. Diffusion (propagation) methods are central in this study. We used the random walk-based personalised PageRank (Page et al.,1999), previously used in similar tasks (Jiang et al.,2017), as implemented in igraph (Csardi and Nepusz,2006). The remaining diffusion-based methods were run on top of the regularised Laplacian kernel (Smola and Kondor,2003), computed through diffuStats (Picart-Armada, Thompson, et al.,2017). We included the classical diffusion raw, a weighted approach version gm that assigns a bias term to the unlabelled nodes, and two statistically normalised scores (mc and z), as implemented in diffuStats. The normalised scores adjust for systematic biases in the diffusion scores that relate to the graph topology, in order to provide a more uniform ranking. In the scope of positive-unlabelled learning (Elkan and Noto,2008;Yang et al.,2012), we included the kernelised scores knn and the linear decayed wsld from RANKS (Valentini, Paccanaro, et al.,2014). knn computes each gene score based on the k-nearest positive examples, using the graph kernel to compute the distances. Conversely, wsld uses all the kernel similarities to the positive examples, but applies a decaying factor to downweight the furthest positives. Closing this category, we implemented the bagging Support Vector Machine approach from ProDiGe1(Mordelet and Vert,2011), here bagsvm, which trains directly on the graph kernel to find the optimal hyperplane separating positive and negative genes. Purer ML-based methods were also included. On one hand, networkbased features were generated using MashUp (Cho et al.,2016) and two classical classifiers were fitted to them, based on caret (Kunn,2008) and mlr (Bischl et al.,2016). These are svm, the Support Vector Machine as implemented in kernlab (Karatzoglou et al.,2004), and rf, the Random Forest found in the randomForest package (Liaw and Wiener,2002). On the other hand, we tried the parametric Hopfield recurrent neural network classifier in the COSNet R package (Bertoni et al.,2011;Frasca et al.,2013). COSNet estimates network parameters on the sub-network containing the labelled nodes, extends them to the sub-network containing the unlabelled ones and then predicts the labels.
136 disease gene identification Table 12: List of methods included in this benchmark. Method identifiers are shortened method names used throughout the text. Other columns are self-explanatory. Method Identifier Method Name Method Class Implementation Reference pr PageRank with a uniform prior Baseline igraph (Bioconductor (Gentleman et al.,2004;Huber et al.,2015) package) (Page et al.,1999) random Random Baseline R (see text) randomraw Random Raw Baseline R (see text) EGAD Extending Guilt by Association’ by Degree Baseline EGAD (Bioconductor package) (Ballouz et al.,2017) ppr Personalized PageRank Diffusion igraph (R package) (Jiang et al.,2017) raw Raw Diffusion Diffusion diffuStats (Bioconductor package) (Vandin et al.,2011) gm GeneMania-based weights Diffusion diffuStats (Bioconductor package) (Mostafavi et al.,2008) mc Monte Carlo normalised scores Diffusion diffuStats (Bioconductor package) (Picart-Armada, Fernández-Albert, et al.,2017) zZ-scores Diffusion diffuStats (Bioconductor package) (Picart-Armada, Fernández-Albert, et al.,2017) knn Knearest neighbours Semi-supervised learning RANKS (R package) (Valentini, Armano, et al.,2016) wsld Weighted Sum with Linear Decay Semi-supervised learning RANKS (R package) (Valentini, Armano, et al.,2016) COSNet COst Sensitive neural Network Supervised learning COSNet (R package) (Frasca et al.,2013) bagsvm Bagging SVM (based on ProDiGe1)Supervised learning kernlab (R package) (Mordelet and Vert,2011) rf Random Forest Supervised learning randomForest (R package) + Matlab (features) (Cho et al.,2016) svm Support Vector Machine Supervised learning kernlab (R package) + Matlab (features) (Cho et al.,2016)
8.5 materials and methods 137 Finally, we defined three naive baseline methods: (1)pr, a PageRank with a uniform prior, where input scores on the genes are ignored; (2) randomraw, which applies the raw diffusion approach to randomly permuted input scores on the genes; and (3)random, a uniform re-ranking of input genes without any network propagation. The inclusion of pr and randomraw allowed us to quantify the predictive power of the network topology alone, without any consideration for the input scores on the genes. 8.5.3 Biological networks The biological network used in the validation is of critical importance as current network resources contain both false positive and false negative interactions, possibly affecting subsequent predictions (Huang et al.,2018). Here, we used two human networks with different general properties, one more likely to contain false positive interactions (STRING (Szklarczyk, Franceschini, et al.,2014)), and another more conservative (OmniPath (Türei et al.,2016)), to test the effect of the network itself on network propagation performance. We further filtered STRING (Szklarczyk, Franceschini, et al., 2014) to retain only a subset of interactions. Having tested several filters, we settled upon high-confidence interactions (combined score > 700) with some evidence from the “Experiments” or “Databases” data sources (see Table B in S1Appendix). Applying these filters and taking the largest connected component resulted in a connected network of 11,748 nodes and 236,963 edges. Edges were assigned weights between 0and 1by rescaling the STRING combined score. We did not filter the OmniPath network (Türei et al.,2016). After removing duplicated edges and taking the largest connected component, the OmniPath network contained 8,580 nodes and 42,145 unweighted edges. 8.5.4 Disease gene data We used the Open Targets platform (Koscielny et al.,2016) to select known disease-related genes. In this analysis we defined positive genes as those reported in Open Targets as being the target of any known drug against the disease of interest, from which all the metrics were computed. We decided to use drug targets, including unsuccessful ones, as proxies for disease genes on the basis that genes for which a drug programme has been started, generally with significant investment, are most likely to have strong evidence linking them to the disease. We therefore regard them as a set of high-confidence true positive disease genes. This choice means we potentially miss genes that have strong genetic associations to the disease but are not druggable. In other words, we focus on limiting false positives in our reference set of positives, at the expense of having more false negatives in our set of negatives. Alternatively, genes with a genetic association of sufficient confidence with the disease were also used as an input data stream, in order to assess the predictive power of an indirect source of evidence. Associations were binarised: any non-zero drugs-related association was considered positive, implying that the methods would predict genes on which a drug has been essayed,
138 disease gene identification regardless of whether the drug was eventually approved. Likewise, only genetic associations with an Open Targets score above 0.16 (see Figure A in S1Appendix) were considered positive. We considered exclusively common diseases with at least 1,000 Open Targets associations, of which a minimum of 50 could be based on known drugs and 50 on genetic associations, in order to avoid empty folds in the nested cross-validations. By applying these filters, we generated a list of phenotypes and diseases which we then manually curated to remove, non-disease phenotype terms (e.g. “body weight and measures”) as well as vague or broad terms (e.g.“cerebrovascular disorder” or “head disease”) and infectious diseases. We also decided to exclude cancers from this analysis. Cancer is a complex process starting from the driver mutation(s) causing disruptive processes involving clonal expansions, which are known to carry their own specific and resultant (non-causal) passenger mutations. Also, the fundamental genetic and biological mechanisms underlying cancers (Hanahan and Weinberg,2011) are generally very distinct from other diseases. We considered this might affect the reliability of the seed genes and cancers would therefore deserve a separate benchmark. This left 22 diseases considered in this study (Table 13). Further descriptive material on the role of genes associated with disease within the STRING network can be found in the section “Descriptive disease statistics in the STRING network” from S1Appendix. Table 13: List of diseases included in this study. Disease N(genetic) N(drugs) Overlap P-value FDR allergy 112 57 1 4.22e-01 4.42e-01 Alzheimers disease 208 103 4 1.10e-01 1.42e-01 arthritis 174 188 6 6.08e-02 1.03e-01 asthma 105 80 6 7.77e-05 5.70e-04 bipolar disorder 117 148 3 1.83e-01 2.12e-01 cardiac arrhythmia 75 177 6 9.15e-04 3.36e-03 chronic obstructive pulmonary disease (COPD) 154 116 6 4.18e-03 1.31e-02 coronary heart disease 111 171 4 7.86e-02 1.24e-01 drug dependence 75 143 6 2.96e-04 1.30e-03 hypertension 66 188 2 2.85e-01 3.14e-01 multiple sclerosis 71 167 4 1.83e-02 4.03e-02 obesity 69 194 3 1.06e-01 1.42e-01 Parkinson’s disease 55 145 0 1 1 psoriasis 131 105 7 1.68e-04 9.23e-04 rheumatoid arthritis 138 95 5 5.18e-03 1.42e-02 schizophrenia 410 163 17 5.44e-05 5.70e-04 stroke 90 156 3 1.18e-01 1.44e-01 systemic lupus erythematosus (lupus) 126 109 5 6.30e-03 1.54e-02 type I diabetes mellitus 87 106 3 4.39e-02 8.04e-02 type II diabetes mellitus 130 154 4 9.14e-02 1.34e-01 ulcerative colitis 136 51 7 1.81e-06 3.98e-05 unipolar depression 123 121 4 3.81e-02 7.63e-02 Diseases included in this study, with a minimum of 50 associated genes both in the known drug targets and the genetic categories (see text). The overlap between these two lists of genes showed a degree of dependence between these two Open Targets data streams for some of the diseases. P-values were calculated using Fisher’s exact test and are reported without and with correction for false discovery rate (Benjamini and Hochberg,1995).
8.5 materials and methods 139 8.5.5 Validation strategies Input gene scores We used the binarised drug association scores and genetic association scores from Open Targets as input gene-level scores to seed the network propagation analyses (Fig 35) and test their ability to recover known drug targets. With the first approach (panel (A) in Fig 35), we tested the predictive power of current network propagation methods for drug target identification using a direct source of evidence (known drug targets). In the second approach (panel (B) in Fig 35), we assessed the ability of a reasonable but indirect source of evidence – genetic associations to disease – in combination with network propagation to recover known drug targets. Figure 35: Input gene scores. Two input types were used to feed the prioritisation algorithms: the binary drug scores in panel (A) and the binary genetic scores in panel (B). In both cases, the validation genes were deemed unlabelled in the input to the prioritisers. Cross-validation folds were always calculated taking into account the drugs input and reused on the genetic input. Metrics Methods were systematically compared using standard performance metrics. The Area under the Receiver Operating Characteristic curve (AUROC) is extensively used in the literature for binary classification of disease genes (Lee et al.,2011), but can be misleading in this context given the extent of the class imbalance between target and non-target genes (Saito and Rehmsmeier, 2015). We however included it in our benchmark for comparison with previous literature. More suitable measures of success in this case are Area under the Precision-Recall curve (AUPRC) (Saito and Rehmsmeier,2015) and partial AUROC (pAUROC) (McClish,1989).
242 the r package diffustats xlab("Protein class")+ ylab("Diffusion score")+ ggtitle("Target proteins scores in ’transport and sensing’") ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●● ● ● ● ●● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ●●● ● ● ● ● raw ml gm ber_s ber_p mc z Other Transport and sensing Other Transport and sensing Other Transport and sensing Other Transport and sensing Other Transport and sensing Other Transport and sensing Other Transport and sensing 0.0 2.5 5.0 7.5 0.00 0.25 0.50 0.75 1.00 0.0 0.2 0.4 0.6 0.00 0.05 0.10 0.15 0.20 −0.6 −0.4 −0.2 −0.6 −0.4 −0.2 0.0 0.00 0.05 0.10 0.15 0.20 Protein class Diffusion score method raw ml gm ber_s ber_p mc z Target proteins scores in 'transport and sensing' As expected, all the diffusion scores qualitatively show differences between positive and negative labels, but the quality of class separation will generally depend on the dataset and scoring method. b.4.4 Benchmarking scores with multiple protein functions The package diffuStats is able to perform several screenings at once. To show its usefulness, we will generalise the procedure in the last section but screening all the categories in the yeast graph. First of all, the input data must meet an adequate format - a straightforward approach is to populate a matrix with the input labels (one category per column). # All classes except NA and unlabelled names_classes <- setdiff(names(table_classes), c("U",NA)) # matrix format mat_classes <- sapply( names_classes, function(class) { V(yeast)$Class %in% class
b.4 getting started 243 } )*1 rownames(mat_classes) <- V(yeast)$name colnames(mat_classes) <- names_classes The former 50% known / 50% unknown approach will be kept with the same split, although not all the 12 categories will be totally balanced in the splits now. All the methods will be compared using the area under the ROC curve (AUROC) as a performance index. Please note that diffuStats is equipped with basic performance measures for rankers: the AUROC, the area under the Precision-Recall curve, or AUPRC, and their partial versions. These are available through the helper function metric_fun and can be passed in list format to perf. These measures are based on the precrec R package - further detail can be found in the original manuscript (Saito and Rehmsmeier,2017). list_methods <- c("raw","ml","gm","ber_s","ber_p","mc","z") df_methods <- perf( K = K_rl, scores = mat_classes[known, ], validation = mat_classes[target, ], grid_param = expand.grid( method = list_methods, stringsAsFactors = FALSE), n.perm = 1000 ) ## Using supplied kernel matrix... ## All done ## Using supplied kernel matrix... ## All done ## Using supplied kernel matrix... ## All done ## Using supplied kernel matrix... ## All done ## Using supplied kernel matrix... ## Using supplied kernel matrix... ## X1: permuting scores... ## Permuting... ## X1: computing heatRank... ## All done ## Using supplied kernel matrix... ## X1: permuting scores... ## Permuting... ## X1: computing heatRank... ## All done ## Using supplied kernel matrix... ## All done
244 the r package diffustats This allows plotting of the AUCs over the categories for each method in one step: ggplot(df_methods, aes(x = method, y = auc)) + geom_boxplot(aes(fill = method)) + scale_fill_npg() + theme_bw() + xlab("Method")+ ylab("Area under the curve")+ ggtitle("Methods performance in all categories") ●● ● 0.6 0.7 0.8 0.9 raw ml gm ber_s ber_p mc z Method Area under the curve method raw ml gm ber_s ber_p mc z Methods performance in all categories Scaling up the analysis can be useful for assessing how adequate a diffusion score is in the dataset of interest. These results suggest that, for the current yeast interactome and protein functions, the best priorisations are those obtained through a statistical normalisation, which might motivate its usage in other biological networks. The user can also statistically compare the performance metrics through the function perf_wilcox. This generates a table with (i) the estimates on the differences on performance between the methods in rows and columns, with their confidence intervals and (ii) their associated p-value (Wilcoxon test). Positive and negative estimates respectively favour the method in the row and the column. # Format the data df_perf <- reshape2::acast(df_methods, Column~method, value.var = "auc") # Compute the comparison matrix
b.5 conclusions 245 df_test <- perf_wilcox( df_perf, digits_p=1, adjust = function(p) p.adjust(p, method = "fdr"), scientific = FALSE) ## Warning in wilcox.test.default(x = perf_mat[, met1], y = perf_mat[, met2], : cannot compute exact p-value with zeroes ## Warning in wilcox.test.default(x = perf_mat[, met1], y = perf_mat[, met2], : cannot compute exact confidence interval with zeroes knitr::kable(df_test, format = "latex") raw ml gm ber_s ber_p mc z raw NA 0.018(-0.075,0.068) -0.0084(-0.082,0.025) NA -0.023(-0.036,-0.0097) -0.022(-0.044,-0.0001) -0.027(-0.051,-0.0068) ml 0.55 NA -0.025(-0.047,-0.013) -0.018(-0.068,0.075) -0.042(-0.095,0.047) -0.047(-0.098,0.043) -0.044(-0.1,0.026) gm 0.71 0.01 NA 0.0084(-0.025,0.082) -0.02(-0.049,0.054) -0.02(-0.051,0.057) -0.017(-0.054,0.038) ber_s NA 0.55 0.71 NA -0.023(-0.036,-0.0097) -0.022(-0.044,-0.0001) -0.027(-0.051,-0.0068) ber_p 0.03 0.46 0.54 0.03 NA 0.0015(-0.0076,0.0071) -0.0055(-0.018,0.0079) mc 0.12 0.35 0.55 0.12 0.73 NA -0.0047(-0.014,0.0025) z0.03 0.32 0.46 0.03 0.46 0.34 NA b.5 conclusions The diffuStats package is a new computational tool to compute and compare single-network diffusion scores that are object of active research in several bioinformatics areas. It is an effort to gather a collection of settings in the diffusion process like the graph kernel, the label codification and the choice of a statistical normalisation. The diffuStats package will help the end user in choosing and computing the best performing diffusion scores in the application of interest. b.6 funding This work was supported by the Spanish Ministry of Economy and Competitiveness (MINECO) [TEC2014-60337-R to A.P.] and the National Institutes of Health (NIH) [R01GM104400 to W.T.]. AP. and S.P. thank for funding the Spanish Biomedical Research Centre in Diabetes and Associated Metabolic Disorders (CIBERDEM) and the Networking Biomedical Research Centre in the subject area of Bioengineering, Biomaterials and Nanomedicine (CIBER-BBN), both initiatives of Instituto de Investigación Carlos III (ISCIII). SP. thanks the AGAUR FI-scholarship programme.
246 the r package diffustats b.7 session info Here is the output of sessionInfo() on the system that compiled this vignette: •R version 3.6.2(2019-12-12), x86_64-pc-linux-gnu •Locale: LC_CTYPE=en_US.UTF-8,LC_NUMERIC=C,LC_TIME=es_ES.UTF-8, LC_COLLATE=en_US.UTF-8,LC_MONETARY=es_ES.UTF-8, LC_MESSAGES=en_US.UTF-8,LC_PAPER=es_ES.UTF-8,LC_NAME=C, LC_ADDRESS=C,LC_TELEPHONE=C,LC_MEASUREMENT=es_ES.UTF-8, LC_IDENTIFICATION=C •Running under: Ubuntu 16.04.6 LTS •Matrix products: default •BLAS: /usr/lib/atlas-base/atlas/libblas.so.3.0 •LAPACK: /usr/lib/atlas-base/atlas/liblapack.so.3.0 •Base packages: base, datasets, graphics, grDevices, methods, stats, utils •Other packages: diffuStats 1.4.0, ggplot2 3.1.1, ggsci 2.9, igraph 1.2.4.1, igraphdata 1.0.1, knitr 1.22 •Loaded via a namespace (and not attached): assertthat 0.2.1, backports 1.1.4, BiocManager 1.30.4, BiocStyle 2.12.0, colorspace 1.4-1, compiler 3.6.2, crayon 1.3.4, data.table 1.12.2, digest 0.6.18, dplyr 0.8.3, evaluate 0.13, expm 0.999-4, glue 1.3.1, grid 3.6.2, gtable 0.3.0, highr 0.8, htmltools 0.3.6, labeling 0.3, lattice 0.20-38, lazyeval 0.2.2, magrittr 1.5, MASS 7.3-51.5, Matrix 1.2-18, munsell 0.5.0, pillar 1.4.0, pkgconfig 2.0.2, plyr 1.8.4, precrec 0.10.1, purrr 0.3.2, R6 2.4.0, Rcpp 1.0.1, RcppArmadillo 0.9.800.3.0, RcppParallel 4.4.4, reshape2 1.4.3, rlang 0.4.0, rmarkdown 1.12, scales 1.0.0, stringi 1.4.3, stringr 1.4.0, tcltk 3.6.2, tibble 2.1.1, tidyselect 0.2.5, tools 3.6.2, vctrs 0.2.0, withr 2.1.2, xfun 0.6, yaml 2.2.0, zeallot 0.1.0
References 247 references Allaire, JJ, Romain Francois, Kevin Ushey, Gregory Vandenbrouck, Marcus Geelnard, and Intel 2016 RcppParallel: Parallel Programming Tools for ’Rcpp’, R package version 4.3.20,https://CRAN.R-project.org/package=RcppParallel. Bersanelli, Matteo, Ettore Mosca, Daniel Remondini, Gastone Castellani, and Luciano Milanesi 2016 “Network diffusion-based analysis of high-throughput data for the detection of differentially enriched modules.” Scientific Reports,6, August, p. 34841,issn:2045-2322,doi:10.1038/srep34841. Csardi, Gabor 2015 “igraphdata: A Collection of Network Data Sets for the ’igraph’ Package”, R package version 1.0.1,https://CRAN.R-project.org/ package=igraphdata. Csardi, Gabor and Tamas Nepusz 2006 “The igraph software package for complex network research”, InterJournal, Complex Systems, p. 1695,http://igraph.org. Eddelbuettel, Dirk 2013 Seamless R and C++ integration with Rcpp, Springer. Eddelbuettel, Dirk and Conrad Sanderson 2014 “RcppArmadillo: Accelerating R with high-performance C++ linear algebra”, Computational Statistics and Data Analysis,71 (Mar. 2014), pp. 1054-1063,doi:10.1016/j.csda.2013.02.005. Harchaoui, Zaid, Francis Bach, Olivier Cappe, and Eric Moulines 2013 “Kernel-based methods for hypothesis testing: A unified view”, IEEE Signal Processing Magazine,30,4, pp. 87-97. Lee, Insuk, U Martin Blom, Peggy I Wang, Jung Eun Shim, and Edward M Marcotte 2011 “Prioritizing candidate disease genes by network-based boosting of genome-wide association data”, Genome Research,21,7, pp. 1109- 1121,doi:10.1101/gr.118992.110. Mewes, Hans-Werner, Dmitrij Frishman, Christian Gruber, Birgitta Geier, Dirk Haase, Andreas Kaps, Kai Lemcke, Gertrud Mannhaupt, Friedhelm Pfeiffer, C Schüller, et al. 2000 “MIPS: a database for genomes and protein sequences”, Nucleic acids research,28,1, pp. 37-40. Mostafavi, Sara, Debajyoti Ray, David Warde-Farley, Chris Grouios, and Quaid Morris 2008 “GeneMANIA: a real-time multiple association network integration algorithm for predicting gene function.” Genome Biology,9Suppl 1, S4,issn:1474-760X, doi:10.1186/gb-2008-9-s1-s4.
248 the r package diffustats North, Bernard V, David Curtis, and Pak C Sham 2002 “A note on the calculation of empirical P values from Monte Carlo procedures”, The American Journal of Human Genetics,71,2, pp. 439- 441. Oliver, Stephen 2000 “Guilt-by-association goes global”, Nature,403, February, pp. 601- 603. Paull, Evan O., Daniel E. Carlin, Mario Niepel, Peter K. Sorger, David Haussler, and Joshua M. Stuart 2013 “Discovering causal pathways linking genomic events to transcriptional states using Tied Diffusion Through Interacting Events (TieDIE)”, Bioinformatics,29,21, pp. 2757-2764,issn:13674803,doi:10.1093/ bioinformatics/btt471. R Core Team 2017 R: A Language and Environment for Statistical Computing, R Foundation for Statistical Computing, Vienna, Austria, https :// www.R - project.org/. Saito, Takaya and Marc Rehmsmeier 2017 “Precrec: fast and accurate precision–recall and ROC curve calculations in R”, Bioinformatics,33,1, pp. 145-147. Smola, Alexander J and Risi Kondor 2003 “Kernels and regularization on graphs”, pp. 144-158,doi:10.1007/ 978-3-540-45167-9_12. Suthram, Silpa, Andreas Beyer, Richard M Karp, Yonina Eldar, and Trey Ideker 2008 “eQED: an efficient method for interpreting eQTL associations using protein networks”, Molecular systems biology,4,1, p. 162. Tsuda, Koji, HyunJung J. Shin, and Bernhard Schölkopf 2005 “Fast protein classification with multiple networks”, Bioinformatics, 21, SUPPL. 2, pp. 59-65,issn:13674803,doi:10.1093/bioinformat ics/bti1110. Valentini, Giorgio, Alberto Paccanaro, Horacio Caniza, Alfonso E. Romero, and Matteo Re 2014 “An extensive analysis of disease-gene associations using network integration and fast kernel-based gene prioritization methods”, Artificial Intelligence in Medicine,61,2, pp. 63-78,issn:18732860,doi: 10.1016/j.artmed.2014.03.003. Vandin, Fabio, Eli Upfal, and Benjamin J. Raphael 2010 “Algorithms for detecting significantly mutated pathways in cancer”, Lect. Notes Comput. Sci.,6044 LNBI, 3, pp. 506-521,issn:03029743, doi:10.1007/978-3-642-12683-3_33.
References 249 Wickham, Hadley 2009 ggplot2: Elegant Graphics for Data Analysis, Springer-Verlag New York, isbn:978-0-387-98140-6,http://ggplot2.org. Yen, Luh, Francois Fouss, and Christine Decaestecker 2007 “Graph nodes clustering based on the commute-time kernel”, Advances in Knowledge Discovery and Data Mining, pp. 1037-1045,issn: 03029743,doi:10.1007/978-3-540-71701-0_117. Zoidi, Olga, Eftychia Fotiadou, Nikos Nikolaidis, and Ioannis Pitas 2015 “Graph-Based Label Propagation in Digital Media: A Review”, ACM Computing Surveys,47,3,48:1-48:35,issn:0360-0300,doi:10.1145/ 2700381.
CMETABOLOMICS ENRICHMENT c.1 appendix s1 - graph structure and curation KEGG database Raw KEGG graphKEGG in list format Raw KEGG graph (largest CC) Raw KEGG weighted graph (largest CC) 11 1 1 1 1 2 2 22 2 KEGG graph (curated, inverted weights) 11 1 1 1 1 2 1/ Figure 73: Procedure to obtain the KEGG graph. The KEGG database is read as a collection of lists that contain the annotations. The raw KEGG graph is built through these annotations, where the vertices are KEGG entries from categories in Fig. 74a. We only work with the largest CC of the raw KEGG graph, to which weights are assigned, enabling the curation step that gives place to the KEGG graph. Note that the weights in the definitive KEGG graph are the inverse of the former dissimilarity weights, to be consistent with the diffusive methods. Compounds Reactions Enzymes Modules Pathways Genes (a) (b) Figure 74: (a) Structure of our KEGG graph. Each entry belongs to a level, ranging from 1(pathways) to 5(compounds). (b) MetScape concept of compound-reaction-enzyme-gene network. Their construction is similar in the three lowest levels, while in the upper level they include KEGG genes. The first step to depict current knowledge is to build a graph object from KEGG that enables data enrichment (Fig. 73). The KEGG graph contains various categories (Fig. 74a) and keeps similarities with the networks built through MetScape (Karnovsky et al.,2012), see (Fig. 74b), although our structure is conceived to include biological pathways and modules to obtain This appendix reproduces the supplementary data (Appendices S1to S5) of: Picart-Armada, Sergio, Francesc Fernández-Albert, Maria Vinaixa, Miguel A. Rodríguez, Suvi Aivio, Travis H. Stracker, Oscar Yanes, and Alexandre Perera-Lluna. “Null diffusion-based enrichment for metabolomics data”. PloS one 12, no. 12 (2017). 251
258 metabolomics enrichment where Id is the identity matrix, dis the damping factor and Mis a matrix obtained from the weighted adjacency n×nmatrix Afrom the directed KEGG graph (edges pointing upwards): M=f(t(A)) (56) In this expression, tis the matrix transpose operator and fis the function that normalises each column to sum 1, except when it contains nzeroes; in the latter case it returns a column with nelements equal to 1 n. We apply the function fwith that particularity to be coherent with the R package igraph (Csardi and Nepusz,2006), which considers that terminal nodes resume the random walk uniformly distributed in all the vertices. In Equation 54, the calculation of the PageRank scores uses the binary vector of affected compounds normalised by the amount of affected compounds (Fig. 78): p=G PGi (57) The similarity between the final expression for heat diffusion (see Appendix S2) and PageRank (Equation 54) is remarkable, given the common random walk background for these two methods. The differences between them include the forced upwards directionality of PageRank and the damping factor concept, which allows leaps in the diffusion. The rescaling of the p vector does not affect the null model, as the rescaling factor remains constant in the random trials. All the PageRank scores in our approach have been computed using the standard d=0.85 established in the original publication (Page et al.,1999), a range of damping factors has been swept, going from 0.1(very frequent restarts) to 0.95 (almost no restarts), but results are consistent as a result of the application of our null model, described in Appendix S4. c.4 appendix s4 - null models We retrieve the formulation from Appendix S2to compute the final temperatures in heat diffusion. The same procedure applies to the PageRank approach in Appendix S3, given the similarity between them, but as a proof of concept it will be developed for heat diffusion only. T=RHD ·G By abuse of notation, RHD will contain the columns of the original RHD corresponding only to compounds, thus being a rectangular matrix from now on. Likewise, to have a well-defined matrix-vector product, Gwill only refer to compounds, as they are the only entities that can introduce heat. The vector Gcontains exactly nin ones, corresponding to the affected compounds, and ncomp −nin zeroes, where ncomp is the amount of compounds in the graph. When focusing on a node i, we want to assess whether its temperature Tiwould be expected from a random selection of affected compounds or,
c.4 appendix s4 - null models 259 on the contrary, the affected compounds are more related to node ithan expected. To that end, we define the null distribution of temperatures: Tnull =RHD ·X(58) where Xis the random variable obtained by permuting G. If we define p=nin ncomp , then every Xiis a Bernoulli trial with success probability p.Xi and Xj, for i6=j, are slightly anticorrelated due to the permutation approach. Exact statistical moments of Tnull can be computed: E(Tnull) = RHD ·E(X)(59) being E(X) = p· 1 1 . . . 1 (60) and the same for the covariance matrix Σ(Tnull) = RHD ·Σ(X)·RT HD (61) where Σ(X) = p(1−p)· 1 ρ . . . ρ ρ 1 . . . ρ . . .. . ..... . . ρ ρ . . . 1 (62) being ρ= − 1 ncomp−1. In terms of the elements rij in matrix RHD, we can write µi=nin ncomp ncomp X j=1 rij (63) σ2 i=nin(ncomp −nin) ncomp(ncomp −1) ncomp X j=1 r2 ij −1 ncomp ncomp X j=1 rij 2 (64) In fact, the correlation matrix of the null temperatures gives insights about the combination of the network structure and the null model. Focusing on pathways: if two pathways are strongly correlated, it suggests that both attain high or low temperatures with similar inputs. Thus, this couple of pathways are prone to overlap or to be nearby in the metabolism. Conversely, a strong anticorrelation suggests that warming up a pathway conditions the second pathway to become colder, suggesting that these pathways are dissimilar. (Fig. 79) depicts the correlations matrix for the pathways, where rows and columns have been reordered to illustrate clusters of KEGG pathways. Furthermore, the structure seems to be somehow related to the underlying biology: one of the clusters corresponds to the bulk of human
260 metabolomics enrichment metabolic pathways, whereas the genetic information processing pathways also appear highly correlated. Besides these examples, pathway types do not appear totally shuffled, but as small blocks of pathways sharing a biological role. Cellular Processes Environmental Information Processing Genetic Information Processing Human Diseases Metabolism Organismal Systems +10-1 Correlation Pathway type Figure 79: Correlation matrix between pathways for the heat diffusion process in the KEGG graph (identifiers omitted for clarity). Blue correlations are closer to 1, while red tend to −1. For the calculation of these correlations, a count of 33 compounds was assumed to perform the null model, like in the experimental data. In addition, the biological role of the pathways annotated in KEGG BRITE (Kanehisa et al.,2008) has been drawn in the rightmost bar. (Eqs. 63,64) provide a first approximation (1)to evaluate how high a temperature is. If Tiis remarkably greater than µi=E(Tnull)iin terms of its standard deviation, σi=pΣ(Tnull)ii, then node ishould be reported. The normalised score is zi=Ti−µi σi (65) Another approach (2)to evaluate the relevance of node iis to perform the permutation analysis through Monte Carlo trials. In that case, the random vector Xis drawn nperm times and, for node i, the p-value is approximated as pi=ri+1 nperm+1, where riis the number of trials where Tnulli⩾Ti, see (North et al.,2002) for further details on this estimator. An ensemble solution can be obtained by repeating the procedure nvote times and evaluating each node by majority vote, specifically including it only if it is reported at least bnvote 2c+1times, also allowing a fuzzy representation of the consensus solution. This ensemble approach reduces the variability in the reported solution and also provides confidence measures for each included node. The input can also be subsampled in this approach, although this option has not been explored yet.
c.5 appendix s5 - details on reported solutions 261 c.5 appendix s5 - details on reported solutions The solutions reported in the main body encompass two different scoring functions (heat diffusion and PageRank) and two statistical approaches (z-scores and simulation). Monte Carlo permutations involve consensus solutions among the simulated runs to reduce the variability and enhance robustness and consistency between runs. c.5.1 Solution stratification (Fig. 80) depicts the stratification of all the reported graphs. Solutions tend to keep the same proportions as the original graph, allowing the discovery of relevant nodes in all the categories. This is not only a sign of agreement across the different solutions, but also a necessary behaviour to discover putative nodes in all the categories. The proportion of reported compounds seems lower than the one in the KEGG graph, probably due to the application of the null model, which tends to favour the metabolites in the input and penalise the rest. 0.00 0.25 0.50 0.75 1.00 KEGG graph HD norm HD sim PR norm PR sim Method Proportion of nodes Pathway Module Enzyme Reaction Compound Figure 80: Subgraph stratification by method. Notice the tendency to keep the same distribution as the whole KEGG graph. Compared to the KEGG graph, the proportion of compounds decreases in all the solutions due to the inclination to recover the ones in the input and exclude the rest.
262 metabolomics enrichment c.5.2 Connected component evolution The choice of the number of desired nodes kaffects the number and size of the reported connected components (CC). In general terms, the number of reported CCs seems to lower as the number of reported nodes grows (Fig. 81), as different CCs that contain seed nodes tend to merge. The number of nodes in the largest CC grows monotonically as kincreases and it captures the majority of the nodes in the subgraphs (Fig. 82). 5 10 15 0 100 200 300 Top k nodes Number of CCs diffusion pagerank normality simulation Figure 81: Number of reported CCs for varying k. The discrete nature of the majority vote approach in simulation trials leads to ties, which can be spotted as horizontal lines in the figure - more than a hundred nodes are tied with the best rank in both cases. In general, the number of CCs decreases as the reported subgraphs grow.
c.5 appendix s5 - details on reported solutions 263 0 50 100 150 200 250 0 100 200 300 Top k nodes Nodes in largest CC diffusion pagerank normality simulation Figure 82: Number of nodes in the largest CC reported. If more nodes are reported, the largest CC grows accordingly. For k=300 around 225-250 nodes are in the largest CC in every approach, meaning that approximately 75-83% nodes lie in it. As more than a hundred nodes are tied with maximum rank in the simulated version, lines start at these respective points.
264 metabolomics enrichment c.5.3 Computational cost In sight of further applications of these methodologies, we have performed a benchmark of several implementations of the Monte Carlo approach. For heat diffusion, the temperature calculation (Appendix S2) is achieved through T= −KI−1·G=RHD ·G(66) where Ghas nin ones and ncomp −nin zeroes. We define two strategies to permute the input G:(a) draw nin elements without replacement from the set [1,ncomp], and (b) shuffle the whole vector G. Furthermore, two possible calculations for Tare: (1)explicitly compute RHD and sum the columns indexed by the nin ones, or (2)solve the linear system T= −KI−1·G These strategies have been essayed with growing graph order, from 1,000 to 10,000 nodes. Graphs are randomly generated using the default Barabási- Albert model in igraph (Csardi and Nepusz,2006). Then, 10% of the nodes are randomly selected to be pathways, so the heat flow can be dispelled. We also consider two scenarios, depending on if (I) the input list has a fixed size of 30 compounds, or (II) it scales as 10% of the graph nodes. For each combination of parameters, a benchmark of 30 permutations is run 10 times and the trends are depicted in (Fig. 83). The differences between sampling strategies seem irrelevant compared to the solving method. Solving the linear system seems a good option for small inputs that do not scale with the graph order, but computing the inverse seems to scale better as the vector G becomes less sparse. However, the latter requires a vast amount of memory to store the matrix, so the best implementation will depend on the graph order and memory availability. All the benchmarks have been executed in a desktop workstation (Intel i5 650 at 3.2GHz, 16Gb RAM memory).
c.5 appendix s5 - details on reported solutions 265 Fixed input size Variable input size ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ●●● ● ● ● ● ● ● ● ● ●●●●●●●● ● ● ●●●● ● ● ● ● ● ●● ● ●●●● ● ● ● ● ● ●●● ●● ● ●● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ●● ●● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ●● ● ●● ●●● ● ● ● ● ● ●● ● ●● ● ● ●● ●●● ● ● ●● ● ●● ● ● ●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ●● ● ● ● ● ●● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●●●● ● ● ● ● ● ● ● ● ● ●● ●● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● 10−1.5 10−1 10−0.5 100 100.5 1000 2000 3000 4000 5000 6000 7000 8000 9000 10000 1000 2000 3000 4000 5000 6000 7000 8000 9000 10000 Number of nodes Elapsed time (s) Method ●Inverse Solve Sampling ● ● All Partial Figure 83: Computational cost of several strategies for computing 30 permutations, that is, 30 null temperatures for all the nodes. These simulations have been performed 10 times with each combination of parameters. In the left figure, the input size is kept constant and equal to 30 compounds, whereas in the right it scales with the graph nodes (10%). Two methods to compute the temperatures are compared: direct resolution (solve) and explicit computation of the RHD matrix (inverse). Likewise, two sampling strategies are explored: permute the whole Gvector (all) or just draw the nin compounds (partial).
266 metabolomics enrichment c.5.4 Damping factor influence The damping factor in the PageRank setup (see Appendix S3) is a model parameter that could affect the results if set differently. Although the standard value d=0.85 was used, we analysed the parameter sensibility by sweeping several values of d, computing the z-scores and reporting the top 250 nodes (Figs. 84,85). The normalised scores seem stable in a wide range of choices of d, therefore its choice does not seem a critical issue.
c.5 appendix s5 - details on reported solutions 267 ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ●● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ●● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ●● ●● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● d = 0.1 ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● d = 0.3 ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ●● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ●● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ●● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● d = 0.5 ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ●● ● ● ● ● ● ●● ● ● ●● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● d = 0.55 ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ●● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ●● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● d = 0.6 ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ●● ●● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ●● ● ●● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ●● ● ●● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● d = 0.65 Legend Pathway Module Enzyme Reaction Compound (input: ) Figure 84: Damping factor impact. The normalised z-scores show consistent solutions for a range of damping factors.
274 the r package fella getExcluded(myAnalysis) ## [1] "intruder1" "intruder2" "intruder3" "intruder4" "intruder5" ## [6] "intruder6" "intruder7" "intruder8" "intruder9" "intruder10" Keep in mind that exact matching is sought, so be extremely careful with whitespaces, tabs or similar characters that might create mismatches. For example: input.fail <- paste0(" ", input.full) defineCompounds( compounds = input.fail, data = FELLA.sample) ## Error in defineCompounds(compounds = input.fail, data = FELLA.sample): None of the specified compounds appear in the available KEGG data. d.1.4 Enriching the data Once the FELLA.DATA and the FELLA.USER with the affected metabolites are ready, the data can be easily enriched. Enrichment methods There are three methods to enrich: 1.Hypergeometric test (method = "hypergeom"): it performs the metabolitesampling hypergeometric test using the connections in FELLA’s KEGG graph. This is included for completeness and does not include the contextual novelty of the diffusive methods. 2.Diffusion (method = "diffusion"): it performs sub-network analysis on the KEGG graph to extract a meaningful subgraph. This subgraph can be plotted an interpreted 3.PageRank (method = "pagerank"): analogous to "diffusion" but using the directed diffusion, which matches the PageRank algorithm for web ranking. Statistical approximations For methods "diffusion" and "pagerank", two statistical approximations are proposed: 1.Normal approximation (approx = "normality"): scores are computed through z-scores based on analytical expected value and covariance matrix of the null model for diffusion. This approximation is deterministic and fast. 2.Monte Carlo trials (approx = "simulation"): scores are computed through Monte Carlo trials of the random variables. This approximation requires computing the random trials, governed by the ntrials argument.
d.1 additional file 1: quickstart 275 Enrichment: methods, approximations and wrapper function The function enrich wraps the functions defineCompounds,runHypergeom, runDiffusion and runPagerank in an easily usable manner, returning a FELLA.USER object with complete analyses. myAnalysis <- enrich( compounds = input.full, method = "diffusion", approx = "normality", data = FELLA.sample) ## No background compounds specified. Default background will be used. ## Warning in defineCompounds(compounds = compounds, compoundsBackground = ## compoundsBackground, : Some compounds were introduced as affected but they ## do not belong to the background. These compounds will be excluded from the ## analysis. Use ’getExcluded’ to see them. ## Running diffusion... ## Computing p-scores through the specified distribution. ## Done. The output is quite informative and aggregates all the warnings. Let’s compare an empty FELLA.USER object show(new("FELLA.USER")) ## Compounds in the input: empty ## Background compounds: all available compounds (default) ## ----------------------------- ## Hypergeometric test: not performed ## ----------------------------- ## Heat diffusion: not performed ## ----------------------------- ## PageRank: not performed to the output of a processed one: show(myAnalysis) ## Compounds in the input: 30 ## [1] "C00003" "C00011" "C00027" "C00044" "C00081" "C00083" "C00091" ## [8] "C00111" "C00143" "C00164" "C00288" "C00332" "C00334" "C00479" ## [15] "C00546" "C01352" "C04225" "C04253" "C05258" "C05264" "C05266" ## [22] "C05268" "C05275" "C05280" "C05467" "C14145" "C15979" "C16328" ## [29] "C16329" "C16333" ## Background compounds: all available compounds (default) ## -----------------------------
276 the r package fella ## Hypergeometric test: not performed ## ----------------------------- ## Heat diffusion: ready. ## P-scores under 0.05: 86 ## ----------------------------- ## PageRank: not performed The wrapper function enrich can run the three analysis at once with the option method = listMethods(), or only the desired ones providing them as a character vector: myAnalysis <- enrich( compounds = input.full, method = listMethods(), approx = "normality", data = FELLA.sample) show(myAnalysis) ## Compounds in the input: 30 ## [1] "C00003" "C00011" "C00027" "C00044" "C00081" "C00083" "C00091" ## [8] "C00111" "C00143" "C00164" "C00288" "C00332" "C00334" "C00479" ## [15] "C00546" "C01352" "C04225" "C04253" "C05258" "C05264" "C05266" ## [22] "C05268" "C05275" "C05280" "C05467" "C14145" "C15979" "C16328" ## [29] "C16329" "C16333" ## Background compounds: all available compounds (default) ## ----------------------------- ## Hypergeometric test: ready. ## Top 2 p-values: ## hsa00640 hsa00010 ## 8.540386e-09 9.999888e-01 ## ## ----------------------------- ## Heat diffusion: ready. ## P-scores under 0.05: 86 ## ----------------------------- ## PageRank: ready. ## P-scores under 0.05: 70 The wrapped functions work in a similar way, here is an example with runDiffusion: myAnalysis_bis <- runDiffusion( object = myAnalysis, approx = "normality", data = FELLA.sample) ## Running diffusion... ## Computing p-scores through the specified distribution.
d.1 additional file 1: quickstart 277 ## Done. show(myAnalysis_bis) ## Compounds in the input: 30 ## [1] "C00003" "C00011" "C00027" "C00044" "C00081" "C00083" "C00091" ## [8] "C00111" "C00143" "C00164" "C00288" "C00332" "C00334" "C00479" ## [15] "C00546" "C01352" "C04225" "C04253" "C05258" "C05264" "C05266" ## [22] "C05268" "C05275" "C05280" "C05467" "C14145" "C15979" "C16328" ## [29] "C16329" "C16333" ## Background compounds: all available compounds (default) ## ----------------------------- ## Hypergeometric test: ready. ## Top 2 p-values: ## hsa00640 hsa00010 ## 8.540386e-09 9.999888e-01 ## ## ----------------------------- ## Heat diffusion: ready. ## P-scores under 0.05: 86 ## ----------------------------- ## PageRank: ready. ## P-scores under 0.05: 70 d.1.5 Visualising the results The method plot for data from the package FELLA allows a friendly visualisation of the relevant part of the KEGG graph. Hypergeom In the case method = "hypergeom" the plot encompasses a bipartite graph that contains top pathways and affected compounds. In that case, threshold = 1 allows the visualisation of both pathways; otherwise a plot with only one pathway would be quite uninformative. plot( x=myAnalysis, method = "hypergeom", main = "My first enrichment using the hypergeometric test in FELLA", threshold = 1, data = FELLA.sample)
278 the r package fella My first enrichment using the hypergeometric test in FELLA C00003 C00011 C00027 C00044 C00081 C00083 C00091 C00111 C00143 C00164 C00288 C00332 C00334 C00479 C00546 C01352 C04225 C04253 C05258 C05264 C05266 C05268 C05275 C05280 C05467 C14145 C15979 C16328 C16329 C16333 hsa00640 hsa00010 Diffusion For method = "diffusion" the graph contains a richer representations involving modules, enzymes and reactions that link affected pathways and compounds. plot( x=myAnalysis, method = "diffusion", main = "My first enrichment using the diffusion analysis in FELLA", threshold = 0.1, data = FELLA.sample)
d.1 additional file 1: quickstart 279 My first enrichment using the diffusion analysis in FELLA Propanoate metabolism ... Malonate semialdehyde ... long−chain−3−hydroxyac... 3−hydroxyacyl−CoA dehy... malonate−semialdehyde ... methylmalonate−semiald... acetyl−CoA C−acetyltra... 3−hydroxyisobutyryl−Co... phosphoenolpyruvate ca... (S)−methylmalonyl−CoA ... malonyl−CoA decarboxyl... enoyl−CoA hydratase dodecenoyl−CoA isomera... succinate−−−CoA ligase... succinate−−−CoA ligase... pyruvate:NAD+ 2−oxidor... malonyl−CoA carboxy−ly... Acetyl−CoA:acetyl−CoA ... malonyl−CoA:pyruvate c... Succinate:CoA ligase (... GTP:oxaloacetate carbo... Succinate:CoA ligase (... 3−Oxopropanoate:NAD+ o... 3−Oxopropanoate:NADP+ ... ITP:oxaloacetate carbo... Succinate:CoA ligase (... 3−oxopropanoate:NADP+ ... acetyl−CoA:carbon−diox... (R)−Methylmalonyl−CoA ... propanoyl−CoA:phosphat... (S)−methylmalonyl−CoA ... propanoyl−CoA:electron... propanoyl−CoA:oxaloace... (S)−Methylmalonate sem... 2−oxobutanoate:oxygen ... glycerone−phosphate ph... glycine synthase acetoacetyl−CoA:acetat... Acetoacetate carboxy−l... propanoyl−CoA:carbon−d... (S)−3−Hydroxybutanoyl−... (S)−lactaldehyde:NADP+... propane−1,2−diol hydro... Methylglyoxal + NADPH ... 3−Hydroxypropionyl−CoA... (2S,3S)−2−hydroxybutan... (2S,3R)−3−Hydroxybutan... propanoyl−CoA:electron... (S)−3−Hydroxyhexadecan... (S)−3−Hydroxyhexadecan... (S)−hydroxydecanoyl−Co... (S)−Hydroxydecanoyl−Co... (S)−hydroxyoctanoyl−Co... (S)−Hydroxyoctanoyl−Co... (S)−hydroxyhexanoyl−Co... (S)−Hydroxyhexanoyl−Co... decanoyl−CoA:electron−... cis,cis−3,6−Dodecadien... 3alpha,7alpha,12alpha,... (S)−3−hydroxyacyl−CoA:... (3S)−3−hydroxyacyl−CoA... 3−Hydroxy−OPC8−CoA <=>... 3−Hydroxy−OPC8−CoA + N... 3−Hydroxy−OPC6−CoA <=>... 3−Hydroxy−OPC6−CoA + N... 2−Oxoglutarate dehydro... propanal:NAD+ oxidored... propanoyl−CoA:NAD+ oxi... (S)−Lactaldehyde + NAD... hydroxyacetone:NAD+ 1−... 2−Methyl−trans−aconita... CoA CO2 Acetyl−CoA Hydrogen peroxide GTP ITP Malonyl−CoA Succinyl−CoA Propanoyl−CoA Glycerone phosphate 5,10−Methylenetetrahyd... Acetoacetate HCO3− Acetoacetyl−CoA Propanal Methylglyoxal FADH2 (Z)−But−2−ene−1,2,3−tr... Electron−transferring ... (S)−3−Hydroxyhexadecan... (S)−Hydroxydecanoyl−Co... (S)−3−Hydroxyoctanoyl−... (S)−Hydroxyhexanoyl−Co... trans−Dec−2−enoyl−CoA trans,cis−Lauro−2,6−di... cis,cis−3,6−Dodecadien... 3alpha,7alpha,12alpha−... (3S)−3−Hydroxyadipyl−C... trans−2−Enoyl−OPC8−CoA 3−Hydroxy−OPC8−CoA 3−Hydroxy−OPC6−CoA Categories for each node Pathway Module Enzyme Reaction Compound Input compound PageRank For method = "pagerank" the concept is analogous to diffusion: plot( x=myAnalysis, method = "pagerank", main = "My first enrichment using the PageRank analysis in FELLA", threshold = 0.1, data = FELLA.sample)
280 the r package fella My first enrichment using the PageRank analysis in FELLA Propanoate metabolism ... Malonate semialdehyde ... 3−hydroxyacyl−CoA dehy... phosphoenolpyruvate ca... malonyl−CoA decarboxyl... enoyl−CoA hydratase triose−phosphate isome... dodecenoyl−CoA isomera... methylmalonyl−CoA muta... succinate−−−CoA ligase... succinate−−−CoA ligase... acetyl−CoA carboxylase propionyl−CoA carboxyl... malonyl−CoA carboxy−ly... Acetyl−CoA:acetyl−CoA ... malonyl−CoA:pyruvate c... Succinate:CoA ligase (... GTP:pyruvate 2−O−phosp... GTP:oxaloacetate carbo... Succinate:CoA ligase (... ITP:oxaloacetate carbo... Succinate:CoA ligase (... 3−oxopropanoate:NADP+ ... acetyl−CoA:carbon−diox... (R)−Methylmalonyl−CoA ... 2−oxobutanoate:oxygen ... D−glyceraldehyde−3−pho... glycerone−phosphate ph... beta−D−fructose−1,6−bi... glycine synthase acetoacetyl−CoA:acetat... Acetoacetate carboxy−l... propanoyl−CoA:carbon−d... (S)−3−Hydroxybutanoyl−... 4−Aminobutyraldehyde:N... propane−1,2−diol hydro... 4−aminobutanal:NAD+ 1−... (S)−2−methylbutanoyl−C... (S)−3−Methyl−2−oxopent... (2S,3S)−2−hydroxybutan... (2S,3R)−3−Hydroxybutan... (S)−3−Hydroxyhexadecan... (S)−Hydroxydecanoyl−Co... (S)−Hydroxyoctanoyl−Co... (S)−hydroxyhexanoyl−Co... (S)−Hydroxyhexanoyl−Co... decanoyl−CoA:electron−... cis,cis−3,6−Dodecadien... 3alpha,7alpha,12alpha,... 2−Methyl−1−hydroxybuty... 3−Hydroxy−OPC8−CoA <=>... 2−Oxoglutarate dehydro... propanal:NAD+ oxidored... (S)−Lactaldehyde + NAD... hydroxyacetone:NAD+ 1−... 2−Methyl−trans−aconita... NAD+ CO2 Hydrogen peroxide GTPITP Malonyl−CoA Succinyl−CoA Glycerone phosphate 5,10−Methylenetetrahyd... Acetoacetate HCO3− Acetoacetyl−CoA 4−Aminobutanoate Propanal Methylglyoxal FADH2 (Z)−But−2−ene−1,2,3−tr... (S)−3−Hydroxyhexadecan... (S)−Hydroxydecanoyl−Co... (S)−3−Hydroxyoctanoyl−... (S)−Hydroxyhexanoyl−Co... trans−Dec−2−enoyl−CoA cis,cis−3,6−Dodecadien... 3alpha,7alpha,12alpha−... [Dihydrolipoyllysine−r... trans−2−Enoyl−OPC8−CoA 3−Hydroxy−OPC8−CoA Categories for each node Pathway Module Enzyme Reaction Compound Input compound d.1.6 Exporting the results FELLA offers several exporting alternatives, both for the R environment and for external software. Exporting inside R The appropriate functions to export the results inside R are generateResultsTable for a data.frame object: myTable <- generateResultsTable( object = myAnalysis, method = "diffusion", threshold = 0.1, data = FELLA.sample) ## Writing diffusion results... ## Done. knitr::kable(head(myTable, 20)) .. .and generateResultsGraph for a graph in igraph format:
d.1 additional file 1: quickstart 281 KEGG.id Entry.type KEGG.name p.score hsa00640 pathway Propanoate metabolism - Homo sapiens (human) 0.0036894 M00013 module Malonate semialdehyde pathway, propanoyl-CoA . . . 0.0044683 1.1.1.211 enzyme long-chain-3-hydroxyacyl-CoA dehydrogenase 0.0371099 1.1.1.35 enzyme 3-hydroxyacyl-CoA dehydrogenase 0.0392511 1.2.1.18 enzyme malonate-semialdehyde dehydrogenase (acetylat.. . 0.0069255 1.2.1.27 enzyme methylmalonate-semialdehyde dehydrogenase (Co.. . 0.0165439 2.3.1.9enzyme acetyl-CoA C-acetyltransferase 0.0085923 3.1.2.4enzyme 3-hydroxyisobutyryl-CoA hydrolase 0.0786804 4.1.1.32 enzyme phosphoenolpyruvate carboxykinase (GTP) 0.0700429 4.1.1.41 enzyme (S)-methylmalonyl-CoA decarboxylase 0.0223899 4.1.1.9enzyme malonyl-CoA decarboxylase 0.0002538 4.2.1.17 enzyme enoyl-CoA hydratase 0.0015731 5.3.3.8enzyme dodecenoyl-CoA isomerase 0.0164255 6.2.1.4enzyme succinate—CoA ligase (GDP-forming) 0.0019142 6.2.1.5enzyme succinate—CoA ligase (ADP-forming) 0.0125330 R00209 reaction pyruvate:NAD+ 2-oxidoreductase (CoA-acetylati.. . 0.0885938 R00233 reaction malonyl-CoA carboxy-lyase (acetyl-CoA-forming. . . 0.0000698 R00238 reaction Acetyl-CoA:acetyl-CoA C-acetyltransferase 0.0001037 R00353 reaction malonyl-CoA:pyruvate carboxytransferase 0.0065794 R00405 reaction Succinate:CoA ligase (ADP-forming) 0.0468613 myGraph <- generateResultsGraph( object = myAnalysis, method = "diffusion", threshold = 0.1, data = FELLA.sample) show(myGraph) ## IGRAPH 6bf1c19 UNW- 102 166 -- ## + attr: name (v/c), com (v/n), NAME (v/x), entrez (v/x), label ## | (v/c), input (v/l), weight (e/n) ## + edges from 6bf1c19 (vertex names): ## [1] hsa00640--M00013 M00013 --1.1.1.211 M00013 --1.1.1.35 ## [4] M00013 --1.2.1.18 M00013 --1.2.1.27 hsa00640--2.3.1.9 ## [7] M00013 --3.1.2.4 hsa00640--4.1.1.41 hsa00640--4.1.1.9 ## [10] M00013 --4.2.1.17 M00013 --5.3.3.8 hsa00640--6.2.1.4 ## [13] hsa00640--6.2.1.5 4.1.1.9 --R00233 2.3.1.9 --R00238 ## [16] hsa00640--R00353 6.2.1.5 --R00405 4.1.1.32--R00431 ## [19] 6.2.1.4 --R00432 1.2.1.18--R00705 1.2.1.27--R00705 ## + ... omitted several edges Exporting outside R Results can be saved as permanent files. The data.frame data format can be saved as a .csv file: myTempDir <- tempdir() myExp_csv <- paste0(myTempDir, "/table.csv") exportResults(
282 the r package fella format = "csv", file = myExp_csv, method = "pagerank", threshold = 0.1, object = myAnalysis, data = FELLA.sample) ## Exporting to a csv file... ## Writing pagerank results... ## Done. ## Done test <- read.csv(file = myExp_csv) knitr::kable(head(test)) KEGG.id Entry.type KEGG.name p.score hsa00640 pathway Propanoate metabolism - Homo sapiens (human) 0.0000085 M00013 module Malonate semialdehyde pathway, propanoyl-CoA .. . 0.0010330 1.1.1.35 enzyme 3-hydroxyacyl-CoA dehydrogenase 0.0422528 4.1.1.32 enzyme phosphoenolpyruvate carboxykinase (GTP) 0.0088747 4.1.1.9enzyme malonyl-CoA decarboxylase 0.0005280 4.2.1.17 enzyme enoyl-CoA hydratase 0.0003343 In the same line, the graph can be saved in RData: myExp_graph <- paste0(myTempDir, "/graph.RData") exportResults( format = "igraph", file = myExp_graph, method = "pagerank", threshold = 0.1, object = myAnalysis, data = FELLA.sample) ## Exporting to a RData file using ’igraph’ object... ## Done stopifnot("graph.RData" %in% list.files(myTempDir)) Other formats exported by igraph are also available, internally using their function igraph::write.graph. Check the format argument of ?igraph::write.graph for a list of the supported formats. For example, using "pajek" format:
d.1 additional file 1: quickstart 283 myExp_pajek <- paste0(myTempDir, "/graph.pajek") exportResults( format = "pajek", file = myExp_pajek, method = "diffusion", threshold = 0.1, object = myAnalysis, data = FELLA.sample) ## Exporting to the format pajek using igraph... ## Done stopifnot("graph.pajek" %in% list.files(myTempDir)) This option is toggled if the format does not match any other predefined export option. d.1.7 Session info For reproducibility purposes, below is the sessionInfo() output: sessionInfo() ## R version 3.6.2 (2019-12-12) ## Platform: x86_64-pc-linux-gnu (64-bit) ## Running under: Ubuntu 16.04.6 LTS ## ## Matrix products: default ## BLAS: /usr/lib/atlas-base/atlas/libblas.so.3.0 ## LAPACK: /usr/lib/atlas-base/atlas/liblapack.so.3.0 ## ## locale: ## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C ## [3] LC_TIME=es_ES.UTF-8 LC_COLLATE=en_US.UTF-8 ## [5] LC_MONETARY=es_ES.UTF-8 LC_MESSAGES=en_US.UTF-8 ## [7] LC_PAPER=es_ES.UTF-8 LC_NAME=C ## [9] LC_ADDRESS=C LC_TELEPHONE=C ## [11] LC_MEASUREMENT=es_ES.UTF-8 LC_IDENTIFICATION=C ## ## attached base packages: ## [1] parallel stats4 stats graphics grDevices utils datasets ## [8] methods base ## ## other attached packages: ## [1] magrittr_1.5 igraph_1.2.4.1 KEGGREST_1.24.1 ## [4] org.Mm.eg.db_3.8.2 org.Hs.eg.db_3.8.2 AnnotationDbi_1.46.0 ## [7] IRanges_2.17.5 S4Vectors_0.21.24 Biobase_2.44.0 ## [10] BiocGenerics_0.29.2 FELLA_1.5.3 knitr_1.22
290 the r package fella Figure 87: Internal knowledge representation from KEGG. The scheme outlines the KEGG graph, a heterogeneous network whose nodes belong to a category in KEGG: compound, reaction, enzyme, module or pathway. Lower levels are expected to be more specific entities, while top levels are broader concepts. The enrichment procedure starts from input metabolites and extracts a relevant sub-network from the KEGG graph. Figure extractred from (Picart-Armada et al.,2017) The user should be aware that KEGG is frequently updated and therefore the derived KEGG graph can change between KEGG releases. The metadata from the KEGG version used to build a FELLA.DATA object can be retrieved through getInfo. Enrichment analysis Once the database is ready as a FELLA.DATA object and the input is formatted as a list of KEGG compounds, the enrichment can be performed. The results of the enrichment are stored in a FELLA.USER object, possibly using three methodologies described below. Hypergeometric test For completeness purposes, the hypergeometric test is included in FELLA in the function runHypergeom. As in several ORA implementations, the hypergeometric distribution is used to assess whether a biological pathway contains more hits within the input list than expected from chance given its size. Pathways are ranked according to their p-value after multiple testing correction. Note that the results from this test will differ from a hypergeometric test using the original KEGG pathways, because metabolite-pathway connections are inferred from the KEGG graph. A metabolite is included in a pathway if the pathway can be reached from the metabolite in the upwards-directed KEGG graph, depicted in figure 89. In consequence, metabolites related to the enzymes within a pathway will belong to the pathway, even if they were not in the original definition of the KEGG pathway.
d.2 additional file 2: main vignette 291 Compounds Reactions Enzymes Modules Pathways Figure 88: Network setup for the diffusion process. nput metabolites (in black rings) introduce a unitary flow in the network and only the pathway nodes (blue rings) can leak the flow. The final score of the nodes reflects the “temperature” of a stationary state. Figure extractred from (Picart- Armada et al.,2017). diffusion Diffusion algorithms have been extensively used in computational biology. For instance, HotNet is an algorithm for finding sub-networks with a large amount of mutated genes (Vandin et al.,2011), whereas TieDIE attemps to link a source set and a target set of molecular entities through two diffusion processes (Paull et al.,2013). Other applications include the prioritisation of disease genes (Lee et al.,2011) and the prediction of gene function (Mostafavi et al.,2008). In FELLA, diffusion is a natural way to score all the nodes in the KEGG graph given an input list of metabolites, available using method = "diffusion" in the function runDiffusion. The input metabolites introduce unitary flow in the network. Flow can only leave the network through pathway nodes, forcing it to propagate through the intermediate entities as well (reactions, enzymes and modules), see figure 88. Further details can found in (Picart- Armada et al.,2017). However, the diffusion scores are biased due to the network topology (Picart-Armada et al.,2017) and therefore a normalisation step is required. FELLA offers a normalisation through a z-score (approx = "normality") or through an empirical p-value (approx = "simulation"), both assessing whether the diffusion score of a node is likely to be reached in a permutation analysis, i.e. if the input is random. The normalisation through the z-scores leads to p-scores, defined as: psi=1−Φ(zi) Where psiis the p-score of node i,ziis its z-score (Picart-Armada et al., 2017) and Φis the cumulative distribution function of the standard gaussian distribution. Under this definition, nodes are ranked using increasing pscores. For completeness, two alternative parametric scores have been added. The heavier-tailed t-distribution can be used instead of the gaussian by choosing approx = "t" and supplying the desired degrees of freedom ν.
292 the r package fella Compounds Reactions Enzymes Modules Pathways Figure 89: Network setup for PageRank. Input metabolites (in black rings) are the source of random walks that must climb through the graph levels, up to the pathway nodes. Figure extractred from (Picart-Armada et al.,2017). Similarly, the gamma distribution can be used through approx = "gamma". The p-score is obtained with psi=1−Fi(Ti) Being Tithe raw temperature of node iand Fithe cumulative distribution function of a gamma distribution, adjusted by its shape (µ2 i σ2 i ) and scale (σ2 i µi) parameters. The quantities µiand σ2 iare the mean and variance of the null temperatures and are analytically known from the null model formulation (Picart-Armada et al.,2017). pagerank PageRank (Page et al.,1999) offers a scoring method for the nodes in the KEGG graph, based on a random walks approach. The random walks start at the input metabolites and are forced to explore their reachable nodes, see figure 89. As random walks take into account the direction of the edges, PageRank is applied to the upwards-directed KEGG graph (figure 87) in order to force the walks to reach pathway nodes. Nodes that are frequently visited by the random walks earn a higher PageRank, analogously to the diffusion scores. More details about this particular formulation, implemented in runPagerank, can be found in (Picart-Armada et al.,2017). The PageRank scores are statistically normalised, providing the same options as in the diffusion scores in section D.2.3. Therefore, the argument approx can be set to "simulation" for the permutation analysis, or to "normality", "t" or "gamma" for the parametric alternatives. Enrichment wrapper FELLA contains the wrapper enrich that maps the KEGG ids and runs the desired enrichment procedure with a single call. This can be convenient for producing compact scripts and running quick analyses.
d.2 additional file 2: main vignette 293 Limitations FELLA currently starts the statistical analysis from a list of affected metabolites. Therefore, it inherits a limitation from ORA methods: the need of choosing a cutoff to derive the list of affected metabolites, assuming that the metabolites stem from a differential abundance analysis. Another limitation, shared among network-based models, is the incomplete biological knowledge from which the network is built. The knowledge model in FELLA might also constraint the complexity of the mechanisms that can be found through it. Processes such as genetic and epigenetic events, or the type and directionality of regulatory events, are not considered at the moment. The user should be aware that FELLA neither builds a dynamic model of the biochemical reactions in the metabolism, nor relies on flux balance analysis. Conversely, FELLA is built on a knowledge representation from the biology in KEGG that focuses on offering interpretability to the final user. d.2.4 Case studies The functionalities of FELLA are demonstrated by (1) building a Homo sapiens database and (2) enriching summary metabolomics data from three public datasets. Building the database FELLA requires a database built from KEGG to perform any data enrichment. FELLA contains a small example database as a FELLA.DATA object, accessible via data("FELLA.sample"), but this is a toy example for demonstration purposes, not suited for regular analyses. Therefore, the database for the corresponding organism has to be built before any analysis is run. The first step is to build the KEGG graph from the current KEGG release with the function buildGraphFromKEGGREST. Note that the user can force specific KEGG pathways to be excluded from the graph - the following code removes “overview” metabolic pathways based on KEGG brite. library(FELLA) set.seed(1) # Filter overview pathways graph <- buildGraphFromKEGGREST( organism = "hsa", filter.path = c("01100","01200","01210","01212","01230")) ## Building through KEGGREST... ## Available gene annotations: ncbi-geneid, ncbi-proteinid. Using ncbi-geneid ## Done. ## Building graph...
294 the r package fella ## Filtering 5 pathways. ## Done. ## Pruning graph... ## Current weight: 1 out of 4... ## Current weight: 2 out of 4... ## Current weight: 3 out of 4... ## Current weight: 4 out of 4... ## Done. Once the KEGG graph is ready, the database will be saved locally using buildDataFromGraph. The user can choose which matrices shall be stored using the matrices argument - saving both "diffusion" and "pagerank" might take up to 1GB of disk space. If the user plans on using the z-score approximation, it is advisable to set the normality argument to c("diffusion", "pagerank") in order to speed up future computations. Using the z-scores with a custom metabolite background will require the matrices to be saved as well. Finally, the argument niter controls how many random trials are performed in the estimation of the null distribution of the largest connected component of a k-th order random subgraph. As this is a property of the KEGG graph, it is performed once and reused in each analysis. This finds application when filtering small connected components from the reported sub-network, see section D.2.4. tmpdir <- paste0(tempdir(), "/my_database") # Mke sure the database does not exist from a former vignette build # Otherwise the vignette will rise an error # because FELLA will not overwrite an existing database unlink(tmpdir, recursive = TRUE) buildDataFromGraph( keggdata.graph = graph, databaseDir = tmpdir, internalDir = FALSE, matrices = "diffusion", normality = "diffusion", niter = 50) ## Computing probabilities for random subgraphs... (this may take a while) ## Directory /tmp/RtmpAOhDlw/my_database does not exist. Creating it... ## Done. ## Done. ## Computing diffusion.matrix... (this may take a while and use some memory) ## Done ## Computing diffusion.rowSums... ## Done. When the database is available in local, it can be loaded in an Rsession and assigned to a FELLA.DATA object using the function loadKEGGdata. This
d.2 additional file 2: main vignette 295 should be the only procedure for creating any FELLA.DATA object. The user is given the choice of loading the diffusion and pagerank matrices to ease memory saving. fella.data <- loadKEGGdata( databaseDir = tmpdir, internalDir = FALSE, loadMatrix = "diffusion" ) ## Loading KEGG graph data... ## Done. ## Loading hypergeom data... ## Loading matrix... ## ’hypergeom.matrix.RData’ not present in:/tmp/RtmpAOhDlw/my_database/hypergeom.matrix.RData. Hypergeometric test won’t execute. ## Done. ## Loading diffusion data... ## Loading matrix... ## Done. ## Loading rowSums... ## Done. ## Loading pagerank data... ## Loading matrix... ## ’pagerank.matrix.RData’ not loaded. Simulated permutations may execute slower for pagerank. ## Done. ## Loading rowSums... ## ’pagerank.rowSums.RData’ not present in:/tmp/RtmpAOhDlw/my_database/pagerank.rowSums.RData. Z-scores won’t be available for pagerank. ## Done. ## Data successfully loaded. The contents of the FELLA.DATA object can be summarised as well: fella.data ## General data: ## - KEGG graph: ## *Nodes: 11115 ## *Edges: 34787 ## *Density: 0.0002816029 ## *Categories: ## + pathway [327] ## + module [173] ## + enzyme [1149] ## + reaction [5467] ## + compound [3999] ## *Size: 6.2 Mb ## - KEGG names are ready.
296 the r package fella ## ----------------------------- ## Hypergeometric test: ## - Matrix not loaded. ## ----------------------------- ## Heat diffusion: ## - Matrix is ready ## *Dim: 11115 x 3999 ## *Size: 340.1 Mb ## - RowSums are ready. ## ----------------------------- ## PageRank: ## - Matrix not loaded. ## - RowSums not loaded. The function getInfo provides the KEGG release and organism that generated a FELLA.DATA object: cat(getInfo(fella.data)) ## T01001 Homo sapiens (human) KEGG Genes Database ## hsa Release 93.0+/02-22, Feb 20 ## Kanehisa Laboratories ## 22,498 entries ## ## linked db pathway ## brite ## module ## ko ## genome ## enzyme ## network ## disease ## drug ## ncbi-geneid ## ncbi-proteinid ## uniprot Please note that the database built for this vignette is stored in a temporary folder and will not be persistent. The user should build his or her own database and save it in a persistent location, either in the package installation directory (internalDir = TRUE) or in a custom folder (internalDir = FALSE). Internal databases can be listed using listInternalDatabases. A cautionary note if the user is relying on the internal directory: reinstalling FELLA will wipe existent databases because its internal directory is overwritten. Also, if the database name already exists when saving a new database, the existing database will be renamed by appending _old in order to avoid overwriting.
d.2 additional file 2: main vignette 297 Epithelial cells dataset This example data is extracted from the epithelial cancer cells dataset (Chen et al.,2015), an in vitro model of dry eye in which the human epithelial cells IOBA-NHC are put under hyperosmotic stress. The original study files are deposited in the Metabolights repository (Haug et al.,2012) under the identifier MTBLS214:https://www.ebi.ac.uk/metabolights/MTBLS214. The list of metabolites hereby used reflects metabolic changes in “Treatment 1” (24 hours in serum-free media at 380 mOsm) against control (24 hours at 280 mOsm). The metabolites have been extracted from “Table 1” in the original manuscript and mapped to KEGG ids. mapping the input metabolites The input metabolites should be provided as KEGG compound identifiers. If the user starts from another source (common names, HMDB identifiers), tools like the “compound ID conversor” from MetaboAnalyst can be useful for the ID conversion. compounds.epithelial <- c( "C02862","C00487","C00025","C00064", "C00670","C00073","C00588","C00082","C00043") The first step is to map the input metabolites to the KEGG graph with defineCompounds. This step requires the FELLA.DATA object, loaded in section D.2.4. The user can impose a custom metabolite background with the compoundsBackground argument. By default, all the KEGG compounds in the graph are used. analysis.epithelial <- defineCompounds( compounds = compounds.epithelial, data = fella.data) ## No background compounds specified. Default background will be used. ## Warning in defineCompounds(compounds = compounds.epithelial, data = fella.data): Some compounds were introduced as affected but they do not belong to the background. These compounds will be excluded from the analysis. Use ’getExcluded’ to see them. Notice that defineCompounds throws a warning if any of the input metabolites does not map to the graph. The user can retrieve the mapped and unmapped identifiers through getInput and getExcluded, respectively. getInput(analysis.epithelial) ## [1] "C00025" "C00043" "C00064" "C00073" "C00082" "C00487" "C00588" "C00670" getExcluded(analysis.epithelial) ## [1] "C02862" The status of a FELLA.USER object can be checked by printing the object.
298 the r package fella analysis.epithelial ## Compounds in the input: 8 ## [1] "C00025" "C00043" "C00064" "C00073" "C00082" "C00487" "C00588" "C00670" ## Background compounds: all available compounds (default) ## ----------------------------- ## Hypergeometric test: not performed ## ----------------------------- ## Heat diffusion: not performed ## ----------------------------- ## PageRank: not performed enriching using diffusion Having mapped the compounds, the enrichment can be performed. In this vignette, only the diffusion method in runDiffusion will be applied, although PageRank has an almost identical usage in runPagerank. If the user prefers an explicit permutation analysis, the option approx = "simulation" performs the amount of iterations specified in the niter argument. Conversely, if the desired approximation is the z-score (approx = "normality"), the process does not require permutations. The z-scores are converted to p.scores using the pnorm routine. Likewise, approx = "t" and approx = "gamma" respectively rely on pt and pgamma. Section D.2.3contains further details on the scores. This example applies approx = "normality", a fast option. For a comparison between prioritisations using Monte Carlo trials or the parametric z-score, the user can is referred to (Picart-Armada et al.,2017). analysis.epithelial <- runDiffusion( object = analysis.epithelial, data = fella.data, approx = "normality") ## Running diffusion... ## Computing p-scores through the specified distribution. ## Done. The FELLA.USER object has been updated with the p.scores from the diffusion results: analysis.epithelial ## Compounds in the input: 8 ## [1] "C00025" "C00043" "C00064" "C00073" "C00082" "C00487" "C00588" "C00670" ## Background compounds: all available compounds (default) ## ----------------------------- ## Hypergeometric test: not performed ## ----------------------------- ## Heat diffusion: ready.
d.2 additional file 2: main vignette 299 ## P-scores under 0.05: 282 ## ----------------------------- ## PageRank: not performed At this point, the subgraph consisting of top scoring nodes can be plotted in a heterogeneous network layout. In the presence of signal, this subgraph will exhibit large connected components and contain nodes from all the levels in the KEGG graph. It is also expected that the algorithm gives a high priority to the metabolites specified in the input, although not all of them must necessarily be top ranked. Therefore, the user should expect to find the presence of intermediate entities (reactions, enzymes and modules) that connect the input to relevant KEGG pathways. Note that FELLA can also pinpoint new KEGG compounds as potentially relevant. In this example, the plot is limited to 150 nodes using the nlimit argument from plot. nlimit <- 150 vertex.label.cex <- .5 plot( analysis.epithelial, method = "diffusion", data = fella.data, nlimit = nlimit, vertex.label.cex = vertex.label.cex) ## 282 nodes below the threshold have been limited to 150 nodes.
306 the r package fella Figure 90: Graphical interface: compounds upload Figure 91: Graphical interface: advanced options helper functions FELLA is equipped with helper functions that ease the user experience and avoid direct manipulation of the S4classes. Some of them have been already introduced - a complete enumeration of the exported functions is hereby provided. Functions of the type getease object and slot retrieval, with the following possibilities: getBackground,getExcluded,getInfo,getInput,getName, getPscores.
d.2 additional file 2: main vignette 307 Figure 92: Graphical interface: results Figure 93: Graphical interface: export On the other hand, functions starting by listprovide general purpose data about the package (listMethods,listApprox,listCategories) and a listing of the available internal databases (listInternalDatabases). Finally, functions starting by ischeck if an object belongs to a certain class: is.FELLA.DATA and is.FELLA.USER. Ovarian cancer cells dataset The next example has been extracted from the study on metabolic responses of ovarian cancer cells (Vermeersch et al.,2014). The original files
308 the r package fella can be found in the MTBLS150 study in the Metabolights respository: https: //www.ebi.ac.uk/metabolights/MTBLS150. OCSCs are isogenic ovarian cancer stem cells derived from the OVCAR-3ovarian cancer cells. The abundances of six metabolites are affected by the exposure to several environmental conditions: glucose deprivation, hypoxia and ischemia (column “All” in “Figure 3” from their main manuscript). The common names have been converted to KEGG ids prior to applying FELLA. The analysis is performed using the wrapper enrich that maps the compounds to the internal representation and runs the desired methods. compounds.ovarian <- c( "C00275","C00158","C00042", "C00346","C00122","C06468") analysis.ovarian <- enrich( compounds = compounds.ovarian, data = fella.data, methods = "diffusion") ## No background compounds specified. Default background will be used. ## Warning in defineCompounds(compounds = compounds, compoundsBackground = compoundsBackground, : Some compounds were introduced as affected but they do not belong to the background. These compounds will be excluded from the analysis. Use ’getExcluded’ to see them. ## Running diffusion... ## Computing p-scores through the specified distribution. ## Done. plot( analysis.ovarian, method = "diffusion", data = fella.data, nlimit = 150, vertex.label.cex = vertex.label.cex, plotLegend = FALSE) ## 176 nodes below the threshold have been limited to 150 nodes.
d.2 additional file 2: main vignette 309 ● ● ● ●● ● ● ●● ● ● ● ● ● ●● ● ● ● ● ● ● ●● ● ● ● ● ● ● ●● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ●● ●● ● ● ● ● ● ●● ● ● ●● ● ● ● ● ● ●● ● ● ● ● ●● ● ● ● ● ● ● ●●● ●● ● ● ● ● ● ●● ● ●● ● ● ● ● ● ●● ●● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ●● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● Citrate cycle (TCA cyc... Alanine, aspartate and... Sphingolipid signaling... Lysosome − Homo sapien... Apoptosis − Homo sapie... Antigen processing and... Fc gamma R−mediated ph... Taste transduction − H... Glucagon signaling pat... Type II diabetes melli... Renal cell carcinoma −... Central carbon metabol... Choline metabolism in ... Citrate cycle (TCA cyc... Citrate cycle, first c... Citrate cycle, second ... Glyoxylate cycle Urea cycle Adenine ribonucleotide... Phosphatidylethanolami... Sphingosine degradatio... Ascorbate biosynthesis... Reductive citrate cycl... Dicarboxylate−hydroxyb... Methylaspartate cycle Arginine biosynthesis,... malate dehydrogenase isocitrate dehydrogena... isocitrate dehydrogena... peptide−aspartate beta... phytanoyl−CoA dioxygen... succinate−semialdehyde... succinate dehydrogenas... citrate (Si)−synthase ATP citrate synthase 3−deoxy−D−glycero−D−ga... hexokinase phosphorylase kinase [pyruvate dehydrogenas... choline kinase ethanolamine kinase sphingosine kinase mannose−1−phosphate gu... ethanolamine−phosphate... deoxyribonuclease II fructose−2,6−bisphosph... phosphoethanolamine/ph... mannosyl−oligosacchari... endoplasmic reticulum ... alpha−mannosidase dipeptidyl−peptidase I tripeptidyl−peptidase ... carboxypeptidase C cathepsin X beta−aspartyl−peptidas... cathepsin B cathepsin L cathepsin H cathepsin S legumain cathepsin K cathepsin F cathepsin O cathepsin V caspase−2 caspase−6 caspase−10 cathepsin E cathepsin D N4−(beta−N−acetylgluco... fumarylacetoacetase acylpyruvate hydrolase sphinganine−1−phosphat... fumarate hydratase aconitate hydratase ethanolamine−phosphate... argininosuccinate lyas... adenylosuccinate lyase mannose−6−phosphate is... glucose−6−phosphate is... phosphomannomutase succinate−−−CoA ligase... beta−citrylglutamate s... pyruvate carboxylase acetyl−CoA:oxaloacetat... acetyl−CoA:oxaloacetat... succinate:NAD+ oxidore... Succinate:CoA ligase (... succinyl−CoA:(S)−malat... succinyl−CoA:acetoacet... Succinate:CoA ligase (... isocitrate glyoxylate−... L−aspartate ammonia−ly... succinate−semialdehyde... succinate−semialdehyde... Succinate:CoA ligase (... ethanolamine−phosphate... ATP:D−fructose 6−phosp... D−mannose aldose−ketos... GDP−mannose mannophosp... GDP:D−mannose−1−phosph... GTP:alpha−D−mannose−1−... O−Succinyl−L−homoserin... (S)−malate hydro−lyase... N6−(1,2−dicarboxyethyl... 3−fumarylpyruvate fuma... 2−(Nomega−L−arginino)s... Maleate cis−trans−isom... L−Proline,2−oxoglutara... Citrate:CoA ligase (AD... citrate hydroxymutase citrate hydro−lyase (c... ATP:D−mannose 6−phosph... 4−fumarylacetoacetate ... ATP:ethanolamine O−pho... D−mannose 6−phosphate ... D−mannose−6−phosphate ... (S)−dihydroorotate:fum... S−Adenosyl−L−methionin... CTP:ethanolamine−phosp... Phosphatidylethanolami... succinate:quinone oxid... Sphinganine−1−phosphat... protein−N(pi)−phosphoh... serine−phosphoethanola... 1−(5'−Phosphoribosyl)−... Benzylsuccinate fumara... sphingosine−1−phosphat... alpha 1,2−mannosylolig... (2Z,4E,7E)−2−Hydroxy−6... Phosphoethanolamine ph... Lysoplasmalogen ethano... Plasmenylethanolamine ... citrate:N6−acetyl−N6−h... succinyl−CoA:acetate C... fumarate CoM:CoB oxido... citrate:L−glutamate li... propionyl−CoA:succinat... L−2,3−diaminopropanoat... D−Ornithine + Citrate ... N5−Citryl−D−ornithine ... G10694 + 3 H2O <=> G00... Succinate Fumarate Citrate D−Mannose D−Mannose 6−phosphate Ethanolamine phosphate D−Mannose 1−phosphate beta−D−Fructose 6−phos... The resulting subnetwork2reports several TCA cycle-related entities, also reported by the authors and by previous work (Pollard et al.,2003). It also mentions sphingosine degradation, closely related to the reported sphingosine metabolism in the original work. Enzymes that have been formerly related to cancer are suggested within the TCA cycle, like fumarate hydratase (Lehtonen et al.,2007;Pithukpakorn et al.,2006;Pollard et al.,2003)succinate dehydrogenase (Ni et al.,2008;Pollard et al.,2003) and aconitase (Singh et al., 2006). Another suggestion is lysosome - lysosomes suffer changes in cancer cells and directly affect apoptosis (Kirkegaard and Jäättelä,2009). Finally, the graph contains several hexokinases, potential targets to disrupt glycolysis, a fundamental need in cancer cells (Kaelin and Thompson,2010). Malaria dataset The metabolites in the last example are related to the distinction between malaria and other febrile ilnesses in (Decuypere et al.,2016). The study files can be found under the MTBLS315 identifier in Metabolights: https://www. ebi.ac.uk/metabolights/MTBLS315. Specifically, the list of KEGG identifiers has been extracted from the supplementary data spreadsheet, using all the possible KEGG matches for the “non malaria” patient group. compounds.malaria <- c( "C05471","C14831","C02686","C06462","C00735","C14833", "C18175","C00550","C01124","C05474","C05469") analysis.malaria <- enrich( compounds = compounds.malaria, data = fella.data, 2This analysis is subject to KEGG release 83.0, from August 17th, 2017. Posterior KEGG releases might alter the reported sub-network
310 the r package fella methods = "diffusion") ## No background compounds specified. Default background will be used. ## Warning in defineCompounds(compounds = compounds, compoundsBackground = compoundsBackground, : Some compounds were introduced as affected but they do not belong to the background. These compounds will be excluded from the analysis. Use ’getExcluded’ to see them. ## Running diffusion... ## Computing p-scores through the specified distribution. ## Done. plot( analysis.malaria, method = "diffusion", data = fella.data, nlimit = 50, vertex.label.cex = vertex.label.cex, plotLegend = FALSE) ## 171 nodes below the threshold have been limited to 50 nodes. ● ● ● ● ● ● ● ●● ● ● ●● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ●● ● ● ● Linoleic acid metaboli... Sphingolipid metabolis... Aldosterone synthesis ... Prostate cancer − Homo... C21−Steroid hormone bi... C21−Steroid hormone bi... 3alpha−hydroxysteroid ... steroid 11beta−monooxy... corticosterone 18−mono... Delta4−3−oxosteroid 5b... sphingomyelin synthase sphingomyelin phosphod... galactosylceramidase glycosylceramidase UDP−alpha−D−galactose:... CDP−choline:N−acylsphi... Sphingomyelin cholinep... Sphingomyelin ceramide... Acyl−CoA:sphingosine N... Cortisol:NAD+ 11−oxido... Cortisol:NADP+ 11−oxid... 11beta,17alpha,21−Trih... steroid,reduced ferred... 17alpha,21−dihydroxy−5... Corticosterone,reduced... 18−Hydroxycorticostero... D−galactosyl−N−acylsph... Galactosylceramide + U... Digalactosylceramide g... Urocortisone:NAD+ oxid... Urocortisone:NADP+ oxi... Urocortisol:NAD+ oxido... Urocortisol:NADP+ oxid... 3alpha,11beta,21−Trihy... 3alpha,11beta,21−Trihy... linoleate:oxygen (8R)−... (8R,9Z,12Z)−8−hydroper... ceramide:phosphatidylc... Cortisol <=> 11beta−Hy... (8R,9Z,12Z)−8−hydroper... Sphingomyelin Cortisol 18−Hydroxycorticostero... Galactosylceramide 17alpha,21−Dihydroxy−5... 11beta,17alpha,21−Trih... 3alpha,11beta,21−Trihy... 8(R)−HPODE In this case, the depicted subnetwork3contains the modules C21-Steroid hormone biosynthesis, progesterone => corticosterone/aldosterone and C21-Steroid hormone biosynthesis, progesterone => cortisol/cortisone, related to the corticosteroids as a main pathway reported in the original text. This is part of the also reported Aldosterone synthesis and secretion; aldosterone is known to show changes related to fever as a metabolic response to infection (Beisel,1975). Another plausible hit in the sub-network is linoleic acid metabolism, as erythrocytes infected by various malaria parasytes can be enriched in linoleic 3This analysis is subject to KEGG release 83.0, from August 17th, 2017. Posterior KEGG releases might alter the reported sub-network
d.2 additional file 2: main vignette 311 acid (Fitch et al.,2000). In addition, the pathway sphingolipid metabolism can play a role in the immune response (Maceyka and Spiegel,2014;Seo et al., 2011). As for the enzymes, 3alpha-hydroxysteroid 3-dehydrogenase (Si-specific) and Delta4-3-oxosteroid 5beta-reductase are related to three input metabolites each and might be candidates for further examination. d.2.5 Conclusions The FELLA R package provides a simple, programmatic and intuitive enrichment tool for metabolomics summary data. Starting from a list of metabolites, FELLA not only pinpoints relevant pathways but also intermediate reactions, enzymes and modules that links the input metabolites to the pathways. The reported entries have a network structure focused on interpretability and new hypotheses generation, giving a richer perspective than classical pathway enrichment tools. This comprehensive layout can also suggest potential enzymes and new metabolites for further study. Finally, FELLA comes equipped with a graphical user interface that promotes its usage to a wider audience and offers interactive sub-network examination. Funding This work was supported by the Spanish Ministry of Economy and Competitiveness (MINECO) [BFU2014-57466-P to O.Y., TEC2014-60337-R and DPI2017- 89827-R to A.P.]. O.Y., A.P. and S.P. thank for funding the Spanish Biomedical Research Centre in Diabetes and Associated Metabolic Disorders (CIBERDEM) and the Networking Biomedical Research Centre in the subject area of Bioengineering, Biomaterials and Nanomedicine (CIBER-BBN), both initiatives of Instituto de Investigación Carlos III (ISCIII). SP. thanks the AGAUR FI-scholarship programme.
312 the r package fella d.2.6 Session info Here is the output of sessionInfo() on the system that compiled this vignette: •R version 3.6.2(2019-12-12), x86_64-pc-linux-gnu •Locale: LC_CTYPE=en_US.UTF-8,LC_NUMERIC=C,LC_TIME=es_ES.UTF-8, LC_COLLATE=en_US.UTF-8,LC_MONETARY=es_ES.UTF-8, LC_MESSAGES=en_US.UTF-8,LC_PAPER=es_ES.UTF-8,LC_NAME=C, LC_ADDRESS=C,LC_TELEPHONE=C,LC_MEASUREMENT=es_ES.UTF-8, LC_IDENTIFICATION=C •Running under: Ubuntu 16.04.6 LTS •Matrix products: default •BLAS: /usr/lib/atlas-base/atlas/libblas.so.3.0 •LAPACK: /usr/lib/atlas-base/atlas/liblapack.so.3.0 •Base packages: base, datasets, graphics, grDevices, methods, parallel, stats, stats4, utils •Other packages: AnnotationDbi 1.46.0, Biobase 2.44.0, BiocGenerics 0.29.2, FELLA 1.5.3, IRanges 2.17.5, knitr 1.22, org.Hs.eg.db 3.8.2, S4Vectors 0.21.24 •Loaded via a namespace (and not attached): assertthat 0.2.1, backports 1.1.4, BiocManager 1.30.4, BiocStyle 2.12.0, biomaRt 2.40.1, Biostrings 2.51.5, bit 1.1-14, bit64 0.9-7, bitops 1.0-6, blob 1.2.0, compiler 3.6.2, crayon 1.3.4, curl 3.3, DBI 1.0.0, digest 0.6.18, evaluate 0.13, GO.db 3.8.2, GOSemSim 2.10.0, grid 3.6.2, highr 0.8, hms 0.5.0, htmltools 0.3.6, httr 1.4.0, igraph 1.2.4.1, KEGGREST 1.24.1, lattice 0.20-38, magrittr 1.5, Matrix 1.2-18, memoise 1.1.0, pillar 1.4.0, pkgconfig 2.0.2, plyr 1.8.4, png 0.1-7, prettyunits 1.0.2, progress 1.2.2, R6 2.4.0, Rcpp 1.0.1, RCurl 1.95-4.12, rlang 0.4.0, rmarkdown 1.12, RSQLite 2.1.1, stringi 1.4.3, stringr 1.4.0, tcltk 3.6.2, tibble 2.1.1, tools 3.6.2, vctrs 0.2.0, xfun 0.6, XML 3.98-1.20, XVector 0.23.2, yaml 2.2.0, zeallot 0.1.0, zlibbioc 1.29.0
d.3 additional file 3: gilt-head bream study 313 d.3 additional file 3: gilt-head bream study d.3.1 Introduction This vignette contains a case study of the effects of environmental contamination on gilt-head bream (Sparus aurata) (Ziarrusta et al.,2018). Fish were exposed over 14 days to oxybenzone and changes were sought in their brain, liver and plasma using untargeted metabolomics. Samples were processed using Ultra-performance liquid chromatography mass-spectrometry (UHPLC-qOrbitrap MS) in positive and negative modes with both C18 and HILIC separation. The mortality of exposed fish was not altered, as well as the brain-related metabolites. However, liver and plasma showed perturbations, proving that adverse effects beyond the well-studied hormonal activity were present. The enrichment procedure implemented in FELLA (Picart-Armada et al., 2017) was used in the study for a deeper understanding of the dysregulated metabolites in both tissues. Building the database At the time of publication, the KEGG database (Kanehisa, Furumichi, et al.,2016) –upon which FELLA is based– did not have pathway annotations for the Sparus aurata organism. It is common, however, to use the zebrafish (Danio rerio) pathways as a good approximation. KEGG provides pathway annotations for it under the organismal code dre, which will be used to build the FELLA.DATA object. library(FELLA) library(igraph) library(magrittr) set.seed(1) # Filter the dre01100 overview pathway, as in the article graph <- buildGraphFromKEGGREST( organism = "dre", filter.path = c("01100")) tmpdir <- paste0(tempdir(), "/my_database") # Make sure the database does not exist from a former vignette build # Otherwise the vignette will rise an error # because FELLA will not overwrite an existing database unlink(tmpdir, recursive = TRUE) buildDataFromGraph( keggdata.graph = graph, databaseDir = tmpdir, internalDir = FALSE, matrices = "none",
314 the r package fella normality = "diffusion", niter = 100) We load the FELLA.DATA object to run both analyses: fella.data <- loadKEGGdata( databaseDir = tmpdir, internalDir = FALSE, loadMatrix = "none" ) Given the 11-month temporal gap between the study and this vignette, small changes to the amount of nodes in each category are expected (see section 2.4Data handling and statistical analyses from the study). Please see the Note on reproducibility to understand why. fella.data ## General data: ## - KEGG graph: ## *Nodes: 10821 ## *Edges: 32013 ## *Density: 0.0002734209 ## *Categories: ## + pathway [163] ## + module [171] ## + enzyme [1021] ## + reaction [5467] ## + compound [3999] ## *Size: 5.9 Mb ## - KEGG names are ready. ## ----------------------------- ## Hypergeometric test: ## - Matrix not loaded. ## ----------------------------- ## Heat diffusion: ## - Matrix not loaded. ## - RowSums are ready. ## ----------------------------- ## PageRank: ## - Matrix not loaded. ## - RowSums not loaded. Note on reproducibility We want to emphasise that each time this vignette is built, FELLA constructs its FELLA.DATA object using the most recent version of the KEGG database. KEGG is frequently updated and therefore small changes can take place in the knowledge graph between different releases. The discussion on our findings was written at the date specified in the vignette header and using the KEGG release in the Reproducibility section.
d.3 additional file 3: gilt-head bream study 315 d.3.2 Enrichment analysis on liver tissue Defining the input and running the enrichment Table 1from the main body in (Ziarrusta et al.,2018) contains 5KEGG identifiers associated to metabolic changes in liver tissue and 12 in plasma. Our first enrichment analysis with FELLA will be based on the liver-derived metabolites. Also note that we use the faster approx = "normality" approach, whereas the original article uses approx = "simulation" with niter = 15000 This is not only intended to keep the bulding time of this vignette as low as possible, but also to demonstrate that the findings using both statistical approaches are consistent. cpd.liver <- c( "C12623", "C01179", "C05350", "C05598", "C01586" ) analysis.liver <- enrich( compounds = cpd.liver, data = fella.data, method = "diffusion", approx = "normality") ## No background compounds specified. Default background will be used. ## Running diffusion... ## Computing p-scores through the specified distribution. ## Done. All the metabolites are successfully mapped: analysis.liver %>% getInput %>% getName(data = fella.data) ## $C12623 ## [1] "trans-2,3-Dihydroxycinnamate" ## [2] "(2E)-3-(2,3-Dihydroxyphenyl)prop-2-enoate" ## ## $C01179 ## [1] "3-(4-Hydroxyphenyl)pyruvate" "4-Hydroxyphenylpyruvate" ## [3] "p-Hydroxyphenylpyruvic acid" ## ## $C05350
322 the r package fella ## [13] dre01040--1.14.19.3 dre01212--1.14.19.3 dre00564--1.1.5.3 ## [16] dre00564--2.3.1.15 dre00564--2.7.8.29 dre00564--2.7.8.5 ## [19] dre00564--3.1.1.32 dre00591--3.1.1.32 dre00592--3.1.1.32 ## + ... omitted several edges Examining the pathways Figure 3from the original study is a holistic view of the affected metabolites found in plasma, based on literature and on an analysis with FELLA. The 11 metabolites are depicted within their core metabolic pathways. We will check whether FELLA is able to highlight them, by first showing the reported metabolic pathways: tab.plasma[tab.plasma$Entry.type == "pathway", ] ## KEGG.id Entry.type KEGG.name ## 1 dre00062 pathway Fatty acid elongation - Danio rerio (zebrafis... ## 2 dre00260 pathway Glycine, serine and threonine metabolism - Da... ## 3 dre00564 pathway Glycerophospholipid metabolism - Danio rerio ... ## 4 dre00591 pathway Linoleic acid metabolism - Danio rerio (zebra... ## 5 dre00592 pathway alpha-Linolenic acid metabolism - Danio rerio... ## 6 dre01040 pathway Biosynthesis of unsaturated fatty acids - Dan... ## 7 dre01212 pathway Fatty acid metabolism - Danio rerio (zebrafis... ## p.score ## 1 1.000000e-06 ## 2 1.934171e-06 ## 3 1.080598e-05 ## 4 2.639328e-02 ## 5 1.000000e-06 ## 6 1.000000e-06 ## 7 2.448355e-05 And then comparing against the ones in Figure 3: path.fig3 <- c( "dre00591",# Linoleic acid metabolism "dre01040",# Biosynthesis of unsaturated fatty acids "dre00592",# alpha-Linolenic acid metabolism "dre00564",# Glycerophospholipid metabolism "dre00480",# Glutathione metabolism "dre00260" # Glycine, serine and threonine metabolism ) path.fig3 %in% V(g.plasma)$name ## [1] TRUE TRUE TRUE TRUE FALSE TRUE All of them but Glutathione metabolism are recovered, showing how FELLA can help gaining perspective on the input metabolites.
[Document text truncated for crawler view.]