Paquete de R para análisis de redes epistáticas

2 ESCUELA TÉCNICA SUPERIOR DE INGENIERÍA INFORMÁTICA Grado en Ingeniería de la Salud Paquete de R para análisis de redes epistáticas R Package for the analysis of epistatic networks Realizado por Samuel F. Benavides García Tutorizado por Oswaldo Trelles Salazar Óscar Torreño Tirado Departamento Arquitectura de computadores UNIVERSIDAD DE MÁLAGA Málaga, septiembre de 2016 Fecha defensa: El Secretario del Tribunal 3 4 Resumen Las mutaciones o cambios en la dotación genética de los organismos, son uno de los mecanismos básicos de la evolución de las especies. Estas pueden encontrarse en zonas codificantes o zonas no codificantes del genoma, siendo las primeras -las encontradas en zonas codificanteslas que suscitan más interés ya que en muchos casos pueden tener efectos negativos asociados con la prevalencia de enfermedades en el organismo. Entre estos efectos negativos, las mutaciones también pueden estar relacionadas con la aparición de enfermedades asumiendo que dichos cambios aumentan la propensión del individuo a padecer una enfermedad. Hoy en día se acepta que muchas de las enfermedades con origen en mutaciones son causadas por mecanismos epistáticos, esto es una interacción entre varias mutaciones que tienen efecto en su conjunto sobre una enfermedad. Estas relaciones epistáticas o de alto orden presentan dificultades para ser estudiadas debido a la complejidad computacional que presentan. Por lo tanto, es necesario diseñar nuevos métodos software de estudio de relaciones Genotipo-Fenotipo que permitan un análisis exhaustivo epistático, en este caso por parejas, complementando el análisis con teoría de grafos. Este trabajo pretende contribuir a la investigación actual identificando relaciones Genotipo-Fenotipo, especialmente las que requieren análisis epistático, creando un paquete de software para R, uno de los lenguajes de programación más usados en biología computacional y en biomedicina, ofreciendo una interfaz de uso simplificada para completar el análisis de redes epistáticas creando un grafo de Polimorfismos de Nucleótidos Únicos o SNPs según sus siglas en inglés. A partir de dichos grafos el usuario puede obtener un grafo de genes y permitiendo al usuario obtener los genes o SNPs con mayor centralidad. La herramienta a desarrollar también permite al usuario almacenar sus resultados intermedios y finales para poderlos procesar posteriormente con otros programas externos. La herramienta a desarrollar permite su gestión con programas externos como Galaxy, un WMS (Workflow Management System) que permite la gestión de flujos de trabajo, y Cytoscape, un software especializado en la visualización y gestión de grafos. Ambas herramientas son comúnmente utilizada por el público objetivo. El resultado de este trabajo pretende hacer posible la mejora en los estudios de esta área de investigación, permitiendo a más usuarios con menor conocimiento tecnológico realizar estudios GWAS expandiendo su uso y favoreciendo el incremento de resultados y conocimiento. Palabras clave: epistasis, grafos, R, GWAS, Cytoscape, genotipo-fenotipo, Galaxy 5 Abstract Mutations or changes in the genetic endowment are the engine of evolution. They can be found in coding or non-coding regions of the genome, those found on coding regions are usually the most interesting ones as in many cases they can lead to changes in the behaviour of the effects on the organism. Between negative effects, mutations have been traditionally related to diseases assuming that these changes favor a certain disease. Nevertheless nowadays it is well known that many diseases are caused by epistasis which means that interactions between several mutations cause the disease. Due to the complexity it is problematic to study these epistatic relationships or high order relationships. As a result, new methods for studying genotype-phenotype need to be designed. These methods use, in this specific case, pairwise exhaustive epistatic analysis complement it with graph theory. The main purpose of this study is to contribute to the current state of the art with a software capable to identify genotype-phenotype relationships, especially those that were described previously. A software package written in R, which is one of the most used programing languages among software developers in biology and biomedicine, is created and offers a simple interface to complete an epistatic network analysis. The package creates a graph of SNPs turning it into a genes allowing the user to choose the top genes or SNP and then gives the possibility to open the graph in an external software. The aim of the tool is to be able to allow its management with external programmes like they are Galaxy, for creation and management of workflows, and Cytoscape, a software for visualization and management of graphs. Both tools are commonly used by the target population. The result of this project intends to provide all necessary tools to allow more users to perform GWAS expanding the use of it and favouring the increment of results and knowledge. Keywords: epistasis, graph, R, GWAS, Cytoscape, geno-pheno, Galaxy 6 Index Section Page 1. Introduction 9 1.1. Motivation 9 1.2. Objectives 10 1.3. Methodology 11 2. State of the art 12 2.1. General review 12 2.2. Postgwas 14 2.3. Postgwas 16 2.4. Other software 20 3. Analysis and design 24 4. Development 30 4.1. Introduction 30 4.2. createGraphFromTable 30 4.3. fillInfo 30 4.4. replaceSNPwithGene 31 4.5. getTopDegree 32 4.6. openGraph 33 4.7. workflow implementation in Galaxy 35 5. Use case and validation 37 5.1. Dataset 37 5.2. Replication 37 5.3. Use case 41 6.1 Conclusions 45 6.2 Conclusiones (versión en español) 46 7. References 47 7 8 1. Introduction 1.1. Motivation As consequence of new sequencing technological breakdown, commonly referred as ‘next generation sequencing’[1] the price of sequencing nucleotide chains has decreased leading to an increase of the popularity of these technologies on the sequencing market. The price decreasement is related also to the fact of the frequent usage of these methods. Nowadays one of the main applications of sequencing data (genes, genomes) is biomedical research, and it is also starting to be used for clinical purposes. Furthermore in last years, as a result of the decrease of its price, the technology is starting to be used often in hospitals in order to identify pathologies. Figure 1 Cost to sequence a human genome expressed in thousand of USD 1 Once the genome is sequenced or alternative methods such as SNP-microarray [2] are used to detect mutations, these mutations are taken into account to diagnose a patient or predict a disease. In order to make that possible it is necessary to perform a GWAS (Genome Wide Association Study) which is used principally to relate mutations to diseases, in other words to infer relationships genotype-phenotype. This type of studies has identified several mutations related to diseases such as asthma [3], breast cancer [4] or hypertension [5] among many others. Traditionally a single isolated SNP (Single Nucleotide Polymorphism) [6] has been associated to a phenotypical expression. However this one-to-one correlations only cover a small part of the full spectrum of diseases, as many of the phenotype instances may be caused by epistasis, which means interactions between two or 1 https://www.genome.gov/27541954/dna-sequencing-costs-data/ 9 Considering all previously mentioned it seems necessary that whole code will be written in R, even though R can run compiled C code. Using R opens the possibility of using external libraries such as ‘igraph  ’ to build and manage part of the graphs to be created in the package, not only this but it makes the code fully open to the user. Another point mentioned before is the name of the package. It may seem as an unimportant decision but it is the first contact with the user and as the user may not be looking for a specific package the name of the package should be descriptive. Many packages incorporate the letter ‘R’ in upper case into their names relating the package with the R technology. Considering this a good descriptive name would be ‘epiR  ’, however, such naming convention is not currently in use neither it was used in any previous existing package in the CRAN or in Bioconductor as mentioned in the CRAN policy. 2.3. Postgwas This R package offers several different functions for the interpretation of the results obtained by GWAS, including a method for identifying a gene giving a SNP and a flaking distance or linkage disequilibrium. It is also aimed to create a graph of SNPs to obtain a gene enrichment network. A more detailed view of the methods contained in the package can be seen in Figure 4. 16 Figure 4 Diagram of functions in the postgwas package [16] The package works focused in a dataframe which is the main enviroment, this dataframehas a function to create it “setPostgwasBuffer  ”, a function to get a clone of the buffer “getPostgwasBuffer  ” and a function to clear it “clearPostgwasBuffer  ”. From this buffer all functions are applied, “getGenotypes  ” is a function developped to be used mainly by other functions, this retrieves when possible the genotype of a given set of SNPs, “snp2gene  ” is one of these functions and adds to the given dataframe a field named “gene” with the gene where the SNP is located, two specifications are created, one is “snp2gene.prox  ”, which is the one arousing more interest, this functions adds a flaking distance where a SNP can be located to be placed on that specific gene; other specification is “snp2gene.ld  ” which takes into account the linkage disequilibrium to assign a gene to a SNP. 17 There are two plotting functions in this package, one of them is “regionalPlot  ”, it creates regional plots for a given window size around a given list of SNPs. These regional plots contain one or more p-value graphs, gene information and LD between SNPs from the GWAS. Other plotting function is “manhattanPlot  ” it plots a manhattan plot, with SNPs exceeding a user-defined threshold being highlighted and annotated with closest genes. Four other important functions are “gene2p  ”, it calculates aggregate p-values for genes assigned to multiple SNPs, by taking into account the dependency between SNPs. Another important function is “getInteractions  ” which uses biomart to determine pathway annotations for each gene, then the common pathways between genes form interactions. A function to visualize previous network is “gwas2network  ” which displays the given network and highlights all given genes and their connecting edges using transparency effects. The last function is “gwasGOenrich  ” provides an interface for GO term enrichment, it returns a dataframe with id, terms and P-value columns containing identifier, description and enrichment p-value of all GO terms tested. This might be one of the strongest references as it carries most of the functions we want to implement. The first use of this package will be essential in order to measure its performance and evaluate the user’s experience and will be used to compare our developed software. Firstly the package was removed in December 2013 from the CRAN for no response 3 nor update was received after a request to clean up warnings. Therefore it must be downloaded externally in a zipped file and installed manually instead of doing it directly from the CRAN repository. Normally when such a situation occurs it is because there is no package update adapting it to the rules of the CRAN. Surprisingly no further work has been done with the package, as pointed out in the CRAN. The files are available in github, where you can encounter a misleading information in the help section which says that the package can be downloaded from the CRAN, being this version same as the last version available from the CRAN’s archive this version is 1.11. There is a summary function in the package named ‘postgwas  ’, which has among its documentation a usage example. When running it not only there was no result but some errors were raised. Several attempts to run it were made and during none of them there was a positive result. 3 https://cran.r-project.org/web/packages/postgwas/index.html 18 After carefully observing the messages obtained we can get to the conclusion that it is a problem using an external service as the error is raised every time a external service was invoked, it is discarded a possible mistake on the input data as it was used the sample data provided by the package. Therefore none of the functions in this package could be used in future implementations as the lack of update together with the impossibility to run the code and its removal from the CRAN make this source code an unstable basis for further developed functions. 2.4. Other software Although Postgwas is one of the packages with higher similar functionality than our package we want to create, there are many other softwares for their use in GWAS studies. There exist a webpage which sort many kind of software for their use in omics science named “Omic Tools”[30], it offers as well tools for GWAS but these tools are mainly focused to single association studies. There are many methods for detecting epistasis in a wide range of interfaces, I will describe only the few most used although there are many others, PLINK is one of the most well-known method, it takes a regressive approach and it’s interface is the command line, it is widely used and has a wide range of input formats, it supports bit-level based parallelism; BOOST [31] takes a regressive approach as well and has a command line interface, it performs an exhaustive search but has statistical power limitations; BEAM [32] takes a bayesian approach and also has a command line interface, it incorporates posterior information but does not support parallelism natively; MB-MDR [33] takes a data mining approach and is used through the command line, it is widely used and updated regularly but does not offer native support for Windows OS; FaST-LMM [34] takes a linear mixed model approach and is used also through the command line, it performs an exhaustive search but is computationally expensive. A deeper analysis of the previous software will be done in this paragraph. PLINK includes a regression-based epistasis test. However it takes quite a lot of time to perform a test. A faster implementation was developed later on including native parallelism support [35], while this method is much faster it is only applicable to quantitative traits and cannot be used in case/control studies. BOOST is able to evaluate a small set of SNPs in a standard desktop computer. Bayesian methods are one of the alternatives to exhaustive methods like previous, this not require all 19 interactions to be tested, an example of this is BEAM which divides SNP in three different categories based on inner probabilities. Another way for detection of epistasis is machine learning and data mining MB-MDR is an example of this, it categorizes genotypes into high-risk and low-risk reducing data to a single dimension, however as it is an exhaustive method as well, it is generally applied to studies where the dataset has previously been reduced. Epistatic analysis present computational issues and as methods become more complex and data sets increase in size, computational Normally after the epistasis detection there are used specific software to obtain results, this software input normally is generic, that means that is not specially designed for epistatic analysis. There is no similar package for management of epistatic networks after the epistatic analysis, many GWAS software found is based on single association studies, and those found which work with epistasis performs only the epistatic detection. In the area of workflows, although there are many options, one interesting WMS is Galaxy as it is a web based tool with a great community on the area of the life science as can be seen in Table 2. Being web based makes easy to create and share workflows and tools by creating an instance of a workflow in a server so anyone can both download and access the designed tool, this makes Galaxy a perfect tool for the reproducibility of any study following a workflow implementation. Galaxy runs completely on the server therefore there is no specific requirement from the user’s computer to run any program. It is also a easy to use software allowing non-computational users to run a program with no need of dealing with the command line. 20 Date Posts Threads Questions Answers Users 2014/07/09 8223 2978 2959 2908 2517 2014/10/13 9230 3252 3220 3252 3555 2015/01/23 10298 3547 3505 3646 4700 2015/07/24 13275 4379 4303 4616 6342 2016/01/29 15742 5035 4958 5469 7909 2016/08/15 18823 5935 5826 6463 9731 Table 2. Statistics from the past two years of the Galaxy Biostar forum showing users, posts, threads, questions and answer. 4 There are three ways to use galaxy, first one is trough the official usegalaxy.org portal which offers plenty of tools thought for bioinformatics bust doesn’t allow the user to add it’s own modules, other way to use Galaxy is accessing to any of it’s many servers which have specialized modules depending on their expertise, the 5 last way how to use Galaxy is installing a local instance, where you can add all necessary modules, this offers the opportunity to publish it through a public server. The main view from a Galaxy page can be seen in Figure 5 and it shows the groups of modules on the left side where the epiR group with all necessary modules will be. 4 https://wiki.galaxyproject.org/GalaxyProject/Statistics 5 https://wiki.galaxyproject.org/PublicGalaxyServers 21 Figure 5 Main view of the galaxy home, on the left side there are the tools, on the right a history with all datasets and results obtained and in the upper part there is the menu. On the area of the enrichment analysis David [36] which is the version developed by the NCBI or WebGestalt [37] developed by the zhang lab, both are interesting software for this purpose, and they are web-based but the second one offers more options at the time of performing the enrichment analysis. Considering the interface simplicity, as seen in Figure 6, and the variety of analysis WebGestalt is recommended. Figure 6 Main Webgestalt interface offering all possible parameters to the user 22 3. Analysis and Design Taking into account the aim of this project the functionality that we want to implement must be chosen carefully, obtaining the requirements for our package. The package has to offer a quick connection between the results and a software which is able to visualize the result graph, either a SNPs one or a genes one. Thus the first idea is to connect the package to Cytoscape, is a de-facto standard network visualization application in life science field. One option could be that a function directly launches Cytoscape with the required graph (this is the simplest option) and the other option would be to save the graph in a Cytoscape readable format and plot the top clusters in an image so the user can see the biggest clusters. While the first option seems to be the best, it requires to add Cytoscape to the system path, which is something that could cause some troubles to the users, on the other hand the second option offers low quality as plotted images cannot be zoomed. The final decision was to invoke Cytoscape directly from the function, loading the required graph automatically (as seen in Figure 7) and saving the graph in .gml format which is a format that can be imported by Cytoscape. This leads to a conflict with the CRAN repository. The conflict was solved by implementing the package in two versions. The first one is meant to be published in the CRAN, which will plot an image instead, and another one is meant to be manually installed, this last version would launch Cytoscape. Figure 7 Cytoscape main view after loading a graph, in the lower left corner there is a miniature of the graph and a toolbar is placed on the uppest side. 23 Considering the big computational effort to perform an epistatic analysis it seems to be obvious that it has to be performed in a HPC infrastructure as computing time could be in the range of years if performed on an average user PC. Therefore the input for the developed R package must be the result from the epistatic analysis, typically provided to biologists and clinicians in plain text edge list format as seen in Figure 7. This input has also the advantage that as it is an edge list in plain text it is easily editable by experts applying their own knowledge (e.g., known SNPs which are not related to the disease, removal of certain relations as known from low information, etc). Another advantage is that the user can create his/her own dataset based on his/her experience or simply the user may want to perform an experiment without any previous need of obtaining biological samples. A pipeline design will be extremely suitable for this project as it is based on a workflow, therefore the functions must share inputs and outputs and for this a WMS is a good option. At this point we have a basic scheme about the workflow to be developed. Now we present a more detailed idea about the components: As input we have an edge list which basically consist on the edges of a graph set in pairs as seen in figure 7. The edge list is made by pairs of RSIDs that correspond with the SNP’s identificators. Second, we create a graph of SNPs and as well we will need a tool to get a gene from a given rsid and a flanking distance. As other previously existing functions like snp2gene  from the package postgwas are currently out of service the gene graph will be generated from the SNP graph. The input format is a plain text with two columns separated by tabulators where any of the columns has to contain a valid rsid as shown in figure 8. This file may not have a header. SNP1 SNP2 rs1487013 rs1593270 rs12532565 rs17105765 rs12356602 rs11842523 rs129128 rs11842523 Figure 8 Example of an input edge list 24 To evaluate and visualize the graphs we need two tools, one to get the most relevant nodes from a graph according to a selected metric and another to visualize the graph. This issue has been already discussed with the conclusion that there will be one function plotting the top clusters and another opening cytoscape with a given graph. The decision about which metric to use to evaluate the graph still needs to be made. To be able to make a proper decision about this matter firstly it is necessary to figure out how the final graph will look like and for that we have to in the first place think about what we are actually generating. We generate a graph from a set of SNPs which could be obtained in different ways. The most common relation will be between two SNPs but there will be many other higher order relations. Therefore our graph will be composed of many isolated clusters of SNPs related between themselves. All previously exposed will be taken into account when deciding the metric from one of the exposed on Table 3, regarding that the graph will highly probable be composed by several disconnected networks, both betweenness and closeness centrality cannot be used because the distance between disconnected nodes from a graph is infinite thus the only remaining option is degree centrality. Metric Description Algorithm Reference Degree centrality The number of connections that a node has in a network. (v)eg(v)C D=d [38] Closeness centrality Inverse of the average length of the shortest paths to/from all the other nodes. (x) C=1 (y,x)∑ y d [39] Betweenness centrality The number of network shortest paths that pass through a specific node. (v)CB= ∑ s=v=t∈V/ / σst σ (v) st [40] Table 3 Most used graph theory metrics The final design of the workflow may look like in Figure 9. 25 If input object is a graph Get degree centrality from the graph Sort the degree centrality list Cut the list to the first N nodes Return List Else Error message Return End if Figure 12 Pseudocode of the function getTopDegree 4.6. openGraph This functions has two variants, the main one takes a graph, saves it and open it in Cytoscape, the variation from this function plot the top minimum value between six and the number of clusters as seen in Figure 13. The saving format is .gml as it is the format that share the package igraph and Cytoscape, there is no need to write the extension as if it is not or it not correct the function will place it automatically so it can be recognized by Cytoscape the pseudocode of both options is described in the following Figure 14. Figure 13 Plot representation from a sample graph using openGraph function 32 If input is a graph If Cytoscape is on the system path If filename is not ending on .gml Add .gml to filename End if Save graph Invoke to Cytoscape with the file path Else Send error message End if Else Send error message End if If input is a graph If Cytoscape is on the system path If filename is not ending on .gml Add .gml to filename End if Save graph Initialize i to the min of 6 or the size Divide the plotting area into i/2:2 For first i clusters Plot cluster on corresponding area End for Else Send error message End if Else Send error message End if Figure 14 Pseudocode of both versions from the function opengraph 4.7. Workflow implementation in Galaxy As part of the planned work for this project the port of the package to Galaxy is additionally one important part of it as it makes simpler to run experiments with no need to face difficulties of using the UNIX shell, although to port the R code to galaxy we need firstly to make it run in the shell. Other programing languages such as C or Java can run directly in the shell after being compiled, unfortunately this is not the case of R as is an interpreted language therefore it must be run using it specific tools, in this case the RScript command. To make it work in the shell, the software should be adapted considering that there is not an enviroment where to store variables, therefore the first thing we must do is change the return values for functions to save the data, in this case the function save from the package base  is enough to save R objects into a binary file, if what we want is to write a plain text file, as it happens in functions like getTopDegree  we must use the function write  . Those are not the only changes that had to be done, the messages should be written through stdout or stderr so Galaxy can recognize it, if a message is sent through stderr then it will directly finish the module with an error status, also to adapt it to receive parameters through the console a header must be added, including a checking for the number of parameters, to the function file. 33 In order to insert the code in Galaxy a wrapper is needed, this is created by inserting all parameters with their class together with the command to run the software previously adapted to run in a console. Outputs need to be detailed in a XML file as well including their class. Some complex classes which are not currently defined could be added to use it by adding it to a specific file for it. It is very important that the input and output classes from different modules must be the same in order to allow them to be connected by the users. After the porting is done the workflow must be created, a main workflow including creation of both graphs and both analysis and visualization of graphs has been designed as shown in Figure 15, alternatively there will be another workflow analyzing and visualizing only the gene graph as it is the final intention of the workflow. The user will anyhow have access to the modules to run them in an isolated way. The user can create its own workflow even changing the order or using different modules to complement the results, this facility is one of the advantages of porting the package to a WMS, allowing the reutilization of code. 34 Figure 15 View of the complete workflow in the Galaxy workflow editor canvas. 35 5. Use case and validation 5.1. Dataset The dataset used to perform this analysis is the one obtained during the study of Clarkson Disease [44], it consists of samples from 12 individuals with Clarkson disease, and 18 healthy control samples. Four of the patients were female, eight were male. In the control group, eight were female and ten were male. All of the patients were Caucasian, as it is the typical demographic of the disease, as were 16 of the 18 individuals in the control group. Genotyping was carried out with Affymetrix Genome-Wide Human Array 6.0 SNP chips. Prior to commencing the epistatic analysis of the dataset, single SNP association was carried out using the software package PLINK. Epistatic analysis was then carried out using the linear mixed model based FaST-LMM [45]. After more than three days of computing the exhaustive epistatic analysis an edge list of rsid pairs is obtained, this was shared by the Bitlab group so results can be replicated. The dataset consists of an edge list of 695 SNP pairs in plain text with header. It will be provided digitally. 5.2. Software benchmarking This section aims to replicate previously published results [44] to demonstrate the software is properly working First step is to create a SNP graph, this will be done by using the function createGraphFromTable. To analyze and compare it to the published results we will use functions getTopDegree and openGraph. The list of the top 20 ranked SNPs is shown in Table 13. Running the function createGraphFromTable required less than a second as well as getTopDegree. This is a reasonable time as other methods depending on external services may take longer execution time depending on the connection of the user. Then the association of the SNPs with the genes is performed using the function fillInfo and the result is shown in the Table 9. 36 Rank SNP Gene symbol Chromosom e Rank SNP Gene symbol Chromosome 1 rs11842523 PCDH8 13 11 rs7719321 MAP1B 5 2 rs4593336 NT5E 6 12 rs701170 PGBD5 1 3 rs17451360 SLC38A11 2 13 rs10191604 MYO3B 2 4 rs150533 LRRC16A 6 14 rs7420094 MMADHC 2 5 rs6927384 LRRC16A 6 15 rs12142665 PGBD5 1 6 rs1380237 MAP1B 5 16 rs7302874 CCDC64 12 7 rs6878132 MAP1B 5 17 rs11904398 MYO3B 2 8 rs7703322 MAP1B 5 18 rs12175817 GCLC 6 9 rs7704592 MAP1B 5 19 rs9382212 GCLC 6 10 rs7716699 MAP1B 5 20 rs10041715 RIOK2 5 Table 9 Result of the most significant SNPs and their related genes Comparing results with those obtained from the original publication the first fourteen elements match perfectly with the first fourteen shown in the original table from the publication as seen in Table 10. The last changes, where we can see three of the last six ranked in the publication, do not mean that results are incorrect. Differences in the results are due to the measure used and the size of the graphs. Degree centrality measure result is expressed in integers and therefore there can be many different nodes with same value. This problem increases when we are measuring small size graphs. In this case the results from rank fourteen have all degree centrality 6 being a total of 15 different nodes with this same score. Therefore amplifying the top to 29 would have all SNPs cited on the original study as can be observed on Table 11. 37 Rank SNP Gene symbol Chromosome Rank SNP Gene symbol Chromosome 1 rs11842523 PCDH8 13 11 rs7719321 MAP1B 5 2 rs4593336 NT5E 6 12 rs701170 PGBD5 1 3 rs17451360 SLC38A11 2 13 rs10191604 MYO3B 2 4 rs150533 LRRC16A 6 14 rs7420094 MMADHC 2 5 rs6927384 LRRC16A 6 15 rs10079905 RIOK2 5 6 rs1380237 MAP1B 5 16 rs10065590 RIOK2 5 7 rs6878132 MAP1B 5 17 rs10079905 RIOK2 5 8 rs7703322 MAP1B 5 18 rs10476724 RIOK2 5 9 rs7704592 MAP1B 5 19 rs12175817 GCLC 6 10 rs7716699 MAP1B 5 20 rs9382212 GCLC 6 Table 10 Result provided in the publication of reference Rank SNP Gene symbol Chromosome Rank SNP Gene symbol Chromosome 15 rs12142665 PGBD5 1 23 rs10476724 NA 5 16 rs7302874 CCDC64 12 24 rs12652555 RIOK2 5 17 rs11904398 MYO3B 2 25 rs1560323 NA 5 18 rs12175817 GCLC 6 26 rs6556989 RIOK2 5 19 rs9382212 GCLC 6 27 rs6556992 RIOK2 5 20 rs10041715 RIOK2 5 28 rs10248777 SNX10 7 21 rs10065590 RIOK2 5 29 rs10238880 SNX10 7 22 rs10079905 RIOK2 5 Table 11 Result of the SNPs with Degree centrality score six. Those showing NA did not find any gene within the flanking distance. In order to validate the results obtained from the gene graph only the first 12 ranked genes will be displayed rather than the 20 in the original study. This is done due to 38 the decreasement of the size of the gene graph, being from rank 12 the score two or less, that makes a total of 42 genes with a score of 2. As we can see in the comparison except of a switch between two genes, both scored 4 (SLC38A11 and TLDC1) and further changes in the order between those genes with score 3. The results are identical to those published as can be seen comparing Table 12 and Table 13. Rank Gene Chromosome Rank Gene Chromosome 1 PCDH8 13 11 EMP1 12 2 PGBD5 1 12 RIPK4 21 3 LRRC16A 6 13 CRYBA4 22 4 CCDC64 12 14 CHRD 3 5 TLDC1 16 15 TUSC1 9 6 SLC38A11 2 16 ETAA1 2 7 ASIC2 17 17 C16orf47 16 8 PSD3 8 18 GALR1 18 9 LHPP 10 19 TLE3 15 10 MAMLD1 X 20 MYO3B 2 Table 12 Results obtained on the reference paper Rank Gene Chromosome Scor e Rank Gene Chromosome Score 1 PCDH8 13 10 7 ASIC2 17 3 2 PGBD5 1 9 8 EMP1 12 3 3 LRRC16A 6 8 9 LHPP 10 3 4 CCDC64 12 5 10 MAMLD1 X 3 5 SLC38A1 1 2 4 11 PSD3 8 3 6 TLDC1 16 4 12 RIPK4 21 3 Table 13 Results obtained using the package epiR to replicate the results 39 To provide some biological meaning to the results obtained, an enrichment analysis must be performed. In this case we use WebGestalt, setting the enrichment type to GO enrichment, and then setting the test adjustment to Benjamini & Hochberg (BH) and at last defining the significance level as 0.1. We perform the analysis, obtaining the result after few seconds, stating that the two most affected functions are morphogenesis of epithelium as can be seen in Figure 16. When we go into depth by clicking on the boxes we can see that the genes related to these functions are PCDH8 and RIPK4, PCDH8 is the gene where more mutations has been detected so it makes sense that that gene is the cause of the illness, even more when according to some studies [46] state that the cause of the disease may be apoptosis of endothelial tissue, that means the death of the cell itself. This may be caused by a malfunction during the morphogenesis, therefore the regulatory mechanism from the cell would activate the cell’s apoptosis. Figure 16 Result obtained from the enrichment analysis 5.3. Use case The developed package will be used to perform a simple experiment to deduce whether the mutations affecting this biological function are located in the left flank or on the other hand on the right flank, To do this we will modify the flanking distance. We will perform an identical analysis but this time we will set the left flanking distance to zero and right to 500,000 bp and vice-versa after performing this experiment 40 together with the enrichment analysis we would be able to determine whether the left flank is more decisive on this disease or on the other hand if the right one is more decisive. In order to perform this we only have to run once more the function replaceSNPwithGene from the previous study reassigning there the flanking distance, this avoid the internet connection and speed up the process. After running the function we obtain the results shown in Table 14. Left Flank Right Flank Rank Gene Score Rank Gene Score 1 PCDH8 10 1 PCDH8 11 2 LRRC16A 8 2 LRRC16A 8 3 PGBD5 6 3 PGBD5 6 4 GALNT2 4 4 CCDC64 4 5 SLC38A11 4 5 SLC38A11 4 6 TLDC1 4 6 EMP1 3 7 CCDC64 3 7 GALNT2 3 8 DNAH17 3 8 DNAH17 3 Table 14 Results obtained from performing the analysis with 45 mpb The gene graphs obtained is shown in Figure 17 and in Figure 18, there we can visualize that PCDH8 is the gene most connected in the network, this provides a higher significance to that gene. 41 [32] Zhang, Y., & Liu, J. S. (2007). Bayesian inference of epistatic interactions in case-control studies. Nature genetics, 39(9), 1167-1173. [33] Calle, M. L., Urrea, V., Malats, N., & Van Steen, K. (2010). mbmdr: an R package for exploring gene–gene interactions associated with binary or quantitative traits. Bioinformatics, 26(17), 2198-2199. [34] Lippert, C., Listgarten, J., Davidson, R. I., Baxter, J., Poon, H., Kadie, C. M., & Heckerman, D. (2013). An exhaustive epistatic SNP association analysis on expanded Wellcome Trust data. Scientific reports, 3, 1099. [35] Schüpbach, T., Xenarios, I., Bergmann, S., & Kapur, K. (2010). FastEpistasis: a high performance computing solution for quantitative trait epistasis.Bioinformatics, 26(11), 1468-1469. [36]Huang, D. W., Sherman, B. T., & Lempicki, R. A. (2009). Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources.Nature protocols, 4(1), 44-57. [37] Zhang, B., Kirov, S., & Snoddy, J. (2005). WebGestalt: an integrated system for exploring gene sets in various biological contexts. Nucleic acids research, 33(suppl 2), W741-W748. [38] Freeman, L. C. (1978). Centrality in social networks conceptual clarification. Social networks, 1(3), 215-239. [39] Okamoto, K., Chen, W., & Li, X. Y. (2008, June). Ranking of closeness centrality for large-scale social networks. In International Workshop on Frontiers in Algorithmics (pp. 186-195). Springer Berlin Heidelberg. [40] Freeman, L. C. (1977). A set of measures of centrality based on betweenness. Sociometry, 35-41. [41] Ashburner, M., Ball, C. A., Blake, J. A., Botstein, D., Butler, H., Cherry, J. M., ... & Harris, M. A. (2000). Gene Ontology: tool for the unification of biology.Nature genetics, 25(1), 25-29. [42] NCBI2R, M. S. NCBI2R-An R package to navigate and annotate genes and SNPs. R package version,1(4). [43] Yates, A., Akanni, W., Amode, M. R., Barrell, D., Billis, K., Carvalho-Silva, D., ... & Girón, C. G. (2016). Ensembl 2016. Nucleic acids research,44(D1), D710-D716. [44] Upton, A., Trelles, O., & Perkins, J. (2015). Epistatic analysis of clarkson disease. Procedia Computer Science, 51, 725-734. [45] Listgarten, J., Lippert, C., & Heckerman, D. (2013). FaST-LMM-Select for addressing confounding from spatial structure and rare variants. Nature Genetics,45(5), 470-471. [46] Assaly, R., Olson, D., Hammersley, J., Fan, P. S., Liu, J., Shapiro, J. I., & Kahaleh, M. B. (2001). Initial evidence of endothelial cell apoptosis as a mechanism of systemic capillary leak syndrome. CHEST Journal, 120(4), 1301-1308. 48