Building pathway graphs from BioPAX data in R.

Biological pathways are increasingly available in the BioPAX format which uses an RDF model for data storage. One can retrieve the information in this data model in the scripting language R using the package rBiopaxParser, which converts the BioPAX format to one readable in R. It also has a function to build a regulatory network from the pathway information. Here we describe an extension of this function. The new function allows the user to build graphs of entire pathways, including regulated as well as non-regulated elements, and therefore provides a maximum of information. This function is available as part of the rBiopaxParser distribution from Bioconductor.

Introduction

Biological pathways represent signalling and/or metabolic events involving protein and non-protein molecules. They are increasingly used in gene and protein expression studies to provide an aggregate score for gene sets encoding for defined biological events [1]. Several pathway databases, either curated or not, have adopted the BioPAX [RRID:SCR_009881] (Biological Pathway Exchange) language as a standard for pathway representation using the RDF (Resource Description Framework) data model [2].

The structure of BioPAX is founded upon groupings, called classes, for physical entities and interactions with hierarchical networks of their sub-classes. Interactions between physical entities are represented such that conjoint interactions may form a specific pathway with defined, but different types of interactions between the involved physical entities. The BioPAX format is being actively developed, with BioPAX level 2 format focusing on metabolic pathways and BioPAX level 3 introducing full support for signalling pathways.

SPARQL (Simple Protocol And RDF Query Language) is a query language able to retrieve and manipulate data stored in RDF. Pathway information is often combined with statistical data analysis using tools such as R [3]. The rBiopaxParser [RRID:SCR_002744] [4] is an R package to retrieve data stored in a BioPAX RDF format. It comes with several options that are useful to probe the data and extract specific information from it, for example participants of a pathway, stoichiometric conditions to be fulfilled for an interaction, etc.

One such option is the pathway2RegulatoryGraph (P2RG) function that converts a pathway into a graphical structure. This is extremely useful for visual representation and subsequent graph-based network analysis. The P2RG function returns the parts of a pathway that are regulated (activated or inhibited) by proteins or protein complexes; this is important to understand the role of regulated proteins in pathways. Here we present an adaptation of P2RG, denoted pathway2Graph (P2G) which can be used to build a graph of the entire pathway, including the regulated as well as the non-regulated elements. This new function expands P2RG and can be used to investigate all different types of processes and connections of pathways instead of only studying the regulated elements of pathways. P2RG retrieves regulatory interactions, such as inhibitions and activations (shown in Figure 1 as continuous edges). The new P2G additionally, retrieves protein modifications, such as translocations or complex formation, which are shown as discontinuous edges in Figure 1.

Figure 1.

Hypothetical pathway.

This cartoon of a pathway shows examples of nodes and edges that could be encountered in a BioPAX database. The nodes are proteins, complexes or other physical entities and the edges are interactions between the nodes, that represent either interactions among proteins or protein modifications. The solid edges are those detected by the P2RG function and the solid and dashed edges are detected by the P2G function.

P2G is specifically aimed at retrieving results from Reactome BioPAX level 3. In this paper we describe detailed information on this function which, we believe, will help rBiopaxParser users to better understand the graphs generated from pathway information. We have verified P2G results by directly querying the original BioPAX data using SPARQL.

Methods and results

The classes of PhysicalEntity and Interaction that are used in Reactome v51 to represent information on pathways are shown in Figure 2. This graph was generated using the tool RDF2Graph [5] on the Reactome Level 3 RDF file. The nodes in Figure 2 represent classes and the edges show the possible relationships, called predicates, these classes could have in the database. As depicted in Figure 2, the node Pathway could have one or more PathwaySteps that consist of different types of Interaction sub-classes. All the Interaction nodes shown in Figure 2 describe interactions between PhysicalEntities, hence are connected to them by particular types of predicates as indicated in the edge labels. The Interaction classes are interconnected because they can be dependent on each other. The Control interaction and its sub-classes ( Catalysis and Modulation) represent signalling events. They regulate BiochemicalReaction and Degradation interactions which mostly represent metabolic reactions.

Figure 2.

Interplay of classes in Reactome BioPAX.

This figure shows a network of the Interaction and PhysicalEntity classes that are a part of any pathway in Reactome v51 BioPAX level 3. Nodes are classes and the directed edges are links between them in the database. The green nodes are the Pathway and PathwayStep classes, the blue nodes are Interaction classes and orange nodes are PhysicalEntity classes.

To create a regulatory graph, the P2RG function starts with the Control, Catalysis and Modulation interactions that are either activating or inhibiting other interactions. This method provides a graph with plenty of information on the regulatory components of the pathway. The nodes of this graph are physical entities like Proteins or SmallMolecules and the directed edges are either activation or inhibition events. An example of such a reconstruction is shown in Figure 1, where P2RG is able to retrieve the black nodes and the continuous edges. However, interactions can be missed if they are not regulated by the Control interactions and could result in the loss of valuable information in the graphical representation of the pathway.

The new function P2G can start with any type of interaction in order to obtain a graph with all possible physical entities involved in the pathway. Similar to the result of the P2RG function, the P2G function gives a graph with nodes that are physical entities, but the edges are not strictly activation or inhibition events. The directed edges could represent several types of events like translocation of a protein or cleavage of DNA, these are shown as dis-continuous edges in the cartoon in Figure 1. The P2G function recognizes the continuous and the dis-continuous edges and thus retrieves the black as well as the red nodes shown in Figure 1. In some cases there is more than one documented connection between the same physical entities. In this case only the first connection is used as an edge in the final pathway graph.

Comparison of two methods: P2G vs P2RG

The Reactome database (v51) categorizes pathways into 27 branches. Here we worked with pathways that have more than one interaction, which resulted in 1,666 pathways. Using P2RG, graphs for 1,548 pathways were retrieved. By using the new P2G function, we were able to retrieve information on all 1,666 pathways. The highest number of pathways were obtained, using either method, in the “Disease” category (P2RG: 3,396 pathways, P2G: 4,888 pathways). In 85% of the cases, pathways retrieved using P2G consisted of more physical entities (nodes) than those retrieved using P2RG. 19% of the P2G retrieved pathways have at least twice the number of nodes, and 60% have at least twice the number of interactions between nodes (edges) as compared to the P2RG version, Figure 3 is an example of this difference. Total numbers of nodes and edges in major Reactome categories are given in Table 1. Missing information causes the appearance of disconnected graphs when reconstructing pathways. By using the new P2G function, the percentage of disconnected pathways is reduced by 9%. Additionally, P2G also has the option of only retrieving the largest connected component, for example with this option, in Figure 3. A only the top left part of the graph will be retrieved and the two other disconnected parts discarded. The pathways have directed edges because most of the interactions have direction. Edges without a direction are represented as bidirectional edges in the output of P2G.

Figure 3.

Graphs of the pathway ‘Apoptosis induced DNA fragmentation’.

Both graphs were extracted from the same BioPAX file. A) Graph recovered using the new P2G function; B) Graph recovered using P2RG function. In both panels blue nodes are proteins or protein complexes, white nodes are non-protein entities. Black encircled nodes are found in both graphs and red encircled nodes are only detected with the new P2G function. Names of the nodes are in Table 2.

Table 1.

Numbers of nodes and edges.

The number of nodes and edges of ten different pathways (Reactome Categories) are indicated as obtained after application of P2RG and P2G on the same set of BioPAX RDF information.

Reactome Categories	P2RG Nodes	P2RG Edges	P2G Nodes	P2G Edges
Binding and Uptake of Ligands by Scavenger Receptors	0	0	68	56
Cell-Cell communication	13	14	142	142
Disease	3,396	5,878	4,888	12,159
Gene Expression	652	900	1,110	2,450
Immune System	1,431	2,233	2,419	5,045
Membrane Trafficking	86	121	181	382
Metabolism	3,082	5,922	3,479	11,289
Signaling Pathways	2,069	3,274	3,430	7,131
Steroid hormones	72	147	81	333
Transcription	281	420	623	1,324

Table 2.

Node names and locations of the “Apoptosis induced DNA fragmentation” pathway.

The first column has the names of the nodes in the pathwayas depicted in Figure 3. The second column has the actual name of the node and the third column the cellular location of the node. All this information is represented as given in Reactome version 51. The nodes shown with a black outline in Figure 3 are shown here in bold font.

Node	Name	Location
Protein8776	DFFB	Cytosol
Protein8777	DFFA	Cytosol
Complex4232	DFFA : DFFB	Cytosol
Complex4233	Importin alpha : Importin beta	Cytosol
Complex4234	DFF : associated with Importin alpha : Importin beta	Cytosol
Complex4235	DFF : associated with Importin alpha : Importin beta	Nucleoplasm
Complex4169	Active CASP3	Cytosol
Complex4238	Active CASP3	Nucleoplasm
Complex4236	DFFA : DFFB	Nucleoplasm
Complex4239	Caspase cleaved DFFA	Nucleoplasm
Complex4240	Caspase cleaved DFFA : DFFB	Nucleoplasm
Protein8779	DFFB	Nucleoplasm
Protein8784	DFFA fragment	Nucleoplasm
Protein8785	DFFA fragment	Nucleoplasm
Protein8783	DFFA fragment	Nucleoplasm
Complex4241	DFFB homodimer	Nucleoplasm
PhysicalEntity567	DFFB homodimer/homooligomer	Nucleoplasm
Complex2061	Histone H1 bound chromatin DNA	Nucleoplasm
Complex4242	DFFB associated with chromatin	Nucleoplasm
Protein8786	HMGB1/HMGB2	Nucleoplasm
PhysicalEntity109	DNA	Nucleoplasm
Complex4243	HMGB1/HMGB2 – bound chromatin	Nucleoplasm
Complex4244	DFF cleaved DNA	Nucleoplasm

As an example, we discuss here the ‘Apoptosis induced DNA fragmentation’ pathway, which is in the “Programmed Cell Death’ category ( Figure 3). When the information in the BioPAX file is reconstructed with the P2RG function, the pathway has seven nodes ( Figure 3.B); with the P2G function the same pathway has 16 more nodes ( Figure 3.A). Detailed information on these nodes, as retrieved with P2G and P2RG, is given in Table 2 and demonstrates the additional information retrieved by P2G. The node ‘Complex4169’, which is found in the cytosol, translocates to the nucleus where it is called ‘Complex4238’. However, this information is only available from the P2G function because the node ‘Complex4169’ does not regulate any other interaction or node. The presence of extra nodes in the P2G retrieved graph ( Figure 3) also visualizes that ‘Complex4240’ breaks up into its’ individual components after being cleaved by Caspase-3 (‘Complex4238’). This extra information is very useful for researchers analysing the phenomena represented by the pathway. In case P2G retrieved pathways graphs are used for analysis (e.g, differential gene expression analysis) the presence of these extra nodes may improve biological interpretation of experimental data.

Conclusion

P2G is a useful addition to the rBiopaxParser package because it retrieves all the components of a pathway from the database and provides complete graphical information for both signalling as well as metabolic pathways. The P2G function ( pathway2Graph) is currently available in the rBiopaxParser package in the Bioconductor 3.4 release.

Data availability

The data referenced by this article are under copyright with the following copyright statement: Copyright: © 2016 Benis N et al.

The input data for this package is the BioPAX format of any pathway database. We used the Reactome database which is freely available for download in different formats from the website www.reactome.org. A subset of this database is given as Supplementary file 1.

Software availability

Software available from: The function pathway2Graph is available in the latest version of the R package rBiopaxParser and can be installed from Bioconductor.

Latest source code: https://github.com/frankkramer-lab/rBiopaxParser/tree/2.12.0

Archived source code as at the time of publication: http://dx.doi.org/10.5281/zenodo.61618 [6]

Software license: GPL-2

The authors have addressed the (minor) issues I had with the article previously, thanks!

I have read this submission. I believe that I have an appropriate level of expertise to confirm that it is of an acceptable scientific standard.

The authors of this paper describe a new function provided by the rBiopaxParser library, which is a R based system for parsing BioPax documents. BioPax is coded in RDF, which is a linked data format that describes the subject matter using graph triples in the form of a Subject, Predicate and Object. Their pathway2Graph and the old gathwayRegulatoryGraph essentially attempt to do graph transformations, taking the RDF graph and converting it to a graph that is suitable for analysis.

The paper’s primary metric for demonstrating improvements in their method is measuring the number of non-zero pathways retrieved after extracting data from Reactome. In addition, I tested code and their supplied supplemental data.

```

> library(rBiopaxParser)

> biopax=readBiopax("19c8ac7b-96b2-4db4-a78a-c2defed535ae.owl")

> a <- pathway2RegulatoryGraph(biopax, "Pathway1020")

> b <- pathway2Graph(biopax, "Pathway1020")

> nodes(a)

[1] "CASP3"                        "CASP3(176-277)"

[3] "DFFA(1-?)"                    "DFFB"

[5] "DFFA(118-331)"                "DFFA(1-117)"

[7] "DFFA(225-331)"                "DFFA(118-224)"

[9] "Histone H1"                   "DFF40 homodimer/homooligomer"

[11] "HMGB1/HMGB2"

> nodes(b)

[1] "Protein8776"       "Protein8777"       "Complex4232"

[4] "Complex4233"       "Complex4234"       "Complex4235"

[7] "Complex4169"       "Complex4238"       "Complex4236"

[10] "Complex4239"       "Complex4240"       "Protein8779"

[13] "Protein8784"       "Protein8785"       "Protein8783"

[16] "Complex4241"       "PhysicalEntity567" "Complex2061"

[19] "Complex4242"       "Protein8786"       "PhysicalEntity109"

[22] "Complex4243"       "Complex4244"

```

From the old version (a) to the new extraction (b), for this particular pathway the number of nodes went from 11 to 23. Interestingly, the names of the vertices became less descriptive in the newer method, going from gene names like ‘CASP3’ to strings extracted from the RDF urls, like ‘Protein8776’.

It would seem that the authors claim, that the number of extracted elements from the BioPax files is accurate. But what is slightly confusing is the significance of this change. Was the original method faulty and this is a bug fix? Or is the method being updated to deal with the BioPax standard as it changes from version 1 to 2 to 3?

This seems like a minor but necessary change to the library to make it have a better extract information from BioPax. More explanation about the nature of this change, and how the new parsing strategy improved how the library was able to deal with all of the various types of relationship classes in the BioPax format would make a big difference in helping to illustrate the improvements mentioned in this paper.

I have read this submission. I believe that I have an appropriate level of expertise to confirm that it is of an acceptable scientific standard, however I have significant reservations, as outlined above.

The authors have developed a new function that allows the user to build a regulatory network in a graph format based on pathway information. In the version that the authors developed, the output graph includes regulatory and non-regulatory interactions and allows the viewer to more fully comprehend the underlying network. Nodes in the network represent classes, while the edges show the relationships among these classes. An example of the approach is provided with the Reactome database. The authors’ approach, P2G, was used to analyze data for 1,666 pathways. P2G returned more nodes than was retrieved by the earlier version in 85% of the cases. The software will be available as part of the Bioconductor 3.4 release. This will likely be a valuable addition to the Bioconductor package that will provide scientists with a means for generating graphical and intuitive networks from gene expression and metabolic data. The manuscript is clear with an appropriate title and abstract. The article is clearly written and the conclusions are based on the data.

I have read this submission. I believe that I have an appropriate level of expertise to confirm that it is of an acceptable scientific standard.

Thank you for your comments.

The authors have developed a new function in the rBiopaxParser package to generate figures from BioPAX formatted biological pathway data. This new function, called pathway2Graph (P2G) replaces an older function called pathway2RegulatoryGraph (P2RG). P2G includes more interaction terms than P2RG, and therefore generates more complete interaction graphs. That said, as written the changes are mysterious - why did the original function limit the interaction types when generating graphs and what is advantageous about including more interactions? Addition of the P2G function is a small improvement to the rBiopaxParser package that will be useful, but the discussion of P2G's advantages and disadvantages should be expanded.

Specific points:

The language throughout is highly technical, potentially compromising its readability to end users.

The edge labels in Figure 1 are very difficult to read. As these are discussed in the text, their font size and/or weight should be increased.

It would increase the readability of this software tool article if the authors described a (biological) scenario where the P2G function would be uniquely useful compared to P2RG. It’s obvious that including more interaction types will lead to more complete graphs, but in what scenario would this be useful for a user?

Why were all interaction types not originally included in P2RG? The advantages and any potential disadvantages of including all interaction types should be discussed more.

Currently, the difference between P2G and P2RG is rather minor (1,666 vs 1,548) – might this difference change in the future if more edges are added through the interaction types that are unique to P2G?

The work is technically sound and presents a useful extension to the rBiopaxParser package, but the paper describing this work will be useful to a broader audience if changes similar to those above are incorporated.

I have read this submission. I believe that I have an appropriate level of expertise to confirm that it is of an acceptable scientific standard, however I have significant reservations, as outlined above.

Thank you for your comments. In the new version we have expanded the Introduction with non-technical specifics that should explain the basic differences between the functions to a broader audience. Figure 2 (Figure 1 in the previous version) has been changed to increase visibility of edge labels. In order to explain the differences between the outputs of the two functions we have added a section in the Methods and Results section where we describe a particular pathway (‘Apoptosis induced DNA fragmentation’) in terms of the extra information gained by using the new function. A new table (Table 2) with more biological information on this pathway has also been added.

This articles describes the addition of a new function to the extant rBiopaxParser R library. This new function converts a BioPAX-formatted pathway of gene or protein interactions into a graphical structure that is human-viewable. This function supercedes an earlier function which performed the same task, but was unable to convert all nodes and egdes in the pathway. This resulted in a loss of information. This loss is now remedied, however, with the new and improved function.

I installed the package and tested it and it seemed to work flawlessly, although I didn't try anything fancy. I imagine that the new function is a significant advance for people working on these pathways routinely as it must have been frustrating to have missing data.

My only reservation about the paper is that I was unclear on the point of the new function until I'd read most of the paper. Perhaps the authors could extend the last paragraph of the introduction to be more clear about why the new function was necessary and how it has improved the R package. Although this request is a small change, the article is confusing as it is and making this would be a big improvement.

I have read this submission. I believe that I have an appropriate level of expertise to confirm that it is of an acceptable scientific standard, however I have significant reservations, as outlined above.

Thank you for your comments. The Introduction has been expanded with an image (now Figure 1) to emphasize the differences between the two functions.