Design criteria

The following features were considered core to the design of Graphia:

• Data and operating system agnostic. Import data from any source saved in standard file formats. The software should run on all major desktop operating systems and modern hardware configurations.
• Fast and scalable. Support the rapid loading of data, fast computation of graph layout and analysis algorithms, high quality data visualisations. Deliver smooth and responsive graphical rendering of millions of data points (node/edges) on standard desktop hardware.
• Dynamic rendering. Visualise in real time changes to the graph structure associated with alterations in input parameters or additional data.
• 3D graph visualisation. Provide a navigable and immersive environment in which to explore and interpret large and complex graph topologies.
• Correlation graphs as an essential function. Rapidly convert any numerical or categorical data table into a correlation graph, supporting pattern finding and data mining.
• Attribute handling and visualisation. Visualise attributes (metadata) associated with nodes and edges using colour, size and text to distinguish between attribute values.
• Advanced analysis capabilities. Support a wide range of analytical algorithms and approaches that empower a user to explore, query and interpret data, such as the k-NN algorithm [22] for edge pruning, and the MCL [23] and Louvain [24] graph clustering algorithms.
• Extensible. Provide extensible architecture through use of a plugin system to allow the core to be extended or adapted for specific application areas or data types.
• User Interface. Provide a simple and intuitive user interface (UI) that is easy to navigate, featuring a graph display area supplemented with a table listing selected nodes and associated attributes and data values. The UI should provide easy access to menus providing functionality and display active transformations and visualisations (Fig 1).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 1. Graphia user interface displaying a correlation graph.

(A) Graph display area, showing correlation graph with a cluster selected (unselected nodes faded). (A1) Display context menu options (right click). (B) Node (row) attribute display area. (B1) Table of selected nodes and their attributes, imported and calculated within the tool. (B2) Data plot area, in the case shown here a mean histogram of selected node values. (B3) Visualisation of column annotations. (B4) Data plot context menu options for changing plot (right click). (1a) Add Transform button, (1b) Active transforms, (2a) Add visualisation button, (2b) Active visualisations, (3) General toolbar, (4) Attribute parameter selection, (5) Display of graph metrics (number of nodes, edges, components), (6) Plot/table function toolbar.

https://doi.org/10.1371/journal.pcbi.1010310.g001

Code architecture

Graphia is written in C++17 and is built upon Qt version 5, the cross-platform widget toolkit. For graphics, the industry standard OpenGL is used. The minimum driver support required is version 3.3 core profile, but more modern extensions will be used if they are available. Various open source libraries are employed, mostly for loading external data formats. These libraries and their associated licenses are enumerated in the About dialog of the application, accessible from the Help menu.

Graphia is architected so that loading, and data type specific user interfaces are confined to plugins. These are independent modules to the core application and can be removed or added without affecting any base functionality.

At the highest level the code is organised hierarchically into four separate directories:

• app—the core application code
• plugins—the existing bundled plugins
• shared—code used by both the core and plugins (this includes interface headers)
• thirdparty—any library code not authored locally

These are further divided into subdirectories dealing with specific areas of functionality. Graphia has been developed using standard object-orientated best practices. Continuous integration is employed to prevent portability build regressions, using recent versions of the compilers GCC, clang and MSVC. In addition, static analysis tools such as clang-tidy and cppcheck are used to identify potential problems early. CMake is used as a build system, and is set up for Linux, Windows and macOS compilation.

Implementation

Dynamic graph layout and rendering.

Graph layout is an iterative process. Many programs only display the results of a layout algorithm after it has run a defined number of iterations. With Graphia, the layout is shown live, such that graphs ‘unfold’ in real time. However, the true power of dynamic graphs is realised when a transformation operation is performed or following the addition of new data. These changes are immediately reflected in the appearance of the graph. As graphs change dynamically, there is a need to identify the graph components and map how they interact when construction parameters are adjusted. The ability to quickly identify and move between components is a unique feature of Graphia. Components are rendered in a concentric pattern, arranged large-to-small. Smaller components may be filtered away using a transform. Most existing network analysis tools render graphs in 2D. The layout algorithm implemented in Graphia is innovative in that it applies current force-directed layout techniques, but in a dynamic fashion. Graphia renders graphs in 3D or 2D, making use of modern graphics hardware to display extremely large graphs efficiently with options for node/edge shading, relative node sizing and spacing (Fig 2A–2F). When graphs are relatively small or there is a need to share images by a conventional medium, i.e. a document, 2D graph visualisations have advantages. However, 2D visualisations are limiting when there is a need to display and interact with large graphs with complex topologies.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. Different graph visualisation options.

(A) 3D perspective view, smooth shading (the default), with visualisation of node categorical attribute (MCL cluster). (B) 3D orthographic view, flat shading (no perception of distance—all nodes same size, unless sized by attribute value). (C) 3D perspective view, smooth shading. (D) 2D view, smooth shading. (E) 2D view, flat shading. (F) compressed 2D layout, flat shading, showing node overlap view. Visualisation of (G) Betweeness centrality values, (H) Eccentricity values, (I) PageRank values. G-I are continuous (numerical) attributes, so a colour spectrum and size gradient is used for node display (2D, smooth shading). Betweenness and eccentricity are calculated for both nodes and edges, therefore visual encoding is applied to both.

https://doi.org/10.1371/journal.pcbi.1010310.g002

Case Study 1: Visualisation of phylogenetic trees

Hierarchical data structures are often represented by tree graphs and used in biology to represent relationships between species, strains, samples or genes. While trees are an intuitive way of visualising such relationships, when the number of branches on the tree become large, the ability to display such graphs at a local or global level is challenging. Here we show two examples of taxonomic trees visualised by Graphia representing the different levels of phylogeny, from a central node representing the class of organisms, up through branches representing the order, family, genus, with species and subspecies being the leaves of the tree. The examples described are taxonomic trees for all mammals and insects (Fig 3), as defined by the NCBI Taxonomy database [25]. The taxonomic tree of mammals consists of 9843 nodes and 9862 edges and is shown in 2D with nodes coloured by type, i.e. what level of the taxonomic tree they represent (Fig 3A). A small section associated with apes is highlighted (Fig 3B). When the graph is loaded using the WebSearch plugin, selection of a single node automatically initiates a search for the name of the selected node Fig 3C. Shown in Fig 3D is a taxonomic tree of all insect species, consisting of 275,328 nodes and 275,528 edges displayed in 2D. Graphia provides a unique environment to search, cluster and explore such data, the third dimension greatly assists in visualising the structure of such large graphs.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. Visualisation of taxonomic trees.

(A) A taxonomic tree of all mammals was downloaded from the NCBI’s Taxonomy database, with nodes coloured according to type, i.e. subspecies (blue), species (pink), genus (orange), etc. The graph comprised of 9,843 nodes and 9,862 edges and is shown with a 2D layout. (B) Zoomed-in view of the area in square shown in A, with a single node selected (Western gorilla). (C) Right click on a node provides the ability to search the web for the node identifier either via Google or through a predefined database selected in the Network tab of Options dialogue. (D) Taxonomic tree of all insects from the NCBI’s Taxonomy database, nodes coloured by Louvain cluster. The graph consists of 275,328 nodes and 275,528 edges and in this respect represents a large graph where visualisation in 2D is challenging.

https://doi.org/10.1371/journal.pcbi.1010310.g003

Case Study 2: Analysis of single cell transcriptomics data

Single cell RNA sequencing (scRNA-Seq) generates gene expression profiles for thousands of individual cells in a single assay. The approach is an unbiased way of identifying cell types in a mixed population of cells, both known and uncharacterised, in addition to the genes that define them. This has seen a surge of interest in dimensionality reduction methods—in particular t-SNE [26] and UMAP [27]—as users seek to optimise the visualisation of results (Fig 4A and 4B). These methods are constrained by issues associated with representing the underlying data structure, i.e. the relationships between cell groupings. When 10s of thousands of cells are analysed, a 2D plot space is limiting.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. Analysis of cell and gene associations in scRNA-Seq data.

The structure of scRNA-Seq data is commonly represented using approaches such as (A) t-SNE and (B) UMAP as shown here for immune cells derived from the Tabula Muris dataset. However, the distance between data points and groups of data points is difficult to interpret. (C) Graphia enables the construction of cell-to-cell networks built on a similarity parameter. Here, the 48 most significant PCA values for each cell were first calculated and this PCA profile used to construct a correlation network. The plot bottom left of C, shows the PCA profiles of cells in the two largest cell clusters. To better show graph structure, a k-NN (k = 10) transformation was applied and outlier cells removed (r < 0.85 and node degree < 10, nodes coloured white). The graph comprises of 12,498 nodes (cells) and 143k edges. Cell clusters have been annotated as the cell types defined by the authors. (D) Shows a gene correlation network generated from these data by first calculating the average expression of genes within cell clusters and then calculating a correlation matrix from these values. (E) Plots show the average expression profile (y-axis) of a selection of gene-clusters across the aggregated cell-clusters (x-axis). The label gives the cluster number, e.g., C1, the number of genes within the cluster (966) and the association of the genes with a given biology or cell type.

https://doi.org/10.1371/journal.pcbi.1010310.g004

An alternative approach is to treat single cell data as a graph, where nodes represent cells or genes, and edges the similarity between them. There are a number of measures that may be used to calculate the distance between cells or genes and here we discuss our currently favoured approach. The Tabula Muris dataset [28] includes a scRNA-Seq data from 20 different mouse tissues. For the purpose of this study we selected only data from tissue immune cells, as annotated by the authors. Preprocessing and quality control was performed as per the Tabula Muris pipeline (https://github.com/czbiohub/tabula-muris), producing a normalized dataset of 14,466 cells from 12 tissues. Principal component analysis (PCA) was conducted to reduce the gene profile of a cell to principal components (PC), with the 48 most significant PCs being considered (adj. P-value <0.05), based on Jackstraw permutations [29]. This file (cells as rows, PCs as columns) was loaded into Graphia and a graph was generated based on the Pearson correlation coefficient between the PC profile of cells. An initial network was generated applying the k-NN algorithm (k = 15), and outlier cells were removed. Outliers were defined as cells with poor correlations and connectivity with other cells of the same cell type (r < 0.85 and node degree < 10). Applying these filters produced a network that better separated known cell-types. Subsequently, the filtered 12,498 cell-to-cell network was clustered using the Louvain clustering algorithm [24] with a granularity of 0.8, identifying 36 cell clusters (Fig 4C).

Following the identification of cell clusters, it is generally of interest to identify gene markers and expression modules associated with different cell types. Due to the inherent noise in scRNA-Seq data, it is of limited use to construct gene networks based on the correlation between expression values of individual cells. One alternative is to construct gene correlation networks based on an aggregated gene expression value across cells from each cluster, i.e. the similarity between clusters of cells instead of individual cells. Accordingly, a matrix of the average gene expression values across the clusters defined above was loaded into Graphia. A network of strongly correlated genes (r > 0.85) was generated and gene coexpression modules identified using the Markov clustering (MCL) algorithm [23] with a granularity setting of 1.7 (Fig 4D). Graphia enables a dynamic and rapid exploration of these clusters, allowing a user to understand where a given cluster sits within the context of the entire graph, the identity of genes present within a given cluster and the profile of all or some of those genes across samples, in this case cell clusters (Fig 4E).

Case Study 3: Exploration of bacterial pangenome structure

Whole genome sequencing is now routine in many fields. One common use is in the characterisation of microbial species and public databases already hold tens of thousands of genome sequences for the best studied organisms. Graph-based methods and tools that support the visualisation and analysis of such data are well established. For instance, Bandage is a software package now widely used to visualise de novo assembly graphs of bacterial genomes [8]. PPanGGOLiN [30] and Panaroo [31] are graph-based pangenome clustering tools for the analysis of genomic diversity within a bacterial species (i.e. its pangenome), which can then be used to statistically classify genes according to their occurrence in the genomes.

Comparative analyses of bacterial sequences have revealed a high degree of genetic diversity between isolates of the same organism, leading to the concept of “core” genes present in all isolates and “accessory” genes present only in some isolates. The distribution and organisation of accessory genes has a significant impact on an organism’s ability to adapt to different hosts and niches, virulence and drugs. Fig 5 shows a graph generated from a previously published dataset of whole genome sequencing data of 778 Staphylococcus aureus isolates [32]. Genomes were annotated using Prokka (v1.13) [33] and their pangenome defined using the analysis pipeline PIRATE (v1.0.3, default settings used, 95% sequence similarity threshold) [34]. In this visualisation, nodes represent individual genes/gene variants and edges the syntenic relationships between them. When the visualisation is enhanced through making node/edge size and colour proportional to the number of genomes in which a gene is present or two genes are syntenic across the dataset, core regions of the genome become easy to identify (Fig 5A). These can also be collapsed to single nodes to simplify the graph. Similarly, areas of high variability within the pangenome are obvious (Fig 5B). Graphia can be used to identify specific nodes, for example the path of a single genome (RF122) can be shown in the context of the wider pangenome (Fig 5C) or used to explore smaller local variations, e.g. the quorum sensing locus agrABCD, which shows four variants when defined at this similarity threshold (Fig 5D and 5E). In principle, displaying variation between sequences as a network is applicable to any such data. This has relevance to “pan-reference” genomes for more complex species such as humans, as well as for showing other variation, such as clustering repeat regions [35] or alternative isoforms of transcripts [36].

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 5. Visualisation of the pangenome of Staphylococcus aureus.

(A) The full pangenome of 778 isolates. Nodes represent individual orthologous genes as identified by PIRATE. Node size is determined by the number of genomes in which a gene has been identified. Edges denote where two genes are syntenic, and their thickness is determined by the number of times this syntenic connection is observed across isolates. Syntenic stretches of core genes have been collapsed for clarity using the “Contract Edges” transform, and low confidence nodes and edges (n < 3) have been removed. Coloured by Weighted Louvain Clustering (granularity = 1.0). (B) Highly variable region (boxed area in A) with a high density of “phage-like” genes. Nodes and edges are sized and coloured by frequency. (C) Genes highlighted are all found in single S. aureus isolate, RF122. (D) the agrABCD locus coloured by gene-association clustering. Frequently, an alternative allele is not identified as being the same gene, but their position is strongly indicative of shared function. (E) the same locus coloured by gene identity.

https://doi.org/10.1371/journal.pcbi.1010310.g005

Case Study 4: Analysis of human genome variation

Genome variation also occurs at the level of individual DNA base pairs, so called single nucleotide variants (SNVs). Genotypes can be scored in an individual in terms of their allele dosage, i.e. both being the same as the reference nucleotide (0), heterozygous (1) or homozygous for the variant (2). Calculation of the correlation across a range of these positions in a population of individuals creates a relationship matrix.

Fig 5 shows various views of data from the 1000 genomes project [37]. Here 23,675 SNVs from chromosome 22 were used to generate graphs of the relationships between individuals in this cohort and the SNVs themselves. These data represent the genomes of 2,504 individuals, who were selected to represent 26 distinct ethnic populations from five continents. The average correlation between the SNV profile of individuals is low and the graph shown in Fig 5A is constructed using a threshold of r = 0.238, a value at which most of the genomes in the cohort form one connected component. The k-NN algorithm [20] was applied to reduce the number of edges (from 400.4k to 5996, opening up the local structure) such that only the strongest three relationships between individuals were maintained. The topology of the graph is clearly strongly influenced by the ethnicity of individuals with discrete clusters being observed for all the five continental populations and in some cases individuals from certain countries or ethnicities showing a local grouping within this overall structure. Also visible from the graph are a number of closely related individuals (Fig 6Ai) and a number of instances where an individual does not co-occur with their annotated population, for example there are a number of South Americans in with the Africans, and vice versa (Fig 6B). Transposing the matrix to analyse the similarity between the profile of SNVs across the 2,504 individuals, at the threshold used here (r = 0.75), 11,600 SNVs formed 2,467 separate graph components of more than one node (Fig 6C). After clustering the graph using the Louvain algorithm [22], many of the clusters contain nearby SNVs; likely haplotype blocks, some of which were clearly associated with a given a population. Graph analysis represents an improved approach compared to e.g. PCA plots, visualising genetic associations between individuals and genetic variants.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 6. Analysis of single nucleotide genome variant data.

(A) The graph shown was constructed from data from the 1000 genomes project based on the correlation (r threshold = 0.238) between the allele dosages at 23,675 SNVs from chromosome 22. Nodes represent the 2,504 individuals included in the study and edges the three most significant correlations with their neighbours (k-NN was applied where k = 3). In most cases, individuals’ group with others from the same continent although there are instances where this does not appear to be the case. Visualisation of edge weights (Ai) also highlights cases where individuals would appear to be closely related. (B) Colouring of nodes by the attribute ‘population’ provides a higher resolution to the graph and populations showing a high degree homogeneity have been labelled. (C) Transposing the data upon import demonstrates SNVs whose pattern across the genome covaries. Clustering of these data shows many to represent haplotype blocks and inspection of their profile across genomes, demonstrates some SNV clusters to be associated with a given ethnic grouping, e.g. cluster 3 (Africans) and cluster 14 (East Asia), whilst others little obvious association with ethnicity, e.g. cluster 6. Plots show the average score of SNV’s within a cluster (y-axis, 0,1,2), across the 2,504 individuals ordered by continent and then population (x-axis).

https://doi.org/10.1371/journal.pcbi.1010310.g006