MetaNetVar: Pipeline for applying network analysis tools for genomic variants analysis

Tools used in the pipeline

As much as possible, the MetaNetVar pipeline uses existing tools for network analysis. We only considered tools that are freely available, with no license restrictions. We describe briefly each tool that is integrated into the MetaNetVar pipeline. Tools vary in scope, and some include additional functions that include network analysis.

FunSeq2 (Version 2.1.2)

FunSeq2 is an existing tool for prioritizing variants using several different approaches, including network-based analysis^4,5. FunSeq2 identifies hub genes and provides the measure of centrality for those hub genes. Inference of the network analysis results requires further processing of the program’s output. We chose to include FunSeq2 in our pipeline because of its capability to identify functionally important, non-coding variants in the context of biological networks.

NetworkX (Version 1.10)

NetworkX is a network analysis framework available in a Python language software package. It allows for “the creation, manipulation, and study of the structure, dynamics, and functions of complex networks”⁶. It contains many standard graph algorithms and accepts and outputs 13 file formats, where nodes can be anything and edges can hold any type of data. NetworkX was used to calculate the degrees and betweenness centrality of nodes (genes) and to create XML format and static PNG figures of subnetworks containing the input genes. The degrees and the betweenness centrality gives you a measurement of how important the gene is in the network. A network analysis framework similar to NetworkX is CytoscapeJS⁷. We chose to include NetworkX in our pipeline because of our experience with Python.

HotNet2 (Version 1.0.1)

HotNet2 is an algorithm for detecting “significantly altered subnetworks in a large gene interaction network”^8–10. The algorithm uses heat diffusion kernel to capture the local topology of the interaction network. The subnetworks in genome scale interactions that have non-random mutations are identified using this approach. The limitation of HotNet2 are the challenges in getting the scripts running straight out of the box, along with the long computational time involved in the preliminary influence matrix creation process. We chose to include HotNet2 in our pipeline because of our experience with Python.

dmGWAS (Version 3.0)

dmGWAS_3.0 is an existing tool for overlaying gene-level summaries of case-control association p-values onto an existing network (in this case, we use the network extracted from GeneMania detailed below) and then identifying subnetworks that are particularly enriched for strong associations using a greedy algorithm^11,12. Unlike the previous version of dmGWAS (2.0), dmGWAS_3.0 also has the ability to incorporate differential gene coexpression data (in other words, the difference in gene co-expression between cases and controls) as weights for the edges in the network, but for the sake of simplicity we did not make use of this new functionality in our pipeline. Due to this choice, we discovered that the dense-module search output (ResList.RData) took the format of the previous version dmGWAS_2.0¹³ and could not be manipulated using the tools referenced in the current documentation. Therefore, we created our own short script to extract out the basic statistics and subnetwork nodes associated with each input gene present in the network (see below: ModuleStrengthSummaryByGene.txt and Top1000ModuleScores.txt). We later discovered that some of the old tools capable of manipulating dmGWAS_2.0 output (ResList.RData) were preserved in the current code package and could be used for further data exploration by a motivated user by loading the output file (ResList.RData) into R and installing the requisite packages (dmGWAS, igraph), although some of the tools did not appear to be fully functional anymore (such as the subnetwork plotting capability in simpleChoose()).

Overall, the primary limitation that we observed for dmGWAS was computing time, so we adapted the existing code to make use of parallel computing using the BiocParallel package in R¹⁴.

To utilize dmGWAS_3.0, it is first necessary to convert the input file containing the case-control association p-values for each SNP to a gene-level summary. Within our pipeline, we complete this conversion using VEGAS, an existing command line (Linux/Unix) based tool recommended within the dmGWAS documentation^15,16. It should be noted that by default, VEGAS uses the HapMap2 CEU (Central Europeans, Utah) population to estimate patterns of linkage disequilibrium for each gene.

VEGAS is written in Perl but also makes use of two R packages (corpcor and mvtnorm) and depends on functions provided by PLINK, a commonly-used whole genome analysis toolset^17,18. The output of VEGAS requires further processing before input into dmGWAS. We found that several of the VEGAS gene-level p-values were rounded to either 0 or 1, which was incompatible with dmGWAS, so we substituted the minimum p-value present in our test file (1e-06) for 0 and 0.999 for 1.

Table 1 provides a summary of the tools used in this pipeline, including notes about their advantages and disadvantages.

FunSeq2

FunSeq2 utilizes multinet¹⁹, which is an integrated network consisting of regulatory interactions from the ENCODE regulatory network^20,21, phosphorylation interactions from the SignaLink database^22–24, protein-protein interactions from BioGRID (release 3.1.83)²⁵, and metabolic interactions from KEGG²⁶. While there are options for users to bring in their own network for use with FunSeq2 analysis, this pipeline prototype uses the pre-packaged multinet.

NetworkX & dmGWAS

The NetworkX and dmGWAS are libraries and do not include particular network data.

We paired GeneMania with NetworkX and dmGWAS. The GeneMania network is a protein-protein interaction network^27,28. Two genes are connected if they are found to interact in a protein-protein interaction study. The network was created from various protein-protein interaction databases, including BioGRID and Pathway Commons^29,30. We used version 2014-10-15 of Homo_sapiens.COMBINED network.

HotNet2

HotNet2 uses mutation data to prioritize subnetworks by identifying significantly mutated subnetworks in genome scale interaction networks. In our pipeline, we have used the 2012 version of HINT (High-quality INTeractomes) a database of high-quality protein-protein interactions³¹.

Example data

As a sample input for our pilot, we searched NCBI dbGaP for a sample study that provided a real world list of variants. We used data from a clinical study of age-related macular degeneration³² with dbGAP identifier phs000182.v3.p1.

As an additional input example, we used data from ClinVar³³. ClinVar is a database of interpretations of clinical significance of variants for reported conditions, hosted by the National Center for Biotechnology Information (NCBI). It includes germline and somatic variants of any size, type, or genomic location with interpretations from several sources (such as clinical testing laboratories, research laboratories, or locus-specific databases). It includes a link of variants to phenotypes. For this example, we identified variations submitted by LabCorp and extracted disease-variant pairs, for diseases with 30+ variants. The example dataset is provided on the MetaNetVar GitHub page.

FunSeq2

The parsed input file required for the FunSeq2 analysis, the PHP script that generates this parsed input file from the original dbGaP association data, and the output files (using default parameters) are provided in the GitHub project page. An example file generated from a filtered list of SNPs from the ClinVar database, for the Cardiomyopathy phenotype is also provided for testing purposes and can be found at http://github.com/NCBI-Hackathons/Network_SNPs/blob/master/test/sample_output/funseq2/cardiomyopathyfunseqoutput.

NetworkX

We took the NetworkX library and created a script that we refer to as SNPsNet. This script generates one output file containing the degrees and betweenness centrality measure of genes that are input into the pipeline, as well as creating two directories (see Figure 2). The two directories contain figures of subnetworks with the input genes and the XML format of the subnetworks. With these results, the user can prioritize the input variants or genes by sorting how important each gene is based on degreeness or centrality, as well as visualizing the subnetwork. Since NetworkX is not primarily a visualization tool, the XML file can be input into several other tools to better visualize the graph.

Figure 2. NetworkX outputs a file containing the degreeness and centrality of each gene, as well as two directories containing subnetwork graph figures for each input gene (.png) and its XML format.

HotNet2

The influence matrix for HINT was pre-computed and then used in the current version of our pipeline. Influence matrix creation is a one-time process for a given network and, if required, advanced users may use custom influence matrices with MetaNetVar by modifying the path to the input influence matrix file and corresponding gene index file. For evaluation of MetaNetVar, we generated heat scores from a test mutation file. The .json file containing heat scores on each gene, which was used in subsequent steps, may be accessed at https://github.com/NCBI-Hackathons/Network_SNPs/blob/master/heat_score.json.

The final step of weighted graph generation uses the influence matrix for HINT, the HINT index file, and the heat score .json file, to remove edges with weight less than the delta value, and extract the resulting connected components. Two output files were generated: components.txt (available at https://github.com/NCBI-Hackathons/Network_SNPs/blob/master/components.txt) and results.json (available at https://github.com/NCBI-Hackathons/Network_SNPs/blob/master/results.json)

dmGWAS

The sample association file from the age-related macular degeneration dataset (phs000182) was parsed down to a two-column text file containing only SNP identifiers (“rs numbers”) and case-control association p-values. This file was fed into VEGAS and a gene-level summary file was created, which was further parsed into a simple two-column text file containing gene identifiers (gene symbol) and “weight” (an integrated p-value for the gene). Within dmGWAS, this input was overlaid onto a network provided by GeneMania to produce a network of weighted nodes from which particularly “dense” subnetworks are identified (full output: ResList.RData). Finally, our program summarizes the data into two easily navigable tab-delimited .txt files which can be viewed within accessible programs such as Microsoft Excel (ModuleStrengthSummaryByGene.txt, Top1000ModuleScores.txt). Figure 3 and Figure 4 demonstrate example output files.

Figure 3. An example of the first output summary file produced by dmGWAS in our pipeline: ModuleStrengthSummaryByGene.txt.

This file provides the Normalized Module Score for each gene included in the network (“Zn”, where a larger value indicates the gene is more enriched for significant case-control associations), and the gene-level summary case-control association p-value provided by VEGAS. It is ordered by percentile rank to allow comparison across different network analysis programs.

Figure 4. An example of the second output summary file produced by dmGWAS in our pipeline: Top1000ModuleScores.txt.

This second output provides similar information as the first output file, but expands it to include the list of genes (nodes) present in each gene of interest subnetwork. Only subnetwork output for the top 1000 seed genes is provided (as determined by percentile rank).

Limitations

The current version of the pipeline is set to use data from dbGaP and ClinVar out-of-the-box. However, advanced users could tweak the provided scripts to make it run using other input formats. Some of the components of the pipeline use processes that are parallel and compute-intensive in nature. Using the provided working implementation of the pipeline through Amazon Web Services requires some computing skills.