Discovering gene networks

Ant´onio Jos´e Santos Freitas Gon¸calves Discovering Gene Networks Department of Computer Science Faculty of Sciences at University of Porto November 2012 Ant´onio Jos´e Santos Freitas Gon¸calves Discovering Gene Networks Dissertation presented to Faculty of Sciences at University of Porto in partial fulfillment of the requirements for the degree of Master in Computer Science Scientific Advisor: V´ıtor Manuel de Morais Santos Costa, PhD Scientific Co-Advisor: Irene May Lin Ong, PhD Department of Computer Science Faculty of Sciences at University of Porto November 2012 to my parents and girlfriend without whom, none of this would have been possible... 5 Acknowledgements First and foremost, I wish to thank my scientific advisor: V´ıtor Santos Costa. I learned a lot from him. I am grateful for his support, guidance and patience especially over this past year. A special word of thanks to Irene Ong who also supported and helped a lot with the guidance of this work. A very special thank and appreciation to my parents, who always helped, believed and encouraged me in every way they could so that I could have a proper education and also succeed in my professional life. Also, a very special thanks and appreciation to my girlfriend Raquel, who never let me give up and always believed in me, and gave me her inconditional support every day. Thanks also to Jeffrey Lewis who provided data and consulted. A thanks to all my Masters teachers at DCC FCUP who inspired me to pursue knowledge discovery and research. I also wish to thank the Fundac˜ao para a Ciˆencia e a Tecnologia. This work was supported with the FCT grant PTDC/EIA-EIA/100897/2008. A word of appreciation to my friends and colleagues for the support. Last, but not least, a thanks to my crazy cats Nero and Sasha who were always able to make me laugh when my work was not going very well at that time. 7 Abstract Transcriptional regulation plays an important role in every cellular decision. Gaining an understanding of the dynamics that govern how a cell will respond to diverse environmental cues is difficult using intuition alone. In this work, logic-based regulation models based on state-of-the-art work on statistical relational learning are introduced, and evaluated on time-series gene expression data of the Hog1 pathway. Results show that plausible regulatory networks can be learned from time series gene expression data using a probabilistic logical model. Hence, network hypotheses can be generated from existing gene expression data for use by experimental biologists. Keywords:-Gene Regulation, Network/Pathway Analysis, Statistical Relational Learning Palavras-Chave: Regula¸c˜ao de Genes, An´alise de Redes/Caminhos, Aprendizagem Relacional Estat´ıstica 9 4.1 On the left a Decision tree; In the middle, a corresponding truth table; On the right a corresponding BDD . . . . . . . . . . . . . . . . . . . . 58 4.2 ReductionofaBDD ............................ 59 4.3 Both OBDDs are reduced, but they equivalent. Although only the left one is ordered [x, y, z]. ........................... 60 4.4 Execution of the reduce algorithm..................... 61 4.5 Two arguments for a call apply(+, Bf, Bg)................. 62 4.6 The recursive call structure for apply for example in Fig.4.5. . . . . . . 63 4.7 The result of apply(+, Bf, Bg)........................ 63 4.8 A simple directed graph, where each edge has a probability of being true. 64 4.9 A BDD that computes the total probability of the path ae. . . . . . . . 66 5.1 AsimpleBDD................................ 72 5.2 Implementation Method . . . . . . . . . . . . . . . . . . . . . . . . . . 77 5.3 Same-step Correlation Generated Hog1 Promoter Network . . . . . . . 78 5.4 Learning Curves for Test Data Mean Square Error, resp for VV (∆ ⇒ ∆), VL (∆ ⇒L) and LL (L⇒L). .................... 79 5.5 Speedup graph based on Table 5.5 results . . . . . . . . . . . . . . . . . 81 5.6 Speedup graph based on Table 5.6 results . . . . . . . . . . . . . . . . . 82 5.7 Speedup graph based on Table 5.7 results . . . . . . . . . . . . . . . . . 83 16 Acronyms AIK - Akaike’s Information Criterion BDD - Binary Decision Diagram BN - Bayesian Network COD - Coefficient of Determination CTL - Computation Tree Logic CUDD - CU Decision Diagram DBN - Dynamic Bayesian Network DNA - Deoxyribonucleic acid dsRNA - Double-stranded Ribonucleic acid EM-algorithm - Expectation–Maximization algorithm GE - Gene Expression GRN - Gene Regulatory Network GP - Genetic Programming Lac - Lactose LMS - Least Mean Square LNA - Linear Noise Approximation LTM - Linear Transcription Model MAPK - Mitogen-activated protein kinases Mg2+- Magnesium 17 miRNA - Micro Ribonucleic acid mRNA - Messenger Ribonucleic acid MWSLE - Minimum Weight Solutions to Linear Equations Na+- Sodium NaCl - Sodium Chloride ncRNA - Non-coding Ribonucleic acid OBDD - Ordered Binary Decison Diagram ODE - Ordinary Differential Equation PBN - Probabilistic Boolean Network PDE - Partial Differential Equation PIN - Protein Interaction Network PRM - Probabilistic Relational Model PSN - Protein Signaling Network RNA - Ribonucleic acid RNAseq - Whole Transcriptome Shotgun Sequencing rRNA - Ribosomal Ribonucleic acid SDE - Stochastic Differential Equation siRNA - Small interfering Ribonucleic acid SRL - Statistical Relational Learning SVD - Singular Value Decomposition Tet - Tetracycline tmRNA - Transfer-messenger Ribonucleic acid tRNA - Transfer Ribonucleic acid YAP - Yet Another Prolog 18 Chapter 1 Introduction If you don’t know where you are going, you’ll end up someplace else. Yogi Berra 19 FCUP Discovering Gene Networks 21 1.1 Introduction Many major cellular decisions involve changes in transcriptional regulation. With the advent of high-throughput technologies and advanced measurement techniques molecular biologists and biochemists are rapidly identifying components of these networks and determining their biochemical activities, but understanding these complex multicomponent networks that govern how a cell will respond to diverse environmental cues is difficult using intuition alone. In this work, our goal is to build a probabilistic logical model that can aid in uncovering the structure and dynamics of such networks and how they regulate their targets. Gaining insight into transcriptional regulation is important not just for understanding the fundamental biological processes, but also for advancing research. A cell responds to environmental changes by detecting molecules that bind to receptors on the surface of the cell and transmits this information to proteins within the cell by activating a cascade of molecular events. These signaling networks are typically studied by measuring molecular events after treatments that stimulate or perturb key elements in the network at the mRNA or gene level. Combining these measurements with phenotypic response enables the study of these treatments on the architecture and function of the underlying signaling networks as well as the relationship between the network behavior and phenotypic response [1]. In this work, to infer the architecture and function of the underlying signaling network of budding yeast, we decided to focus on the pathways activated by MAPK1Hog1 during osmotic stress response. According to previous results [2], this pathway interacts with the general stress (Msn2/Msn4) pathways, so we consider the genes belonging to both pathways in the following study. Despite the challenge of inferring genetic regulatory networks from gene expression data, various computational models have been developed for regulatory network analysis. Examples include approaches based on logical gates [3, 4], and probabilistic approaches, often based on bayesian networks [5]. On one hand, logic gates provide a natural, intuitive way to describe interactions between proteins and genes. On the other hand, probabilistic approaches can handle incomplete and imprecise data in a 1Mitogen-activated protein (MAP) kinases are serine/threonine-specific protein kinases belonging to the CMGC (CDK/MAPK/GSK3/CLK) kinase group. These kinases regulate gene expression, proliferation, differentiation, mitosis, cell survival - among many others. 22 FCUP Discovering Gene Networks very robust way. Our main contribution is in introducing a model that combines the two approaches. Our approach is based on the probabilistic logic programming language ProbLog [6, 7]. In this language, we can express true logical statements (expressed as true rules) about a world where there is uncertainty over data, expressed as probabilistic facts. In the setting of gene expression, this corresponds to establishing: •a set of true rules describing what are the possible interactions existing in a cell; •a set of uncertain facts describing which possible rules are applicable to a certain gene or set of genes. Given time-series gene expression data, we want to choose the probability parameters that best describe the data. Our approach is to reduce this problem to an optimization problem, and use a gradient ascent algorithm to estimate a local solution [8] in the style of logistic regression. We further contribute an efficient implementation to this algorithm that computes both probabilities and gradients through binary decision diagrams (BDD). We evaluate our approach by using it to study expression data on an important geneexpression pathway, the Hog1 pathway [2]. It is well known that under conditions of osmotic stress, the protein kinase Hog1, and the paralogous proteins Msn2 and Msn4 interact to create a response that involves the expression of a large number of proteins. We model these pathways by also incorporating the two transcription factors activated by Hog1:Hot1 and Sko1. 1.2 Thesis Roadmap This thesis is divided into a total of six chapters. In Chapter 1 we present a small introduction to the biological process and the problem of discovering gene networks. In Chapter 2 we explain some important concepts for understanding the work in this thesis, we talk about DNA, RNA, gene regulatory networks and also about expression data, the data we use for our models. Next, in Chapter 3 we talk about related work in the area, presenting some of the techniques that can be used for creating gene regulatory networks. Chapter 4 is a brief introduction to logic programming and also concepts that are necessary for understanding the processes behind the calculations. A short introduction to Problog is also given. In Chapter 5 we explain how to represent FCUP Discovering Gene Networks 23 gene networks with Problog, how we perform our experiments and develop our models and we show the results we obtained. Last, Chapter 6 presents the conclusion we took from our work and we also talk about the future work we intend to perform. Chapter 2 Biology Background Every biologist has at some time asked ’What is life?’ and none has ever given a satisfactory answer. Science is built on the premise that Nature answers intelligent questions intelligently; so if no answer exists, there must be something wrong with the question. Albert Szent-Gy¨orgyi 25 32 FCUP Discovering Gene Networks the adjacent phosphodiester bond to split the backbone [19]. It is known to exist about 100 naturally occurring modified nucleosides, although the specific role of these modifications in RNA are not yet fully understood. Many of the post-transcriptional modifications occur in highly functional regions, which implies that they are important for normal function. Frequently, the functional form of single stranded RNA molecules require a specific tertiary structure which is provided by the secondary structural elements (hydrogen bonds within the molecule). The result of this are the recognizable secondary structure known as hairpin loops,bulges and also internal loops, as shown in Fig. 2.7. Figure 2.7: RNA Secondary Structures RNA is charged, so in order to stabilize many secondary and tertiary structures, metal ions such as Mg2+ and Na+are needed as loop information is unfavorable due to backbone charge-charge repulsion. These two metal ions can increase the loop flexibility by neutralizing the phosphate charges, causing the loop formation to be less unfavorable. The result of increasing Mg2+ and Na+is a decrease on the energy cost for loop formation [20]. 2.2.2 Synthesis RNA synthesis, or transcription is under the control of the enzyme RNA polymerase. The first step that this enzyme takes is to find the start of the gene on the coding FCUP Discovering Gene Networks 33 strand of the DNA, as the DNA has lots of genes strung out along the coding strand (the enzyme has to pick the right strand and identify the beginning of each gene). This is done by recognizing and binding with one or more short sequences of bases (also known as promoter sequences, Fig. 2.8) ”upstream” of the start of each gene. The transcription process (Fig. 2.9) is composed of the following steps: Figure 2.8: Promoter Sequences •The DNA double helix is unwound due to the helicase activity of the enzyme •Enzyme progression along the template strand following the 3’ to 5’ direction •Synthesis of a complementary RNA molecule Figure 2.9: RNA Synthesis The end of RNA synthesis is indicated by the DNA sequence. RNAs are often modified by enzymes after transcription. There are also a number of RNA-dependent RNA polymerases that use RNA as their template for synthesis of a new strand of RNA. For instance, a number of RNA viruses use this type of enzyme to replicate their genetic material. 34 FCUP Discovering Gene Networks 2.2.3 Types There are several types of RNA, as discussed before: •Messenger RNA (mRNA) - molecule in cells that carry codes from the DNA to the sites of protein synthesis (the ribosomes). Information in DNA cannot be decoded directly into proteins, and thus it is first transcribed, or copied, into mRNA. Each molecule of mRNA encodes the information for one protein (or more than one protein in the case of a bacteria), with each sequence of three nitrogen-containing bases in the mRNA. The mRNA takes the copy of the blueprint to the ribosome where it is used to build the protein. Figure 2.10: Messenger RNA •Ribosomal RNA (rRNA) - these molecules form the structural and functional components of ribosomes, the subcellular units responsible for protein synthesis [21]. This type of RNA is the catalytic component of the ribosomes. rRNA Figure 2.11: RNA abundance in cells FCUP Discovering Gene Networks 35 constitutes approximately 80% to 85% of the total RNA in a cell (shown in Fig. 2.11). In eukaryotes, rRNA synthesis occurs in the nucleolus, a specialized structure within the nucleus. •Transfer-messenger RNA (tmRNA) - molecule of RNA (found in many bacteria and plastids) that has dual functions as both a transfer RNA and a messenger RNA. As a tRNA, it recognizes and binds ribosomes stalled by aberrant mRNAs with the help of its protein partner SmpB. As an mRNA, it adds a degradation tag to protein fragments, targeting them for proteolysis. Two types of tmRNAs are known: single-chain tmRNAs and two-piece tmRNAs. Figure 2.12: Transfer-messenger RNA •Non-coding RNA (ncRNA) - functional molecule that is not translated into protein. ncRNAs have been shown to regulate important biological processes that support normal cellular functions, but relatively little is known about the general structure, function, and transcriptional control of ncRNAs, and even less about their potential functions as a group or as a single entity [22]. •Other types of RNA are the following: Double-stranded RNA (dsRNA) - RNA with two complementary strands, similar to the DNA found in all cells. dsRNA forms the genetic material of some viruses (dsRNA viruses). dsRNA such as viral RNA or siRNA can trigger RNA interference in eukaryotes, as well as interferon response in vertebrates. MicroRNA (miRNA) - small type of dsRNA molecules that regulate translation in eukariotic cells, it is related to RNAi. Unlike other small RNAs, the genes for miRNA are transcribed by RNA polymerase II. MicroRNAs are not only used in eukariotic cells but also by some of the more complex viruses that infect them. Short-interfering RNA (siRNA) - short dsRNA fragments that are known to bound by the the RNA-induced silencing complex1(RISC) [23]. It is class of double-stranded RNA molecules 1RNA-induced silencing complex, or RISC, is a multiprotein complex that incorporates one strand of a small interfering RNA (siRNA) or micro RNA (miRNA). RISC uses the siRNA or miRNA as a 36 FCUP Discovering Gene Networks that play many roles, being the most notable the RNA interference (RNAi) pathway. 2.3 Gene Regulatory Networks Every cell is a complex processor of information, as it is able to integrate and respond to multiple signals in a robust way. The mechanism that cells use to achieve this are remarkable and in a large number. One can make an analogy to electrical circuits, that is, we can decompose the high complexity involved in cellular response into modules that connect to each other by input and output signals [24]. We can take this comparison a little further, as genetic network engineers manipulate living organisms using the biological equivalent of transistors and inverters. A possible description of gene networks is circuits of interconnected functional modules, each consisting of specialized interactions (information flow) between proteins, RNA, DNA and small molecules. The importance of these modules (components) can be exemplified by their Figure 2.13: Example of an architecture of the inducible gene network use in networks that function in cells that have a higher complexity. An advantage of prokaryotic components such as repressors and their operating sites is that they can be transplanted into eukaryotic cells without any loss of their binding or function specificity. A second advantage is the avoidance of unwanted interference with the expression of non-targeted genes [24]. There has been a large development of switch systems based on the Tet or Lac repressor in a large number of organisms ([25], [26], [27], [28]). This development sets the stage for the construction of more complex gene networks. template for recognizing complementary mRNA. When it finds a complementary strand, it activates RNase and cleaves the RNA. FCUP Discovering Gene Networks 37 In order to synthesize proteins that carry out specific functions in the cell, this information is extracted through a process called gene expression. On the other hand, gene regulation describes all the cellular processes which control the expression of proteins. This regulation can occur at the following steps of the DNA information extraction: •Transcription initiation - start of the binding of RNA polymerase to the promoter in DNA •Translation - the mRNA produced by transcription is decoded by the ribosome in order to produce a polypeptide, or a specific amino acid chain, that will later fold into an active protein. •Modifications of mRNAs - modified to remove certain stretches of non-coding sequences called introns. •DNA packing - process in which DNA and associated proteins are formed into a compact, orderly structure. •Post-translational modifications of amino-acid sequences - this modification extends the range of functions of the protein by attaching it to other biochemical functional groups, making structural changes, or changing the chemical nature of an amino-acid. There are many different gene expression patterns in a cell, they depend of factors like the state of the cell, nutrition, environment or even the cell type. The process of gene expression consists of translation and transcription, with the regulation of this expression taking place at different steps, but usually more during transcription (truest when talking about procaryotes like bacteria). As exons and introns are not distinguished in their mRNA molecules, an alternative splicing (one of the more relevant regulation processes taking place after the transcription step [29])takes place, leading to different proteins from the same mRNA molecule. In order to build a model for gene regulation, one needs to specify the interactions (the ones that the dynamic behavior of gene expression can be explained by them) that should be captured by the model. The transcription of a gene can be activated or inhibited by transcription factors which can bind to a site on the DNA near the promoter (the place where the following gene transcription is started) in a way that the transcription factor that binds inhibits the 38 FCUP Discovering Gene Networks Figure 2.14: Gene Regulation: The activation or inhibition of a gene transcription by one or several transcription factors. start of the transcription. In order to activate a transcription, transcription factors can absorb other molecules that bind to that particular site. We can see different possibilities of inhibition and activation in Fig. 2.14. The product of a gene can regulate its own transcription or of another gene, it can also regulate the expression with the help of other products of genes. Regulatory proteins and their binding sites are also modular in that different domains from different proteins can be combined to yield hybrid proteins of novel function. The state in which a living organism is at any certain point of time is not only described by its genome, but also by its set of expressed regulatory genes and its concentration levels of the corresponding gene products. All the possible phenotypic states of a cell correspond to distinct gene expression patterns [30]. Figure 2.15: Processes in a Gene Regulatory Network To sum all up, we can say that gene regulatory networks (GRNs) are interacting DNA-encoded regulatory subsystems in the Genome that have the function of coordinating the inputs from activators and repressors (Transcription Factors) during cell differentiation, development or even in response to environmental causes (natural or FCUP Discovering Gene Networks 39 not). GRNs funtion to specify expression of particular sets of genes for specific times, locations and conditions. 2.4 Expression Data In order to sequence the human genome, methods for measuring the expression levels of single and all genes simultaneously in a genome had to be created. An existing challenge of Computational Biology is the analysis of thousands of measurements of one cell state in order to retrieve the useful information from it. To get these measures, Microarray Technology 2is applied. By using this technology, the researchers can observe the dependency of gene expression on different states of a cell and also on different environmental factors. 2.4.1 Data Measurement The description of any model gets better as more relevant information becomes available. This seems trivial as huge amounts of information exist nowadays, but unfortunately for gene regulation that is not true, some of the data available is only for selected targets, which some of the times are not satisfactorily accessible and most of the data is very noisy. The advance in DNA-microarray technology permits to monitor thousands of genes in one experiment by measuring mRNA concentrations in a cell [31]. Each data point produced by a DNA microarray experiment represents the ratio of expression levels of a particular gene. The result, from an experiment with n genes on a single chip, is a series of n expression-level ratios. Typically, the numerator of each ratio is the expression level of the gene in the varying condition of interest, whereas the denominator is the expression level of the gene in some reference condition. The data from a series of m such experiments may be represented as a gene expression matrix, in which each of the n rows consists of an m-element expression vector for a single gene. The expression measurement is positive if the gene is induced (turned up) with respect to the reference state and negative if it is repressed (turned down). Every spot on the microarray has millions of copies of one probe in order 2An array is an orderly arrangement of samples where matching of known and unknown DNA samples is done based on base pairing rules. An array experiment makes use of common assay systems such as microplates or standard blotting membranes. The sample spot sizes are typically less than 200 microns in diameter usually contain thousands of spots 40 FCUP Discovering Gene Networks Figure 2.16: A general scheme for a microarray. to measure not only the existence of specific mRNA in the test material, but also to measure its amount, making it interesting if time series are produced, such that the difference between two time points can give a hint for genes that are transcribed together or even a hint for some regulatory interactions. Chapter 3 Related Work Knowing is not enough; we must apply. Willing is not enough; we must do. Johann Wolfgang von Goethe 41 48 FCUP Discovering Gene Networks hx1, . . . , xp,¬y1,...,¬yqi. In order to identify a GRN from observed global states they investigate the number of experiments together with the cost of each experiment. This work lead to two important problems obtained from the identification of GRNs: •consistency (checking if a network G0= (V0, F0) coincides or not with the underlying GRN) •stability (Gis stable if there exists a global state consistent with all gene regulation rules) of the network Also, in their work they have derived upper and lower bounds on the required perturbations for the boolean networks. Another approach using boolean networks is the work developed by Ideker et al. [56]) in which they present two methods (Predictor and Chooser) for inferring a genetic network given gene expression measurements. The Predictor method is used to infer hypothetical boolean networks consistent with a profile that was generated previously by exposing the network of interest to a series of genetic or biological perturbations. After the inference of several networks, the Predictor method returns the ones that are most parsimonious (the ones having the fewest number of interactions). The next step is done by the Chooser method that adds an additional perturbation experiment in order to discriminate along the set of previously obtained hypothetical networks. These perturbations are added in a clever way by using a function that is based on entropy to optimally reduce the number of remaining hypothetical networks. One other feature of these two methods is that they can be used interactively and iteratively. Standard boolean networks have a very salient limitation, their inherent determinism. One can look at this from two points of view, the conceptual point of view or the empirical point of view. From the former, one must bear in mind that is likely that the regularity of genetic function and interaction known to exist is not due to hard-wired logical rules, as for the latter, it takes in account that one logical rule per gene may to incorrect results when these rules are being inferred from gene expression measurements. To address these considerations Shmulevich et al. [57] introduced a new model class, the Probabilistic Boolean Networks (PBNs), which have the properties of boolean networks and can cope with uncertainty in the data and model selection. The basic idea of their work is to extend the boolean network in order to accommodate more than one possible function for each node. For this, they have a set Fi={f(i) j}j= 1, . . . , l(i) that corresponds to each node xiand where each f(i) jis a possible function that determines the value of gene xiand l(i) is the FCUP Discovering Gene Networks 49 number of possible functions for gene xi. For inferring they use a method based on the Coefficient of Determination (COD) which produce a number of candidate predictors for each target gene. In this model, the approach was to probabilistically ’create’ good predictors so that each predictor contribution is proportional to its determinative potential, this is close to our own approach. Also in this field of logical models, but on a more formal approach, we have the work by Bernot et al. [58] where they provide a formal way to treat temporal properties of biological regulatory networks, expressed in computational tree logic leading to the possibility of building all the models that satisfy a set of given temporal properties. Their work allows biology to take advantage from all the available formal methods from computer science, temporal properties can be checked against models using Computation Tree Logic (CTL) and model checking. Still in the formal approach, we have the work by Batt et al. [59] that validates models of GRNs by addressing the challenge matching model predictions and experimental data, taking also in account a reliable and efficient comparison between the observations and the predictions. The qualitative modeling and simulation method that they use is based in a refinement of previous work [60]. What is new about their work is that they use model-checking techniques to attend the problem that the state transition graphs that are generated by qualitative simulation may become very large in a way that it may be prohibitive for interesting biological networks. An interesting work involving logical analysis it the one presented by Thomas et al. [61] where they present a logical method for the analysis of the complex dynamics of regulatory networks in terms of feedback circuits. The feature that distinguishes the most their work from others is that the logical method introduced by them is fully asynchronous, that is, current variables are discrete, but time is continuous. Besides this work presented by Thomas et al., they also have published other interesting ones related to logical models, their analysis and creation [62] [63] [64][65]. In recent work by Handorf et al., we are presented with a new method for an automated generation of boolean network models from curated mechanistic network databases. This straight forward method translates into a good gain that is even more important in the context of the fast growing amount of interactions available in these databases. Despite of this great advantage, there is also a drawback that is common to all boolean approaches, interactions strengths and concentrations can not be properly covered by TRUE and FALSE values. 50 FCUP Discovering Gene Networks Further work has been developed in the area of GRNs using logic models [66] [67] [68] and even logic tools to model and analyze signaling pathways [69]. 3.3 Other approaches At the other extreme, a very active field in computing graphical representations of biological networks [70] [71] through literature analysis or through identification of correlations in high-throughput data has emerged. In these graphs, termed protein interaction networks (PINs or interactomes) or protein signaling networks (PSNs), genes and proteins are represented by nodes and potential interactions by edges (links). The edges can be directional or not and signed (inhibitory/activating) or not and typically represent a wide range of interaction modes from direct physical binding to correlated gene expression or integrated database entries. Graphs are an attractive way to summarize diverse relationships among large numbers of biomolecules across multiple organisms, but they are not executable per se and cannot be used to compute input-output relationships. Moreover, network graphs rarely take into account dynamic changes in signaling activities, cell type-specific biochemistry, or contextdependent variations. Despite the difficulty of deciphering genetic regulatory networks from microarray data, numerous approaches to the task have been quite successful. Friedman et al. [5] were the first to address the task of determining properties of the transcriptional program of S. cerevisiae (yeast) by using Bayesian networks (BNs) to analyze gene expression data. Pe’er et al. [72] followed up that work by using BNs to learn master regulator sets. Other graphical approaches are Tanay and Shamir [73], Chrisman et al. [74]). The former implemented a new software platform called Genesis to enable analysis of available transcription profile data sets and target pathways. The latter present a Bayesian framework which combines information from several different sources and makes the correct causal inferences with small sample sizes. The methods above can represent the dependence between interacting genes, but they cannot capture causal relationships. Pe’er et al. [72] ingeniously proposed the use of microarray experiments in which specific genes have been deleted (knockout) in yeast to obtain causality. The use of perturbations such as gene deletion mutants can allow the BN learning algorithm to learn a directed edge that suggests direct causal influence. This approach of combining observational and interventional data delivered promising results. Unfortunately, a complete library of gene knockouts are FCUP Discovering Gene Networks 51 not yet available for organisms other than yeast. The advent of small interfering RNA (siRNA) can be used to reduce the expression of a specific gene in organisms other than yeast, however, siRNA does not guarantee complete silencing of the gene. Ong et al. [75], proposed that the analysis of time series gene expression microarray data using Dynamic Bayesian networks (DBNs) could allow to learn potential causal relationships. DBNs are usually based on discrete models [76] [77] [78] or on continuous time models [79]. Another approach, is the work by Perrin et al.[80]. In their work they used penalized likelihood maximization in EM-algorithms to learn the parameters for a DBN. Dojer et al. [81] also apply DBNs, but this time in the context of perturbation experiments. With the incorporation of this type of data they are able to check that the quality of inferred networks dramatically improves. A different approach is presented by Yeung et al. [82] who propose a scheme to reverseengineer gene networks using Singular Value Decomposition (SVD) to construct a set of candidate solutions and afterwards applies robust regression to identify the solution with the smallest number of connections. Last, we refer the reader to Madeira et.al [83], a work where biclustering algorithms for biological data analysis are described and analyzed. Biclustering algorithms have been used for some time now, the first time the term was used in gene expression data analysis was by Cheng et. al [84]. The main difference between these type of algorithms and the simple clustering algorithms is that the second ones can be applied to columns or rows of the data matrix (each at a time), but the first algorithms can perform clustering on the two dimensions (rows and columns) at the same time - meaning that with biclustering produces a local model while clustering produces a global model. With this we can say that the goal of biclustering algorithms is to identify genes subgroups and conditions subgroups by performing clustering on the rows and columns of the gene expression data matrix at the same time. Contrary to the clustering algorithms, biclustering is able to find sets of genes that have similar activity under a specific set of conditions. Chapter 4 Introduction to Logic Programming Logic is not a body of doctrine, but a mirror-image of the world. Logic is transcendental. Ludwig Wittgenstein 53 FCUP Discovering Gene Networks 55 In this chapter we will talk about Propositional logic (also known as sentential logic, is that branch of logic that studies ways of combining or altering statements or propositions to form more complicated statements or propositions) as a small introduction prior to talking about Binary Decision Diagrams (data structures for representing the semantics of a formula in propositional logic) and Ordered Binary Decisions Diagrams of which we show three algorithms that can be performed on them (Reduce, Apply, Restrict). Next we talk about Problog, giving the fundamental ideas of how it works and how it is implemented. 4.1 Logic The objective of propositional logic is to model human thinking. Starting from declarative phrases (propositions), which can be true (T) or false (F) we construct propositions using connectives such as or (∨), and (∧), not (¬), if...then... (→). Let us consider the following phrases, and an interpretation: •Penguins are birds T •Africa is a continent T •1+1=4F •A triangle has 7 sides F •1<7T From this we may deduce that: •Penguins are birds and Africa is a continent T as it is a conjunction of Tpropositions •A triangle has 7 sides or 1<7T as it is a disjunction of propositions in which one of them is T •not 1+1=4T as it is a negation of a Fproposition 56 FCUP Discovering Gene Networks Connective Symbols Arity Equivalent Symbols Conjunction ∧2 & Disjunction ∨2 + Implies →2⊃ Negation ¬1 ! Table 4.1: Propositional logic symbols In table 4.1 we can see some symbols of propositional logic: Logic programming was created to help programming languages that were more readable and expressive. These programming languages borrow expressive power from mathematical logic. The most popular logical language is Prolog (the first Prolog system was developed in 1972 by Colmerauer with Philippe Roussel). Semantics assign meaning to programs [85]. We have two types of semantics for logic programs: •Declarative semantics - based on the standard model-theoretic semantics of firstorder logic, it describes what we want to use logic programming for. •Operational Semantics - is a way of describing procedurally the meaning of a program [85]. It represents a procedure to satisfy the list of objectives in the context of a given program. The output of this procedure is the truth value from the list with the objectives with their respective instantiation of variables. The Prolog procedure allows for the automatic return (backtracking) in order to examine new alternatives. The operational semantics are based on Herbrand Model, Universe, Base and also Interpretation using an example in order to better explain the subjects. •Herbrand Universe Let Lbe a first order language. The Herbrand Universe UL for Lis the set of all the basic terms that can be obtained from the constants and functions in L. If Ldoes not have any constants, we add a constant for the generation of basic terms. ex. Let us consider the following logic program P1: p(a). FCUP Discovering Gene Networks 57 p(b). q(X):-p(X). The first order language are program P1clauses. Therefore, P1’s Herbrand Universe is UL = {a,b}. •Herbrand Base Let Lbe a first order language. Herbrand Base BL for Lis the set of all the basic atoms which can be obtained from using the predicates of Lwith the basic terms of its corresponding Herbrand Universe as arguments. ex. The Herbrand Base for our logic program P1is BL = {p(a), p(b), q(a), q(b)}. •Herbrand Model Let Lbe a first order language and Sa set of closed formulas of L. A Herbrand Model for Sis the Herbrand Interpretation which is a model for S. ex. Let us consider our logic program P1and also S(the set formed by the program clauses). Rewriting the program we get: p(a). q(b). q(X):-p(X). The domain of Herbrand’s pre-interpretation for this program is {a, b}, we do not have functions. If we consider p, q Herbrand’s Interpretation predicates, this Interpretation we just built is a Herbrand Model for our Program P1. The reason for this is that all program clauses are true in this interpretation. A logic program Pis a set of clauses. If this program has a model then it will have a Herbrand Model. •Herbrand Interpretation Let Lbe a first order language. An interpretation of Lis a Herbrand Interpretation, if the following rules are satisfied: 1. The interpretation domain is the Herbrand Universe, UL 2. The constants in Lare assigned to themselves in UL 3. If fis a n-ary function in L, then we assign fthe (ULn) mapping in UL defined by (t1,· · · , tn)→f(t1,· · · , tn) 64 FCUP Discovering Gene Networks Figure 4.8: A simple directed graph, where each edge has a probability of being true. As a straightforward example of ProbLog, consider the directed graph in Figure 4.8. Each edge has a probability of being true. ProbLog represents edges as probabilistic facts: 0.2::edge(a,b). 0.5::edge(a,c). 0.7::edge(b,c). .... Notice that the probability of each edge being true is independent of all other edges. ProbLog allows one to specify intensional logic programs that describe true relationships that hold in the program. As an example, the next two rules define a path relation between any two nodes: path(N,N). path(N,E) :- edge(N,M), path(M,E). The first rule says that there is a path between node Nand Eif N=E. The second rule defines path recursively: there is a path from Nto E, if from None can reach Mand from Mthere is a path to E. Given a logic program and the set of probabilistic facts, a ProbLog engine allows one to reason about the probability of a relation, which is just the joint probability of all proofs of the relation. More precisely: •The probability of a proof is the product of the probabilities of the facts appearing in the proof. FCUP Discovering Gene Networks 65 •The probability of a relation being true is given by the probability of the union of all proofs. Imagine we want to find out P r(path(a,d)) in the example above. One first proof is obtained by using the second clause, and calling edge(a,b) and path(b,d). The first call is a probabilistic fact. We assume it succeeds. The second call generates edge(b,d), another probabilistic fact, and path(d,d), that succeeds with its first clause. The proof is true if the probabilistic facts are true. Hence, the probability for this proof is the product of the probabilities of its (independent) probability facts: 0.2×0.7. The second proof is obtained by considering the facts edge(a,c) and edge(c,b), and has probability 0.5×0.2. The two proofs are disjoint, hence P r(path(a,d)) = 0.2×0.7 + 0.5×0.2 = 0.24. Computing the total probability is more interesting if different paths have a common edge. As an example, consider P r(path(a,e)). There are three proofs, that we write concisely as paths abde,acde, and abe. Notice that the first proof, abde, shares the edge de with the second proof, acde, and the edge ab with abe. Summing the probabilities of the three paths would count these two edges twice. Kimmig and de Raedt showed that this reduces to the sum-product problem and proposed an effective solution. The idea is that probability can be computed as a sum if the paths do not share edges. This can be obtained by selecting an edge (or fact), and splitting into the case where the edge is true and the case where the edge is false. The process can be repeated recursively until we run out of facts to split. This idea is the same one used to construct binary decision diagrams (BDDs) as shown before. Fig 4.9 shows a BDD that computes the total probability for the path ae. The total probability is obtained by adding the two cases whether edge ab is true or not. The two cases are clearly disjoint, hence: Pr(βae) = P r(ab)×Prl+Pr(ab)×P rr where Prlis the total probability for the left child βc be, the BDD rooted at be, and P rr is the total probability for the right child βcac, the BDD rooted at ac. Notice that this gives rise to a recursive program but, fortunately, the program can be computed in linear time by using dynamic programming and proceeding bottom-up from the nodes 1 and 0: 66 FCUP Discovering Gene Networks ab be bd de 0 1 ac cd ce Figure 4.9: A BDD that computes the total probability of the path ae. Pr(βc de) = 0.8×1+0.2×0 Pr(βc bd) = 0.7×Pr(βc de)+0.3×0 Pr(βc be) = 0.3×1+0.7×P r(βc bd) Pr(βcce) = 0.3×1+0.7×0 Pr(βc cd) = 0.4×Pr(βc de)+0.6×P r(βcce) Pr(βcac) = 0.5×P r(βc cd)+0.5×P r(βc bd) Pr(βc ab) = 0.2×Pr(βc be)+0.8×P r(βcac) As expected, one can observe that the sub-tree βc de is shared by the path using node band node c. The BDD allows us to count this edge only once. Binary decision diagrams provide a very efficient implementation for probability computation over small and medium graphs. Unfortunately, they do not scale to larger graphs with thousands of nodes. In this case, ProbLog implementations rely on approximated solutions, either Monte Carlo methods or often by approximating the total probability by the probability of the best kproofs [7]. FCUP Discovering Gene Networks 67 4.3.1 Learning Problog Programs Arguably, the most fundamental task in learning ProbLog programs is parameter learning, which aims at obtaining the best values for the fact probabilities given a program structure and examples. There are a number of approaches to this problem. LeProbLog [8] provides a very natural and general algorithm, that is well suited for our task. The goal in parameter learning is to find the best set of probability parameters, Θ, given a set of examples E: maxarg Pr(Θ|E) In our setting, we do not have a prior on the parameters and hence we can assume that all parameters have uniform probability. Applying Bayes theorem, we obtain: Pr(Θ|E)∝P r(E|Θ) LeProbLog maximizes Pr(E|Θ) by using gradient ascent. To compute the gradient of a probability of an example e∈Eon a parameter θi∈Θ we can simply differentiate the expression for total probability: Pr =θj×P rl+ (1 −θj)×P rr If i6=j, then the gradient δPr δθjcan be obtained through the chain rule: δP r δθj =θj×δP rl δθj + (1 −θj)×δP rr δθj If i=j, the expression is slightly more complex: δP r δθj =Prl+θj×δPrl δθj + (1 −θj)×δP rr δθj −Prr We refer the reader to [8] for a complete discussion and detailed derivation assuming the parameters follow a sigmoidal function. In practice, this result means that to compute the gradient over some parameter θiwe simply have to follow a bottom-up dynamic programming algorithm, in a similar style as what we had to do to compute the probability. 68 FCUP Discovering Gene Networks It should be noted that this result is obtained when we have the complete BDD. Often, we use k-proofs, where BDD is only a lower approximation to the correct probability, so we may not converge to local maximum. To address this problem, LeProbLog recomputes the best proofs every ksteps during the gradient ascent process. Chapter 5 Implemented Work If we knew what it was we were doing, it would not be called research, would it? Albert Einstein 69 FCUP Discovering Gene Networks 71 5.1 Representing gene networks with Problog ProbLog requires both logical and probabilistic inferencing. Logical inferencing is implemented though the YAP Prolog engine [93]; probabilistic inferencing relies on the CUDD BDD library for constructing the BDDs [94], and on the SimpleCUDD program to compute probabilities and derivatives. BDDs can grow very quickly and can be quite expensive to create, thus the ProbLog implementers decided to used text files to store all BDDs. In practice, we have noticed that accessing BDDs as files and creating a new instance of SimpleCUDD for processing the derivatives associated with each example is expensive and severely limits performance, making processing of large problems almost impossible. On the other hand, in practice quite often we work with limited-size BDDs. We address this problem by manipulating the BDDs as Prolog terms. A CUDD BDD is translated as a list of Prolog terms, where each term in the list is a node in the original BDD, and is as follows: node(Theta, Left, Right, Value) The Left,Right, and Value logical variables represent, respectively, the node’s value, the value of the left child, and the value of the right child. Figure 5.1 shows a three-node fragment of our graph’s BDD. This BDD can be encoded as the following Prolog term: [node(THETA_BE,PBC,1,PBD), node(THETA_BD,PBD,PDE,0), node(THETA_DE,PDE,1,0)] By using this representation, we can implement probability computation and gradient computation in Prolog as a walk over this list of nodes (ordered bottom-up). More precisely, to compute probabilities, we run the following Prolog code at each node: V alue is θ×Left + (1 −θ)×Right Notice that we do not need to explicitly propagate value upwards. Instead, Prolog unification does the propagation implicitly. 72 FCUP Discovering Gene Networks be bd de 0 1 Figure 5.1: A simple BDD Although Prolog code is much slower at performing arithmetic than the SimpleCUDD’s Ccode, the lower interface cost and avoiding the need to create a process per new BDD makes the Prolog based approach a much faster solution. 5.2 Experimental Methodology A cell responds to environmental changes by detecting molecules that bind to receptors on the surface of the cell and transmits this information to proteins within the cell by activating a cascade of molecular events. These signaling networks are typically studied by measuring molecular events after treatments that stimulate or perturb key elements in the network at the mRNA or gene level. Combining these measurements with phenotypic response enables the study of these treatments on the architecture and function of the underlying signaling networks as well as the relationship between the network behavior and phenotypic reponse [1]. To infer the architecture and function of the underlying signaling network of budding yeast, we decided to focus on the pathways activated by MAPK Hog1 during osmotic stress response. According to previous results [2] this pathway interacts with the general stress (Msn2/Msn4) pathways, so we consider the genes that belong to both pathways our study. We obtained time-series gene expression data from [95] for our analysis. The experiments followed the response of actively growing Saccharomyces cerevisiae subjected to an osmotic shock of 0.7 M NaCl. The dose of salt was chosen by the biologists to provide a robust physiological response but also allow high viability and eventual resumption of cell growth. Three biological replicate samples were collected before and FCUP Discovering Gene Networks 73 after NaCl treatment at 30, 60, 90, 120, and 240 min (measuring the peak transcript changes that occurs at or after 30 min) [96]. We focused our attention on the 270 genes of the Hog1 Msn2/4 pathway from [2] for which we have expression data. The reason we limit our focus to the genes of the Hog1 Msn2/4 pathway is due to the fact that this dataset, unlike that from [2] is much more limited. Furthermore, our main goal is to test the ability of ProbLog to learn pathways based on pairwise correlations or relationships that are computed from gene expression data. We performed three types of analysis to estimate values for the edges, which are then used as training examples in ProbLog: 1. We first compared the three time-series to determine data quality and better understand the general patterns in the data. 2. Computed correlations to determine and quantify the main relationships in the pathways (correlations were mapped from -1 to 1 into 0 to 1). 3. Utilized temporal data to determine temporal relationships in the data (we used a sigmoid function to convert activity levels into 0 to 1). In the first and second analyses we computed correlations between normalized gene expression values. To do so, we use the Pearson product-moment correlation coefficient, also known as r: r=Pn i=1(Xi−¯ X)(Yi−¯ Y) pPn i=1(Xi−¯ X)2pPn i=1(Yi−¯ Y)2(5.1) We first compute correlations between genes in the same experiment, where Xand Y range over the temporal data in order to calculate the correlation (r) between each pair of genes, We experimented with both correlating Xtand Yt(same time-step) and correlating Xtand Yt+1 (next time-step). In our first analysis, we further compared the r-coefficient obtained between all correlations in the three biological replicates. Although all three replicates followed the same methodology, variations in initial conditions can significantly affect gene expression. We study this effect in order to understand the inherent variability existing in the data. In our second analysis, we assume that if a pair of genes has high absolute correlation (a strong positive or negative linear dependence between the two variables) it also has a high probability of being connected in the pathway [97], resulting in two possible connected nodes in the gene interaction network. Correlation of genes based on expression 80 FCUP Discovering Gene Networks Gene ∧ ∨ +− Msn2 22 4 2 Msn4 16 23 Hog1 9 7 4 18 Sko1 9 14 10 Hot1 29 10 2 4 Table 5.3: Experiment 3 (VV): proposed parents per gate. Notice that + represents the positive parent, and − the negative parent. Gene ∧ ∨ +− Msn2 23 3 55 1 Msn4 30 3 52 Hog1 7 1 1 Hot1 1 3 107 Sko1 2 1 Table 5.4: Experiment 3 (LV): proposed parents per gate. Notice that + represents the positive parent, and − the negative parent. 5.3.1 Problog Performance As we discussed before in section 5.1, there was a need to improve Problog code in order to face the performance problem created by the BDDs access and creation of a SimpleCUDD instance. After implementing these changes, we experimentally evaluated1how the new method performed. The following tables and graphs show running times of two versions of Problog (original and the new)- the domain (data) is equal for every experiment (the same data and examples), the results are in the following form - hours : minutes : seconds or minutes : seconds. •No BDDs rebuild Just for looking into Table 5.5 and the graph in Fig. 5.5 one can immediately conclude that the execution times using the two versions are completely different, we compressed the calculation times significantly, as we take from hours to minutes with the new 1All of these tests were executed on a machine with an Intel Core i7 CPU 920 2.67GHz with 12 gigabytes of RAM, running Ubuntu (Release 11.04 - Kernel 2.6.38-16-generic). FCUP Discovering Gene Networks 81 Iterations Original Version New Version V V LV LL V V LV LL 50 1:14:48 1:07:50 1:20:02 1:00 1:07 1:36 100 2:26:54 1:58:05 2:47:07 1:53 2:06 2:52 300 8:35:55 6:41:32 9:30:39 5:25 5:55 8:12 400 12:47:46 9:56:08 13:54:20 7:13 7:53 10:56 Table 5.5: Calculation times differences between the two ProbLog versions with no BDDs rebuild 0" 20" 40" 60" 80" 100" 120" 50"Itera.ons" 100"Itera.ons" 300"Itera.ons" 400"Itera.ons" VV" LV" LL" Figure 5.5: Speedup graph based on Table 5.5 results version. Next we will see the comparison when rebuilding the BDDs after 20 or 45 iterations. •With BDDs rebuild From these results we can also see that even with the extra time needed for the BDDs rebuild, we were able to compress the calculations time significantly, again, we passed form several hours to minutes. This calculation time compressing in both situations (with or without rebuild) allows us to run more complex experiments. 82 FCUP Discovering Gene Networks Iterations Original Version New Version V V LV LL V V LV LL 50 1:26:42 1:06:23 1:38:24 1:20 1:26 2:08 100 2:47:49 2:23:11 3:19:33 2:27 2:41 3:54 300 10:32:23 8:07:54 12:23:01 7:22 7:59 11:45 400 15:16:36 11:13:17 16:56:45 9:56 10:47 15:42 Table 5.6: Calculation times differences between the two ProbLog versions with BDDs rebuild after 20 iterations 0" 10" 20" 30" 40" 50" 60" 70" 80" 90" 100" 50"Itera1ons" 100"Itera1ons" 300"Itera1ons" 400"Itera1ons" VV" LV" LL" Figure 5.6: Speedup graph based on Table 5.6 results Iterations Original Version New Version V V LV LL V V LV LL 50 1:21:42 1:03:40 1:27:37 1:14 1:19 1:57 100 2:55:03 2:19:30 3:14:02 2:15 2:24 3:30 300 9:57:14 7:23:38 11:05:50 6:25 6:52 9:51 400 14:11:28 10:49:00 15:04:43 8:21 9:12 13:10 Table 5.7: Calculation times differences between the two ProbLog versions with BDDs rebuild after 45 iterations FCUP Discovering Gene Networks 83 0" 20" 40" 60" 80" 100" 120" 50"Itera.ons" 100"Itera.ons" 300"Itera.ons" 400"Itera.ons" VV" LV" LL" Figure 5.7: Speedup graph based on Table 5.7 results Chapter 6 Conclusions The outcome of any serious research can only be to make two questions grow where only one grew before. Thorstein Veblen 85 FCUP Discovering Gene Networks 87 6.1 Conclusions Learning regulatory networks from gene expression is a hard problem. Data is noisy, there are often hidden/unmeasured data and relationships between genes are highly complex. We present a statistical relational approach to modeling pathways. Our approach allows us to design a coarser and more fine grained model, based on probabilistic gates. We show that the latter model has predictive performance on the time-series data, and recovers important relationships despite having limited data. We were able to recover some important relationships despite the lack of additional data such as knockout, ChIP-chip. 6.2 Future Work We plan to continue improving the model quality and experiment with new data. Specifically, we would like to experiment with implementing a regression based approach, as it naturally fits our framework. In addition, we would like to experiment with different pathways and with proteomic data. Last, but not least, we would like to investigate how to reduce the number of parameters in the model by exploiting strong correlations between gene expression. Some further goals of this work will include: •Expand our models to also learn from RNASeq data, and to compare RNASeq and microarray data. •Generalize the model to include data from ncRNA expression levels, and to support protein data •Combine probabilistic / logic models with differential equations based models. •Include protein-to-protein interactions, promoters, phosphorylation data in our models. •Proteomics data: how to include it in the model? High-throughput sequencing is known to be an effective approach for transcriptome analysis. This methodology called RNA-seq [99], has been used to analyze unknown 88 FCUP Discovering Gene Networks transcript sequences, estimate gene expression levels and study single nucleotide polymorphisms. Recent RNA-Seq experiments [19] have shown great potential for transcriptome profiling. It is well known that sequencing increases the level of biological detail and also that integrative data analysis is also useful. These are some of the reasons arguing for including RNA-Seq data on the proposed work. Non-coding RNA (ncRNA) is the next ’frontier’ for the biological community, this has not been crossed yet, and it continues to present a lot of questions and doubts [11], [12], [13]. If the model is generalized to include data also from ncRNA, it would be a great step forward, and could lead to the discovery of connections between ncRNA’s and other genes, and also the discovery of what they regulate or even co-regulate, which would be a significant contribution to the scientific community regarding the common comprehension and also advance on disease treatments. Proteomic data is a potentially rich, but maybe under-exploited, data source for genome annotation. Peptide identifications from tandem mass spectrometry provide prima facie evidence for gene predictions and can also discriminate over a set of candidate gene models, making this type of data a good candidate for this work model. Gene regulatory networks (GRN) models are difficult to deduce just by using experimental techniques, so, computational and mathematical methods are indispensable. As GRNs are nonlinear, nonlinear differential equation models can model much more complex GRN behavior. The identification of the nonlinear differential equation in these type of models is computationally more intensive and may require more data. Despite this, the range of nonlinear behaviors exhibited by GRNs can be comprehensively realized with nonlinear differential equations. Appendix A Code ################################################ # Delta as a function of 2 other Delta # # Delta(A,t) AND Delta(B,t) -> Delta(C,t+1) # # Delta(A,t) OR Delta(B,t) -> Delta(C,t+1) # # Delta(A,t) -> Delta(C,t+1) # # Delta(A,t) OR_NOT Delta(B,t) -> Delta(C,t+1) # ################################################ next(30,60). next(60,90). next(90,120). next(120,240). delta_ge(E,T1,Z) :- and(X,Y,Z), \+ or(_,_,Z), next(T0,T1), delta(E,T0,X), delta(E,T0,Y). delta_ge(E,T1,Z) :- or(X,_Y,Z), next(T0,T1), delta(E,T0,X). delta_ge(E,T1,Z) :- 89 96 FCUP Discovering Gene Networks next(T0,T1), ge(E,T0,X). delta_not_ge(E,T1,Z) :- and(X,_Y,Z), next(T0,T1), not_ge(E,T0,X). delta_not_ge(E,T1,Z) :- and(_X,Y,Z), next(T0,T1), not_ge(E,T0,Y). delta_not_ge(E,T1,Z) :- or(X,Y,Z), next(T0,T1), not_ge(E,T0,X), not_ge(E,T0,Y). delta_not_ge(E,T1,Z) :- or_not(X,Y,Z), next(T0,T1), not_ge(E,T0,X), ge(E,T0,Y). delta_not_ge(E,T1,Z) :- single(X,Z), next(T0,T1), not_ge(E,T0,X). fal_and(Z) :- and(X,Y,Z), and(X1,Y1,Z), ( X1 \= X ; Y1 \= Y). fal_and(Z) :- and(_,_,Z), or(_,_,Z). fal_and(Z) :- and(_,_,Z), or_not(_,_,Z). fal_and(Z) :- and(_,_,Z), single(_,Z). fal_or(Z) :- or(X,Y,Z), or(X1,Y1,Z), ( X1 \= X ; Y1 \= Y). fal_or(Z) :- or(_,_,Z), and(_,_,Z). fal_or(Z) :- or(_,_,Z), or_not(_,_,Z). FCUP Discovering Gene Networks 97 fal_or(Z) :- or(_,_,Z), single(_,Z). fal_or_not(Z) :- or_not(X,Y,Z), or_not(X1,Y1,Z), ( X1 \= X ; Y1 \= Y). fal_or_not(Z) :- or_not(_,_,Z), and(_,_,Z). fal_or_not(Z) :- or_not(_,_,Z), or(_,_,Z). fal_or_not(Z) :- or_not(_,_,Z), single(_,Z). fal_single(Z) :- single(X1,Z), single(X,Z), X \= X1. fal_single(Z) :- single(_,Z), and(_,_,Z). fal_single(Z) :- single(_,Z), or(_,_,Z). fal_single(Z) :- single(_,Z), or_not(_,_,Z). References [1] Melody Morris, Julio Saez-Rodriguez, Peter Sorger, and Douglas Lauffenburger. Logic-based models for the analysis of cell signaling networks. Biochemistry, 49:3216–3224, 2010. [2] Andrew Capaldi, Tommy Kaplan, Ying Liu, Naomi Habib, Aviv Regev, Nir Friedman, and Erin O’Shea. Structure and function of a transcriptional network activated by the mapk hog1. Nature Genetics, 40:1300–1306, 2008. [3] L. Glass and S.A. Kauffman. A logical analysis of continuous, non-linear biochemical control networks. Journal of Theoretical Biology, 39:103–129, 1973. [4] R. Thomas. Boolean formalization of genetic control circuits. Journal of Theoretical Biology, 42:563–585, 1973. [5] Nir Friedman, Michal Linial, Iftach Nachman, and Dana Pe’er. Using Bayesian networks to analyze expression data. Journal of Computational Biology, 7(3/4):601–620, 2000. [6] Luc De Raedt, Angelika Kimmig, and Hannu Toivonen. Problog: A probabilistic prolog and its application in link discovery. In Manuela M. Veloso, editor, IJCAI 2007, Proceedings of the 20th International Joint Conference on Artificial Intelligence, Hyderabad, India, January 6-12, 2007, pages 2462–2467, 2007. [7] Angelika Kimmig, V´ıtor Santos Costa, Ricardo Rocha, Bart Demoen, and Luc De Raedt. On the Implementation of the Probabilistic Logic Programming Language ProbLog. Theory and Practice of Logic Programming Systems, 11:235–262, 2011. [8] B. Gutmann, A. Kimmig, K. Kersting, and L. De Raedt. Parameter learning in probabilistic databases: A least squares approach. In ECML/PKDD–08, volume LNCS 5211, pages 473–488, Antwerp, Belgium, September 15–19 2008. Springer. 99 100 FCUP Discovering Gene Networks [9] J. Gebert, N. Radde, and G.-W. Weber. Modeling gene regulatory networks with piecewise linear differential equations. European Journal of Operational Research, 181(3):1148 – 1165, 2007. [10] Anirban Ghosh and Manju Bansal. A glossary of dna structures from a to z. Acta Crystallographica Section D, 59(4):620–626, Apr 2003. [11] Jeremy Berg. Biochemistry. W.H. Freeman, New York, 2002. [12] Wolfram Saenger. Principles of nucleic acid structure. Springer-Verlag, New York, 1984. [13] John Butler. Forensic DNA typing : biology and technology behind STR markers. Academic Press, San Diego, 2001. [14] J. D. Watson and F. H. C. Crick. Molecular Structure of Nucleic Acids: A Structure for Deoxyribose Nucleic Acid. nat, 171:737–738, April 1953. [15] Peter Yakovchuk, Ekaterina Protozanova, and Maxim D. Frank-Kamenetskii. Base-stacking and base-pairing contributions into thermal stability of the dna double helix. Nucleic Acids Research, 34(2):564–574. [16] H. Clausenschaumann. Mechanical Stability of Single DNA Molecules. Biophysical Journal, 78:1997–2007, April 2000. [17] Qin Yu and Casey D. Morrow. Identification of critical elements in the trna acceptor stem and tc loop necessary for human immunodeficiency virus type 1 infectivity. Journal of Virology, 75(10):4902–4906, 2001. [18] Miguel Salazar, Oleg Y. Fedoroff, Julie M. Miller, N. Susan Ribeiro, and Brian R. Reid. The dna strand in dna.cntdot.rna hybrid duplexes is neither b-form nor a-form in solution. Biochemistry, 32(16):4207–4215, 1993. PMID: 7682844. [19] Satu Mikkola, Eeva Stenman, Kirsi Nurmi, Esmail Yousefi-Salakdeh, Roger Stromberg, and Harri Lonnberg. The mechanism of the metal ion promoted cleavage of rna phosphodiester bonds involves a general acid catalysis by the metal aquo ion on the departure of the leaving group. J. Chem. Soc., Perkin Trans. 2, pages 1619–1626, 1999. [20] Z. Tan. Salt Dependence of Nucleic Acid Hairpin Stability. Biophysical Journal, 95:738–752, July 2008. [21] Jonathan Pevsner. Bioinformatics and functional genomics. Wiley, 2003. FCUP Discovering Gene Networks 101 [22] Sheetal A Mitra, Anirban P Mitra, and Timothy J Triche. A central role for long non-coding rna in cancer. Frontiers in Genetics, 3(00017), 2012. [23] Michael D Horwich, Chengjian Li, Christian Matranga, Vasily Vagin, Gwen Farley, Peng Wang, and Phillip D Zamore. The drosophila rna methyltransferase, dmhen1, modifies germline pirnas and single-stranded sirnas in risc. Curr Biol, 17(14):1265–72, 2007. [24] Mads Kaern, William J. Blake, and J. J. Collins. The engineering of gene regulatory networks. Annual review of biomedical engineering, 5:179–206, January 2003. [25] Gemma Bell´ı, Eloi Gar´ı, Lidia Piedrafita, Mart´ı Aldea, and Enrique Herrero. An activator/repressor dual system allows tight tetracycline-regulated gene expression in budding yeast. Nucleic Acids Research, 26(4):942–947, 1998. [26] M Gossen, S Freundlieb, G Bender, G Muller, W Hillen, and H Bujard. Transcriptional activation by tetracyclines in mammalian cells. Science, 268(5218):1766– 1769, 1995. [27] S. Nagahashi, H. Nakayama, K. Hamada, H. Yang, M. Arisawa, and K. Kitada. Regulation by tetracycline of gene expression in saccharomyces cerevisiae. Molecular and General Genetics MGG, 255:372–375, 1997. 10.1007/s004380050508. [28] Barbara Ulmasov, John Capone, and William Folk*. Regulated expression of plant trna genes by the prokaryotic tet and lac repressors. Plant Molecular Biology, 35:417–424, 1997. 10.1023/A:1005819007549. [29] Damien Eveillard, Delphine Ropers, Hidde Jong, Christiane Branlant, and Alexander Bockmayr. Multiscale modeling of alternative splicing regulation. In Corrado Priami, editor, Computational Methods in Systems Biology, volume 2602 of Lecture Notes in Computer Science, pages 75–87. Springer Berlin Heidelberg, 2003. [30] A. M Walczak and G. Tkaˇcik. Information transmission in genetic regulatory networks: a review. ArXiv e-prints, jan 2011. [31] M. Schena. DNA microarrays: A practical approach, volume 205 of Practical Approach Series. Oxford Univ. Press., Oxford, 1999. [32] Hidde De Jong. Modeling and simulation of genetic regulatory systems: A literature review. Journal of Computational Biology, 9:67–103, 2002. 102 FCUP Discovering Gene Networks [33] U. Alon, N. Barkai, D. A. Notterman, K. Gish, S. Ybarra, D. Mack, and A. J. Levine. Broad patterns of gene expression revealed by clustering analysis of tumor and normal colon tissues probed by oligonucleotide arrays. Proceedings of the National Academy of Sciences, 96(12):6745–6750, 1999. [34] Raymond J. Cho, Michael J. Campbell, Elizabeth A. Winzeler, Lars Steinmetz, Andrew Conway, Lisa Wodicka, Tyra G. Wolfsberg, Andrei E. Gabrielian, David Landsman, David J. Lockhart, and Ronald W. Davis. A genome-wide transcriptional analysis of the mitotic cell cycle. Molecular Cell, 2(1):65 – 73, 1998. [35] George S. Michaels, Daniel B. Carr, Stefanie Fuhrman, Xiling Wen, and Roland Somogyi. Cluster analysis and data visualization of largescale gene expression data. 1998. [36] X. Wen, S. Fuhrman, G. S. Michaels, D. B. Carr, S. Smith, J. L. Barker, and R. Somogyi. Large-scale temporal gene expression mapping of central nervous system development. PNAS, 95(1):334–339, 1998. [37] Paul T. Spellman, Gavin Sherlock, Michael Q. Zhang, Vishwanath R. Iyer, Kirk Anders, Michael B. Eisen, Patrick O. Brown, David Botstein, and Bruce Futcher. Comprehensive identification of cell cycle–regulated genes of the yeast saccharomyces cerevisiae by microarray hybridization. Molecular Biology of the Cell, 9(12):3273–3297, 1998. [38] Michiel De Hoon, Seiya Imoto, and Satoru Miyano. Inferring gene regulatory networks from time-ordered gene expression data of bacillus subtilis using differential equations. In Pac. Symp. Biocomput, pages 17–28, 2003. [39] Ting Chen, Hongyu L. He, and George M. Church. Modeling gene expression with differential equations. In Pac. Symp. Biocomput, pages 29–40, 1999. [40] Michiel de Hoon, Seiya Imoto, and Satoru Miyano. Inferring gene regulatory networks from time-ordered gene expression data using differential equations. In Steffen Lange, Ken Satoh, and Carl Smith, editors, Discovery Science, volume 2534 of Lecture Notes in Computer Science, pages 283–288. Springer Berlin / Heidelberg, 2002. [41] Erina Sakamoto and Hitoshi Iba. Inferring a system of differential equations for a gene regulatory network by using genetic programming. In Proc. Congress on Evolutionary Computation, pages 720–726. IEEE Press, 2001. FCUP Discovering Gene Networks 103 [42] Thomas G. Kurtz. Solutions of ordinary differential equations as limits of pure jump markov processes. Journal of Applied Probability, 7(1):pp. 49–58, 1970. [43] T. G. Kurtz. Limit theorems for sequences of jump markov processes approximating ordinary differential processes. Journal of Applied Probability, 8(2):pp. 344–356, 1971. [44] Tatsuya Akutsu, Satoru Miyano, and Satoru Kuhara. Inferring qualitative relations in genetic networks and metabolic pathways. Bioinformatics, pages 727–734, 2000. [45] Michael A. Savageau and Eberhard O. Voit. Recasting nonlinear differential equations as s-systems: a canonical nonlinear form. Mathematical Biosciences, 87(1):83 – 115, 1987. [46] Eberhard O. Voit. Symmetries of s-systems. Mathematical Biosciences, 109(1):19 – 37, 1992. [47] Douglas H. Irvine and Michael A. Savageau. Efficient solution of nonlinear ordinary differential equations expressed in s-system canonical form. SIAM Journal on Numerical Analysis, 27(3):pp. 704–735, 1990. [48] Diego di Bernardo, Michael J Thompson, Timothy S Gardner, Sarah E Chobot, Erin L Eastwood, Andrew P Wojtovich, Sean J Elliott, Scott E Schaus, and James J Collins. Chemogenomic profiling on a genome-wide scale using reverseengineered gene networks. Nat Biotechnol, 23(3):377–83, 2005. [49] T. S. Timothy, D. d. Diego, D. David, and J. J. James. Inferring genetic networks and identifying compound mode of action via expression profiling. Science, 301(5629):102–105, July 2003. [50] Mukesh Bansal, Giusy Della Gatta, and Diego di Bernardo. Inference of gene regulatory networks and compound mode of action from time course gene expression profiles. Bioinformatics, 22(7):815–822, 2006. [51] E. P. van Someren, B. L. T. Vaes, W. T. Steegenga, A. M. Sijbers, K. J. Dechering, and M. J. T. Reinders. Least absolute regression network analysis of the murine osteoblast differentiation network. Bioinformatics, 22(4):477–484, 2006. [52] Jesper Tegn´er, M. K. Stephen Yeung, Jeff Hasty, and James J. Collins. Reverse engineering gene networks: Integrating genetic perturbations with dynamical 104 FCUP Discovering Gene Networks modeling. Proceedings of the National Academy of Sciences, 100(10):5944–5949, 2003. [53] Richard Bonneau, David Reiss, Paul Shannon, Marc Facciotti, Leroy Hood, Nitin Baliga, and Vesteinn Thorsson. The inferelator: an algorithm for learning parsimonious regulatory networks from systems-biology data sets de novo. Genome Biology, 7(5):R36, 2006. [54] Patrik D’haeseleer, X. Wen, Stefanie Fuhrman, and Roland Somogyi. Linear modeling of mrna expression levels during cns development and injury. In Pacific Symposium on Biocomputing’99, pages 41–52, 1999. [55] T. Akutsu, S. Kuhara, O. Maruyama, and S. Miyano. Identification of gene regulatory networks by strategic gene disruptions and gene overexpressions. In Proc. the 9th Annual ACM-SIAM Symposium on Discrete Algorithms, pages 695– 702, 1998. [56] T.E. Ideker, V. Thorsson, and R.M. Karp. Discovery of regulatory interactions through perturbation: Inference and experimental design. In Pacific Symposium on Biocomputing, pages 302–313, 2000. [57] Ilya Shmulevich, Edward R. Dougherty, Seungchan Kim, and Wei Zhang. Probabilistic boolean networks: a rule-based uncertainty model for gene regulatory networks. Bioinformatics, 18(2):261–274, 2002. [58] Gilles Bernot, Jean-Paul Comet, Adrien Richard, and Janine Guespin. Application of formal methods to biological regulatory networks: extending thomas’ asynchronous logical approach with temporal logic. Journal of Theoretical Biology, 229(3):339 – 347, 2004. [59] Gr´egory Batt, Delphine Ropers, Hidde De Jong, Johannes Geiselmann, Radu Mateescu, and Dominique Schneider. Validation of qualitative models of genetic regulatory networks by model checking: Analysis of the nutritional stress response in escherichia coli. Bioinformatics, 21:19–28, 2005. [60] Hidde de Jong, Johannes Geiselmann, C´eline Hernandez, and Michel Page. Genetic network analyzer: qualitative simulation of genetic regulatory networks. Bioinformatics, 19(3):336–344, 2003. [61] R. Thomas and M. Kaufman. Multistationarity, the basis of cell differentiation and memory. ii. logical analysis of regulatory networks in terms of feedback FCUP Discovering Gene Networks 105 circuits. Chaos: An Interdisciplinary Journal of Nonlinear Science, 11(1):180– 195, 2001. [62] M. Kaufman, J. Urbain, and R. Thomas. Towards a logical analysis of the immune response. Journal of Theoretical Biology, 114(4):527 – 561, 1985. [63] Ren´e THOMAS, Anne-Marie GATHOYE, and Lucie LAMBERT. A complex control circuit. European Journal of Biochemistry, 71(1):211–227, 1976. [64] R. Thomas. Logical analysis of systems comprising feedback loops. Journal of Theoretical Biology, 73(4):631 – 656, 1978. [65] Ren´e Thomas, Denis Thieffry, and Marcelle Kaufman. Dynamical behaviour of biological regulatory networks—i. biological role of feedback loops and practical use of the concept of the loop-characteristic state. Bulletin of Mathematical Biology, 57:247–276, 1995. 10.1007/BF02460618. [66] J. Demongeot, D. Benaouda, and C. J´ez´equel. “dynamical confinement” in neural networks and cell cycle. Chaos: An Interdisciplinary Journal of Nonlinear Science, 5(1):167–173, 1995. [67] Jacques Demongeot, Marcelle Kaufman, and Ren´e Thomas. Positive feedback circuits and memory. Comptes Rendus de l’Acad´emie des Sciences - Series III - Sciences de la Vie, 323(1):69 – 79, 2000. [68] Jacques Demongeot, Julio Aracena, Florence Thuderoz, Thierry-Pascal Baum, and Olivier Cohen. Genetic regulation networks: circuits, regulons and attractors. Comptes Rendus Biologies, 326(2):171 – 188, 2003. [69] Steven Eker, Merrill Knapp, Keith Laderoute, Patrick Lincoln, , and Carolyn Talcott. Pathway logic: Executable models of biological networks. In Fourth International Workshop on Rewriting Logic and Its Applications (WRLA’2002, volume 71 of Electronic Notes in Theoretical Computer Science. Elsevier, 2002. [70] Hiroaki Kitano, Akira Funahashi, Yukiko Matsuoka, and Kanae Oda. Using process diagrams for the graphical representation of biological networks. Nat Biotech, 23(8):961 – 966, 2005. [71] Hiroaki Kitano. A graphical notation for biochemical networks. BIOSILICO, 1(5):169 – 176, 2003.