enRoute: dynamic path extraction from biological pathway maps for exploring heterogeneous experimental datasets

Biological networks, such as interactions between proteins, biochemical reactions, and signaling processes are commonly depicted in pathway maps. Pathway maps are often hand-crafted and only show the part of the whole known biological network that is immediately relevant for a particular natural process, such as the tyrosine metabolism, or for a particular disease, such as HIV or diabetes. The network described by these pathways is based on published research on the interactions and interdependencies between the various nodes. As a consequence, pathway maps are static and are only valid for the specific processes or disease states they are designed for and fail to adapt to the variation found in real-world data. It is not uncommon, for example, that a de-activation of a node in a cascade invalidates reactions further downstream. For example, the gene PTEN is a part of the phosphoinositide 3-kinase signaling pathway, which regulates cell-growth [1]. If PTEN is mutated it does not fulfill its function and shuts down the pathway, which can lead to tumor growth. Jointly analyzing experimental data and pathways can help in reasoning about and predicting such effects for different conditions. Knowledge about how pathways are modulated by the genetic profile of groups or individual samples can help improving prognosis, treatment, and patient well-being.

Current approaches for visualizing interdependencies between pathways and experimental data do not scale to the now common large and heterogeneous experimental datasets, which often contain hundreds of experiments and multiple data types. We designed enRoute to remedy this. enRoute consists of two views: the pathway view, which shows the whole pathway and hints at interesting paths, and the enRoute view, which visualizes experimental data for parts of the pathway. In the pathway view, shown in Figure 1(a), we show the pathway maps augmented with abstractions of the mapping experimental data. Even though these abstractions are insufficient for an in-depth analysis, they provide an overview and hint at those parts worth investigating in more detail. The enRoute view shows the experimental data for a path that is selected in the pathway view (see Figure 1(a)). The selected path is extracted and juxtaposed with the experimental data, as shown in Figure 1(b). This combined approach successfully addresses the issue of showing large and heterogeneous datasets in the context of networks, the problem of showing multiple groupings of datasets, and it resolves multi-mapping issues that are common in pathway analysis. enRoute is part of Caleydo, an open-source biomolecular visualization framework (http://caleydo.org), which features various other visualization techniques for analyzing tabular and network data.

At the beginning of this paper, we give a brief introduction of the biological background, followed by a detailed analysis of the challenges of visualizing graphs with very large numbers of node-attributes. We continue by reviewing the literature and evaluating how existing approaches address the described challenges. Based on this discussion of the state-of-the-art and its limitations, we present our visualization technique, followed by a validation of our approach in case studies, conducted with experts in molecular biology. In the course of these case studies, we demonstrate how enRoute can be used to analyze large datasets in the context of pathways.

This paper is based on and extends previously published work [2]. In addition to a more detailed description of the original concepts, we extend the previous work with a generalization of enRoute to other pathway databases, a novel method to incorporate mutation status data into enRoute, a method to integrate potentially missing edges in the network, a semi-automatic path selection approach, an improved on-node data mapping in pathway maps, and various other improvements. In addition to the conceptual improvements and extensions, we present two new case studies that show the effectiveness of the enRoute technique.

Life scientists have accumulated intricate knowledge about biochemical and signaling processes in living cells, which has been used to build detailed biological interaction networks called pathways. These pathways summarize the molecular interactions in the biochemical conversions of molecules from source material to complex biomolecules or the signaling from a cell surface receptor via a second messenger to the nucleic transcription machinery. Several initiatives are drawing pathway maps and make these maps available to the scientific community, such as KEGG [3] or Wikipathways [4]. In the databases pathway maps are usually categorized by type (e.g., biochemical conversion or signaling pathways) and by biological purpose (e.g., cellular processes, human disease). Biochemical pathways describe the buildup or breakdown of molecules. For example, in the Glycolysis/Gluconeogenesis pathway of the KEGG database the buildup or degradation of glucose is described in great detail with all biochemical conversions and the connections to the Pentosephosphate pathway and the Citrate cycle pathway. A prominent example of the signaling pathway group is the MAPK signaling pathway, which contains a well-studied signaling cascade leading from an activated cell surface receptor to a phosphorylation cascade of several kinases to the activation of DNA binding complexes, which regulate transcription of genes involved in the proliferation of cells thus enabling the cell to react to growth stimuli from its environment. It is common for pathways to include the spatial organization of cells like cell walls, the Golgi apparatus, or the cell nucleus, thus allowing to depict transport of molecules and signaling through these compartments.

Pathway maps have become a valuable resource for molecular biologists summarizing broad knowledge about molecular interactions and presenting this information in a condensed view highlighting the functional aspects most interesting to the researcher. In the pathway maps nodes represent biological entities like proteins, metabolites, or chemical compounds. Protein nodes are annotated using the gene names from which these proteins are transcribed. These nodes usually contain several isoforms of the same protein and often several proteins of the same family catalyzing the same reaction, leading to extensive multi-mapping of many gene names to a single node. The nodes in a pathway are connected by links that depict biochemical reactions, activating or inhibiting modifications of proteins, or enzymatic reactions. The nodes together with the links between them allow to capture the interaction network of a biological system and present them in a structured graph with node-link diagrams following specific drawing conventions. Most pathway databases contain manually curated pathway maps that are trying to represent the molecular interactions in a visually appealing way while maintaining the scientific content.

In recent years, the advent of high throughput ~ omics technologies has generated large amounts of data. These datasets include genetic and expression data from large consortia like TCGA (http://cancergenome.nih.gov) and ENCODE (http://encodeproject.org/), but also high throughput metabolic screens via nuclear magnetic resonance (NMR) measurements or large scale proteomics by mass spectrometry. Data generated from these analysis includes dynamic data like gene expression levels, metabolite levels, and protein levels in various tissues, cells, and disease states, and also information about the genetic constitution of samples like copy number variation of genes, mutations in genes, or methylation patterns. All this data can only be interpreted in the context of the biological system present in cells, which are captured in the aforementioned pathway maps. The purpose of the techniques presented in this paper is to support the researcher in dynamically mapping biological data onto pathways. This allows researchers to compare, reason, and ultimately explain the complex biological systems and signaling cascades.

As a foundation for the design of our visualization technique and the evaluation of the related work, we have conducted a requirement analysis based on interviews with our collaborators from the Medical University of Graz. Our analysis resulted in five requirements that must be met by a visualization system to successfully enable the joint analysis of pathways and experimental data.

R I: The Scale Requirement - While the scalability of the graph is addressed by the sub-division of the biological network into individual pathway maps, the experimental datasets we consider are quite large. Consequently, a joint pathway and experimental data visualization system must be able to scale to dozens of experimental conditions or groups and hundreds of samples.

R II: The Heterogeneity Requirement - Modern biological studies often include a wide array of complementary but heterogeneous experimental datasets. While, for example, mRNA expression data measures the gene activity, copy number or mutation data can be used to reason about deviating expression values. These heterogeneous datasets need to be presented using different visualization techniques, as they differ in terms of data type. For example, mRNA expression data is numerical, copy number data is a hybrid categorical/numerical dataset, which is often binned into ordinal (ordered categorical) data, and mutation status data is nominal (unordered categorical). In order to analyze these different kinds of data in context of pathways, the visualization system needs to handle all of them simultaneously and also represent each of them using suitable visual encodings.

R III: The Multi-Mapping Requirement - Pathway nodes can represent various different gene products, such as enzymes, proteins, or RNA. In some cases, pathways also summarize a whole gene family into a single node, where different genes produce functionally similar proteins. This is what we call a multi-mapping: One node in a pathway actually represents multiple entities and therefore multiple entries of an experimental dataset can be associated with this node. As understanding this complexity is essential for judging effects of experimental data on a pathway, it is critical to convey multi-mappings adequately.

R IV: The Layout Constraint Requirement - Layouts for pathway maps can be either produced automatically or manually. To make pathway maps easier to understand, a number of drawing conventions have been established. Examples are cycles being drawn in circles or the predominance of orthogonal edges in many popular databases. Also, manually curated pathway maps typically contain rich meta-information indicating, for example, the cell compartments in which specific processes occur. Automatically drawn pathway maps can either try to respect these conventions (see, for example, the algorithm by Lambert et al. [5]), or optimize global layout properties using, for example, algorithms for force-directed layouts. We observe that domain experts prefer manual or at least consistent layouts that adhere to these drawing conventions. A reason for that might be that familiar layouts enable them to recognize and understand a pathway only by its topology. While showing experimental data is easier in automatically generated layouts, as the layout can be adapted to suit the representation, a good visualization technique for joint analysis of experimental data and pathways also needs to work with the large baseline of existing, manually produced pathways.

R V: The Topology-Attribute Coexistence Requirement - We distinguish between two main types of tasks conducted on a pathway: tasks that are based on the topology of the underlying graph, and tasks that are based on the node or edge attributes of the graph [6]. Topology-based tasks are concerned with the connectivity of the graph, e.g., which nodes can be reached from a given node, what are the articulation points of a graph, etc. An example for a topology-based task in pathway analysis is to find all nodes that might be influenced by the inhibition of a node at the beginning of a pathway. Attribute-based tasks are concerned with analyzing the properties of node or edge attributes. Edge attributes in pathways commonly describe the type of a relationship between two nodes, such as biochemical conversion, while mapping experimental data represents the majority of node attributes. An example for an attribute-based task for pathway analysis is to find all nodes in a pathway that are mutated in a large number of the mapping samples.

Visualization techniques for graphs are usually optimized for either topology-based task or for attribute-based tasks, but are rarely suitable for both at the same time. Node-link diagrams, are, for example, well-suited for topology-based tasks, while matrix layouts, where nodes are shown on the sides of a matrix and the cells contain information on whether there is an edge connecting the nodes, are ideally suited for edge-attribute-based tasks [7]. When analyzing pathways and experimental data, however, both types of tasks need to be addressed at the same time. The two central questions an analyst is trying to answer when analyzing both pathways and experimental data are (a) how the experimental data for particular experimental conditions or groups of samples influences the topology of the graph and (b) how effects observed in the experimental data can be explained using the topology of the pathways. Consequently, an effective visualization technique has to enable both: an in-depth analysis of the topology and the pathway attributes.

While there is a wide body of literature on graph drawing and graph visualization, we focus on the discussion of techniques that are either directly relevant for pathways or that can address the scale requirement (R I) with respect to the encoding of node or edge attributes. For a comprehensive review of systems biology visualization refer to the article by Gehlenborg et al. [8]. We identify several techniques that can be used to visualize multiple edge and node attributes in graphs. These are:

The benefit of using on-node mapping is that it makes it easy to address the layout constraint requirement (R IV). Consequently, on-node mapping has been widely used to augment pathways with multiple colored rectangles, each representing a single experiment or an aggregation of multiple experiments [9–11], where the color encodes the value. There are also variations that use color together with selection and animation [12]. The biggest drawback of this approach is its inability to scale (violating R I) as the amount of distinguishable colored rectangles inside a node is severely limited.

An alternative strategy to using multiple colors within each node in one graph is to use multiple graphs where each of them uses a single experiment or a single aggregate of experiments to drive its color-coding. This approach is commonly referred to as small multiples. Small multiples show the same configuration of a plot multiple times while changing one variable ([13]pp.170-175). An example that employs small multiples for automatically layouted pathways is Cerebral [14]. Lex et al. have used small multiples to show differences between experimental data associated with cancer subtypes on top of KEGG pathways [15]. Again, scale (R I) is a limiting factor. Depending on the pathway about four to ten multiples are reasonable.

A technique that can easily address the scale (R I) and the heterogeneity requirement (R II) is using separate linked views for the experimental data and the pathways. This, of course, also preserves the topology (R IV) and can be used to address multi-mappings (R III). Separate linked views use synchronized highlighting (linking & brushing) between the multiple views to communicate relationships. If, for instance, a user selects a node in the pathway, the corresponding experimental values are highlighted in the views depicting the experimental data. Shannon et al. [16] and Barsky et al. [14], for example, use a parallel coordinates plot for experimental data, which is linked to a graph depicting protein interaction and metabolic networks. Streit et al. [17] use heat maps and parallel coordinates to show experimental data related to pathways. The major shortcoming of the separate linked views approach is its failure to address R V, to simultaneously enable topology-and attribute-based tasks. As separate linked views require interaction to show relationships between a single node and its associated data, the joint analysis of the topology and attributes is severely hindered.

Finally, there are methods that adapt the graph layout to be able to show experimental data in pathways. There are numerous systems that calculate an automatic layout for pathways (violating R IV) and choose a node size that enables in-place encoding of experimental data with various visual encodings (e.g., bar [18, 19] and line charts [20]). While this approach scales a little better than simple on-node encoding, it fails to scale to larger numbers of experimental values (R I).

There are also more radical adaption approaches for the graph layout. Schulz et al. [21], for example, use two tables, one for each "side" of a bipartite network and connect the rows in the tables with edges. Each node is represented by one row and there are multiple columns for node attributes. GraphDice by Bezerianos et al. [22] probably takes the most extreme approach by laying out the nodes purely based on their node attributes in a scatterplot while still drawing the edges. Both approaches severely hamper the interpretability of the topology, violating R IV and R V.

A different approach on adapting the graph layout was taken by Meyer et al. with their Pathline tool [23]. Pathline uses a linearized version of a pathway where branches and cycles are conveyed using special visual encodings. Next to the linearized pathway the system shows the Curvemap view, which displays experimental data for both genes and metabolites recorded in time series. While Pathline was the main inspiration for our approach, it suffers from the unconventional pathway layout, which can hinder understanding the graph topology (R V). Also, it currently requires manual creation of the linearized pathways, thereby making it difficult to integrate the large existing databases of pathways.

The approach presented in this paper is implemented as a module of Caleydo, which is an open source visualization framework for biomolecular data [17]. Binary versions for Windows, Linux, and Mac, as well as the full source code and help material are available at http://www.caleydo.org. Both, our visualization technique and the visualization framework are implemented in Java and OpenGL using the freely available JOGL library (available at https://jogamp.org/jogl). In order to highlight selected pathways, we make use of a free implementation of the Bubble Sets technique [29] that we integrated into the Caleydo framework.

Caleydo allows analysts to import data sets with large amounts of experiments, such as mRNA, copy number, and mutation status datasets. We provide pre-packaged and preprocessed TCGA datasets (see http://tcga.caleydo.org), which contain up to 900 samples (depending on the subtype stratification) and ~ 20,000 genes each. As demonstrated in Figure 10, our visualization technique scales well to such large datasets. To cope with the limitations of available screen space, we make use of scroll bars for longer path in enRoute, which is justified by the linear nature of the exploration process along the path.

Our current implementation builds upon the well-established KEGG and Wikipathways databases, which provide information about pathways as image data in combination with an XML-based descriptions of the graph. The description includes the topological information as well as the position and the size of nodes. We use this information to render the various augmentations described on top of the pathways.

Since pathway databases such as KEGG do not provide information about how edges are routed in the images, directly connecting nodes in a path using straight edges can lead to overlaps between edges and nodes that are not part of the path. In Figure 9(a), for example, the edge highlight between node CASP9 and CASP3 collides with the representation of node CASP7, which is not a member of the selected path. In such cases, we route edges around nodes to avoid overlaps. Figure 9(b) shows the Bubble Set with a refined route from CASP9 to CASP3.

The Bubble Sets algorithm renders the overlays semi-transparent on the top of an existing base representation. However, if the base representation uses color-coding, overlaying colored Bubble Sets can lead to a wrong interpretation of the encoded information, as the colors of the Bubble Sets interfere with the color of the augmented content. Additionally, blending decreases the legibility of the original content. To avoid this problem, we cut the area of nodes out of the Bubble Sets.

Pathways can be selected from a drop-down list, which contains all loaded pathways. The full list is reduced on-the-fly to only contain the ones that match a user-specified query string. In addition to the title-based search, analysts can click embedded pathways to replace the current pathway with the full version of the embedded one.

We developed the enRoute visualization technique in close collaboration with a biologist from the Medical University of Graz, who is also an author of this paper. To evaluate the utility of the enRoute visualization technique, we conducted two case studies using different datasets together with this biologist. The first dataset is taken from the Broad-Novartis Cancer Cell Line Encyclopedia (CCLE, http://www.broadinstitute.org/ccle/home), which contains the genetic and pharmacologic characterization of a large panel of human cancer cell lines. The second is a microarray dataset from a model of hepatocellular carcinoma collected at the Medical University of Graz, which can be found through the corresponding original article [30].

Apoptosis regulation in cancer cell lines

The first case study investigates the difference in regulation of the apoptosis pathway in different human tumor cell lines from various organs. Apoptosis is the programmed death of a cell due to internal damage or as a consequence of external stimuli and involves a signaling cascade, which is not mediated by phosphorylation but by targeted degradation of proteins through enzymes called caspases. To initiate apoptosis in a cell, a ligand called TRAIL or TNFa binds to receptors on the outer cell surface. This leads to activation of receptor associated death domain containing proteins, which activate the first of a cascade of caspases, which in turn eventually leads to the cleavage of proteins in the cell subsequently triggering apoptosis. During apoptosis the cell shrinks in volume, exhibits nuclear fragmentation, chromosomal DNA fragmentation, and release of Cytochrome c from the mitochondria, which eventually leads to cell death.

As first step of the analysis, the Apoptosis pathway of the KEGG database is loaded into the the pathway view. The researcher selects TRAIL as the starting point and DEF40 as the endpoint of the signaling cascade to be investigated. As shown in Figure 10(a), the system automatically highlights all possible paths from TRAIL to DEF40 and allows to switch between these alternative paths. The selected path is extracted and displayed in a linear layout in the enRoute view, which shows detailed CCLE mRNA expression, copy number variation, and mutation status data for the genes that map to the nodes of the path, as can be seen in Figure 10(b). In this example, the grouping of the experimental data reflects the source organ of origin of the cell lines. The analyst can now easily identify differences in the copy number variation, gene expression, and mutation status for these groups and relate them immediately to the stages of the linearized path, which would not be possible using other techniques. As indicated by the upward pointing bars for the copy number data, there is a clear amplification of TNFSF10 (an alias for TRAIL) for cell lines originating from ovary, lung, or breast tissue among others. By selecting all samples with an increased copy number of TNFSF10 in the histogram, the system also highlights these cell lines in the mRNA expression plots, as indicated by the arrows in Figure 10(b). The analyst can now confirm that ovary cell lines with a higher copy number of this gene also show a higher expression of TNFSF10 in the mRNA expression plots. Based on this data, one could interpret that these cell lines should be prone to apoptosis, as they have high amounts of the cytokine triggering the cell death. However, given the nature of these cell lines which are derived from human tumors, which are by definition resistant to apoptosis, this conclusion would be misleading. An explanation for this seemingly strange behavior can be found by looking at the copy numbers of the genes TNFRSF10A to TNFRSF10D, which are receptors for TNFSF10 cytokine, further downstream the path. Here, it quickly becomes evident that these genes exhibit a frequent loss of copies, which is easy to see when looking at the many blue bars pointing downwards for these genes. This explains how these cells can still form a tumor, as the loss of the receptor makes the cells immune to the stimulus by the TNFSF10 cytokine, thus preventing apoptosis initiation.

In summary, enRoute enabled the analyst to study copy number variation, gene expression, and mutation status of a large number of samples in parallel in a clearly visualized and linearized sub-path of the complex network initiating apoptosis. The researcher stated that the linearization of the signaling cascade and the clear association of the other data entities to the individual steps of this cascade are very intuitive and greatly facilitated the interpretation of the data and the deduction of a biological interpretation from the dataset.

HCC xenograft models

The second use case demonstrates the ability of the enRoute path extraction feature to aid in a common problem when interpreting pathway maps. In biological reaction systems enzymatic reactions are often carried out not by a single gene or protein but by a family of proteins encoded by several genes. Additionally, enzymatic reactions are often not specific to a single protein family, but can be substituted by other enzymes, which, however, often work with different efficiency. This biological diversity leads to functional nodes in pathways that contain a multitude of genes with historically designated gene names often not indicative of the real function of these genes. It is extremely difficult for the researcher to keep in mind which genes are behind a single node of a pathway and thus understand the mapped biological data. Cross referencing of expression levels to genes and functional nodes can thus only be achieved by resolving the multi-mapping in a node to gene names which can then be mapped to the individual gene expression levels. In enRoute, the researcher can select nodes upstream and downstream of the enzymatic reaction, which is then resolved in a linearized representation of all genes involved. An example is the conversion of all-trans-Retinoate to all-trans-18-Hydroxy retinoic acid, which is a reaction contained in the Retinol (Vitamin A) metabolism pathway shown in Figure 11(a). The single node that represents this reaction contains 16 different proteins or genes that are involved. As depicted in Figure 11(b), enRoute resolves this node into a convenient map of genes, thus making it possible to map the gene expression of each individual gene in all experimental conditions of this experiment. The dataset mapped onto the Retinol pathway was generated by gene expression profiling of normal and cirrhotic human liver, hepatocellular carcinoma, and grafts of three human tumor cell lines (Hep3B, HUH7, and SK-Hep) into immunodeficient mice. Expression was measured in liver samples from patients and in samples from tumor cells grown in culture (TC), as subcutaneous grafts (SC), or as orthotopic grafts (Ortho) in the liver of experimental animals. The task was to find out how well the cell line models correlate to the human disease, with special focus on drug metabolism and oxidative stress. The response of tumors to anticancer treatments is closely linked to the activity of cytochrome p450 enzymes, which metabolize drugs and mediate oxidative stress. When studying the aforementioned conversion of retinoate to retinoic acid by cytochrome p450 enzymes, it can now be seen that the node CYP2A6 actually contains members of the whole CYP2 family of genes and that the expression of these genes is uniform in normal liver and cirrhosis, but very variable in HCC, dividing these cases into a low expressing and a high expressing group. It is immediately visible that the xenograft models are only representative of the low CYP2 expressing group of HCCs. Decomposition of the complex nodes allows the identification and investigation of the individual expression patterns of genes contained in the complex node. It becomes obvious that the expression of CYP2C18 (highlighted in yellow), contained in the node labeled CYP2A6, is higher in the Hep3B model than in the other cell lines.

Additionally, it can be detected that CYP1A1 is highly expressed in HUH7 tissue culture cells. All this information was not visible to the researcher using conventional on-node mapping approaches and was successfully visualized using enRoute. Node decomposition is an integral feature of the enRoute path extraction and thus makes the association of many mapped nodes and their corresponding experimental data readily available to the researcher.

Enabling the joint analysis of biological pathways and large amount of experimental data by means of visualization is a challenging problem. In this paper, we provide a list of requirements that have to be met to support such analysis. We introduce the enRoute visualization technique, which addresses these requirements using a tightly-coupled dual-view approach. Our visualization technique allows experts to select a single path from a pathway map, which is then extracted and shown in linear form in a second view, the enRoute view. This view allows experts to investigate associated experimental data in detail by displaying it in a tabular layout right next to the extracted path. The conducted case studies using two different datasets confirmed the utility of our visualization technique.

Although we currently already support pathways from both the Wikipathways and the KEGG database, we plan to integrate further resources for biological networks. For instance, an integration of the EBI IntAct database [31] will make the analysis of protein interaction networks possible using enRoute.

enRoute allows to select nodes independent of the structure of the pathway graph. However, the boundaries of a single pathway graph also limit the potential of enRoute. Consequently, we plan to extend enRoute to also support the analysis of paths spanning multiple pathways, which introduces new challenges for scalability concerning the topology of the graph. Solving these problems is subject of future research. Another interesting area of research is the selection of pathways, which we have only briefly touched in this paper. We intend to implement common measures for identifying interesting pathways in the context of the available experimental data and integrate them in a unified data and pathway visual analysis workflow.