Collection and input of data.

Reconstruction of causal networks using ViralLink requires four separate input datasets (Fig 1): viral protein-human binding protein interactions, a priori human protein-protein interactions (PPIs), a priori human transcription factor (TF)—target gene (TG) interactions and an unnormalised counts matrix from a gene expression experiment. By default, all data except the transcriptomics counts are provided automatically. However alternative input files can be provided if desired.

The default workflow uses SARS-CoV-2 protein-human binding protein interactions obtained from IntAct [22,23]. This data was reformatted to contain one row per molecular interaction with 2 columns of UniProt IDs: SARS-CoV-2 proteins and human binding proteins. Alternative viral-human PPIs can be provided using the same data format. The workflow assumes all viral-human interactions have an inhibiting action on the human protein, unless a third column named “sign” is present in the input file containing “+” for activatory and “-” for inhibitory interactions. In addition, data is provided with the workflow containing the gene names corresponding to each of the SARS-CoV-2 proteins, to enable easy interpretation of the reconstructed networks.

For a priori human interactions, the workflow obtains and uses integrated collections of PPI and TF-TG interactions from OmniPath and DoRothEA, respectively [24,25]. These interactions are obtained using the ‘OmniPathR’ R package [24,26] to download and filter signed and directed interactions. For DoRothEA, only high and medium confidence level interactions are used (confidence scores A-C). In contrast to importing static input files, this script enables the use of up-to-date interaction data. Alternative interaction data can be used with the workflow provided it has the same format: specifically, it must contain source and target UniProt IDs in the columns ‘to’ and ‘from’ and if the transcriptomics data uses gene symbols, the interaction data must additionally contain gene symbols in the columns ‘source_genesymbol’ and ‘target_genesymbol’. Furthermore, the interactions must be directed and signed with the sign of the interaction given in the column ‘consensus_stimulation’ where the value ‘1’ represents a stimulation and anything else represents an inhibition.

The aforementioned a priori interactions are contextualised using transcriptomics data from any study of interest which compares viral infected to uninfected human cells or tissues. Correspondingly, the workflow requires unnormalised counts data from a transcriptomics experiment (containing UniProt or gene symbols as IDs) and a corresponding mapping table which lists the sample IDs (from the headers of the counts table) in the ‘sample_name’ column and the ‘test’ or ‘control’ status of the sample in the ‘condition’ column. This mapping table is used to carry out differential expression of a test condition (e.g. infected) compared to a control condition (e.g. uninfected). An example expression dataset and mapping table are provided with the workflow.

To process the transcriptomics data, the workflow uses ‘DESeq2’ in R to normalise the counts and to carry out differential expression analysis [27]. Any genes passing the log2 fold change and adjusted p value cut-offs, based on the provided parameters (default 1 and 0.05, respectively), are classed as differentially expressed genes (DEGs). Following removal of all genes with count = 0, normalised log2 counts across all samples are fitted to a gaussian kernel [28]. All genes with expression values above mean minus three standard deviations are considered as expressed genes. Subsequently, context-specific human PPI and TF-TG interactions are generated by filtering only interactions where both interacting molecules are expressed.

File paths to all input datasets and associated parameters (such as desired log2 fold change cut off) are specified in the parameters text file which is read in by the workflow.

Network reconstruction.

The reconstructed causal network contains three layers of interactions, which are obtained, by default, from the three a priori interaction resources:

• Viral proteins interacting with human binding partners: from the SARS-CoV-2 collection in the IntAct database [22,23]
• Intermediary signalling protein interactions: from protein-protein interactions (PPIs) of the OmniPath collection [24]
• Transcription factors (TFs) regulating differentially expressed genes: from a transcriptomics dataset of interest and the DoRothEA collection [25]

A list of all TFs targeting the differentially expressed genes are obtained from the context-specific TF-TG interactions. The human binding proteins of viral proteins are connected to the listed TFs through the context-specific human PPIs using a network diffusion approach called Tied Diffusion Through Interacting Events (TieDIE) [29]. As inputs for the TieDIE tool, the following information is used: 1) The signed, directed and expression based filtered PPIs is used as the input network. 2) Human proteins which are interacting partners of the viral proteins are used as the start nodes. The number of viral proteins bound to each of the human proteins are assigned as the weights of the start nodes. 3) The TFs of the DEGs in the dataset are used as the stop nodes. The weights for each of the TFs in the set of stop nodes were calculated using the following formula (Eq 1) which considers both the log2 fold change of the DEGs as well as the sign (i.e. stimulatory or inhibitory) of the relationship between the TF and the DEG.

(1)

After running TieDIE, a custom R script is used to collate all the data into a final viral-initiated intracellular signalling network (causal network), outputting an edge table representation of the network, with a node table containing additional node annotations. Starting with the interactions output from TieDIE, viral protein-human binding protein interactions are added for each of the present human binding proteins. Similarly, TF-TG interactions (where the TG is a DEG) are added for each of the present TFs, creating a full network with three interaction types: SARS-CoV-2 protein-human binding protein, PPI and TF-DEG. All nodes of the network are added to a node table with annotations including heat values (output from TieDIE), Entrez IDs (obtained in R using the ‘org.Hs.eg.db’ package), gene symbols (obtained from UniProt [30]) and log2 fold change values from the differential expression analysis.

Network investigation.

Following reconstruction of the causal network, ViralLink provides functionality to investigate the results using functional analysis, clustering, centrality measures and visualisation.

Centrality measures. To identify key molecules in the reconstructed network ViralLink uses a betweenness centrality measure—calculating the global importance of a node (in this case a protein) based on the number of shortest paths which pass through them when connecting all node pairs in the network [31]. Nodes with high betweenness centrality play a key role in transduction of signals through the network, and here represent proteins with biological importance in the cellular response to viral infection. Betweenness centrality is calculated for each node in the causal network using the R package ‘igraph’ and output as an annotation in the node table [32]. Alternative centrality measures are available using the ‘igraph’ package and can be integrated into the workflow by the user if required.

Functional analysis. To further investigate important cellular functions and signalling pathways directly affected by the virus of interest, ViralLink carries out functional overrepresentation analysis on different parts of the causal network:

1. The DEGs of the network
2. The upstream human proteins (including human binding proteins, intermediary signalling proteins and TFs)
3. Identified clusters (only those with ≥ 15 nodes are investigated)

Functional overrepresentation analysis is carried out in R using packages ‘ClusterProfiler’ (for Gene Ontology annotations [33]) and ‘ReactomePA’ (for Reactome annotations [34–36]. For analysis of the upstream human signalling proteins and analysis of clusters, all proteins in the context-specific human PPI interactions are used as the background. For analysis of the DEGs, all target genes in the context-specific human TF-TG interactions are used as the background. For Gene Ontology (Biological Process) analysis (except when running the ‘compareCluster’ command), the ‘simplify’ command is used (select_fun = min) to remove redundant functions. All functions with q value ≤ 0.05 are considered significantly overrepresented.

Furthermore, ViralLink also implements a network-based pathway enrichment analysis on the upstream human proteins of the network using the ANUBIX (Adaptive NUll distriButIon of X-talk) algorithm [37]. Here, limitations of overrepresentation analysis, such as the assumption of gene independency, are overcome by integrating knowledge of gene/protein associations using networks. Specifically, upstream human proteins of the reconstructed causal network are tested for enrichment of ANUBIX-provided Reactome and KEGG pathways using the context-specific human PPI interactions as the input network. All functions with q value ≤ 0.05 are considered significantly overrepresented.

An additional R script is provided alongside the workflow which creates subnetworks of the causal network based on functions of interest. These function-specific subnetworks highlight how specific signalling pathways in the infected cell reach (and subsequently affect) specific functions of the DEGs. For example, the subnetwork could be created to show how viral proteins can affect different host toll-like receptor pathways, and how these pathways can ultimately affect DEGs associated with interleukins. In this network the DEG nodes would be replaced with nodes representing the interleukin functions (which must be overrepresented based on the functional analysis). This script requires the output files from the functional overrepresentation analysis, the node and edge tables of the causal network and a file of all UniProt IDs associated with all Reactome functions (which is provided with ViralLink, following download from the Reactome website in April 2020). In addition, the script requires a list of overrepresented DEG functions (Reactome) and a list of upstream signalling functions (Reactome) to visualise. The script outputs an edge table, a node table and a Cytoscape file.

Visualisation. Data visualisation is often an important part of biological network interpretation, providing new insights into the data and visually conveying analysis results [38]. As such, ViralLink has the capability to import reconstructed networks into the open-source Cytoscape network visualisation software [39,40]. Specifically, the workflow employs the ‘RCy3’ R package to interact with Cytoscape programmatically, importing the node and edge tables to create network visualisations and saving the data as a ‘.cys’ file. The causal network, the network clusters (where containing ≥ 15 nodes) and the function-specific networks are visualised in this way. If calculated previously, the causal network nodes are coloured based on their betweenness centrality, however further style and layout customisation must be carried out by the user directly based on the data.

Cluster analysis. Clustering algorithms are commonly used in network biology to investigate the complex structure of molecular interaction networks by extracting groups of densely connected molecules [41,42]. Depending on the number of molecules included, a cluster can represent a molecular complex or a group of molecules which function closely with each other. Cluster analysis can identify subsets of a large network with specific functions and indicate molecules that may have functional redundancy with each other—potentially having implications for drug targeting. ViralLink employs the MCODE clustering method to identify groups of densely connected nodes in PPI networks [41]. To carry out this analysis, ViralLink requires the Cytoscape software [39,40], which is controlled programmatically using the R package ‘RCy3’ with the Cytoscape ‘MCODE’ app (v1.6.1) [43]. MCODE is run using default parameters: degree cut off = 2, haircut = TRUE, node score cut off = 0.2, k-core = 2, max depth = 100. This analysis outputs the data as node annotations in the node table, which are used for the functional analysis and visualisation steps of the workflow.

Collection and input of data.

To demonstrate the application of this workflow for the study of SARS-CoV-2, we created intracellular signalling networks of NHBE cells (from Normal Human Bronchial/tracheal Epithelial cell lines) upon infection with SARS-CoV-2 based on data published by Blanco Melo et al. [11] and viral-human binding protein interactions published by Gordon et al. [4].

Network reconstruction.

The resulting causal network contains 804 nodes (molecules) and 5423 interactions (Fig 2A and S1 and S2 Tables and S1 Data).

Download:

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

Fig 2. Causal network of SARS-CoV-2-infected NHBE cells.

A) Signalling flows from left to right: SARS-CoV-2 proteins/protein fragments (red triangles), human binding proteins (yellow parallelograms), intermediary signalling proteins (blue circles), transcription factors (green rectangles) and differentially expressed genes (grey rhombuses). Where a human protein/gene is acting in multiple layers of the network, it is only visualised once based on the following priority: DEGs, binding proteins, TFs, signalling proteins. B) Results of betweenness centrality analysis, which measures the global importance of nodes (molecules) in the network. Nodes coloured based on their betweenness centrality parameter, with the gene names of the 10 highest scoring (most central) nodes overlaid. DEGs have log2 fold change ≥ |0.5| and adjusted p value ≤ 0.05.

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

Network investigation.

Centrality measures. The 10 most central proteins of the reconstructed causal network (based on betweenness centrality) are involved in a wide range of cellular functions (Fig 2B). Taken together these proteins highlight the propensity for SARS-CoV-2 to affect cell proliferation, apoptosis, cell adhesion, exocytosis and proinflammatory immune responses. These functions are influenced through multiple cellular pathways, most notably MAPK/ERK and PI3K/AKT signalling pathways.

Functional analysis. Functional overrepresentation analysis identified an enrichment of interleukin and interferon related functions among the network DEGs, in line with previously published findings (S1 Fig and S2 Data) [12,46,47]. However, the primary value of the ViralLink output, in comparison to a basic evaluation of SARS-CoV-2 induced DEGs, lies within the upstream signalling proteins of the reconstructed causal network. Overrepresented functions and pathways of the upstream signalling proteins (human binding proteins, intermediary signalling proteins and TFs) included innate immunity-related functions, platelet signaling, PI3K/AKT signalling, MAPK activation, estrogen receptor-mediated signalling, senescence and a number of growth factor receptor-associated functions (such as VEGF signalling, receptor tyrosine kinases, stem cell growth factor signalling (SCF-KIT) and neurotrophin receptor signaling) (S2 Data and S3 Data). Of these functions, only innate immunity-related functions and MAPK regulation were also identified in overrepresentation analyses within the original publication of the infected NHBE cells [11]–thus highlighting the added benefit of applying ViralLink to DEG lists.

Further, network-aware KEGG pathway analysis identified a number of viral-associated pathways (such as viral carcinogenesis and Epstein-Barr virus infection), indicating a similarity between cellular responses to SARS-CoV-2 and other viruses and validating the context-specificity of the reconstructed causal network (S3 Data). Moreover, network-aware analysis identified overrepresentation of the Reactome function ‘negative regulators of DDX58/IFIH1 signalling’. Importantly, DDX58/IFIH1 signalling results from cellular sensing of cytoplasmic viral nucleic acids and leads to anti-viral innate immune responses such as the production of type I interferon. Many different viruses have been shown to repress this system through varying molecular mechanisms, including hepatitis C and influenza A and B [48]. Of the 11 proteins in the causal network associated with this function, two are theoretically capable of binding to SARS-CoV-2 proteins based on interactions published by [4] (TANK Binding Kinase 1, TBK1 and NLR family member X1, NLRX1). Although further investigation is required, this finding indicates a possible mechanistic explanation for the poor induction of type I interferon response seen in SARS-CoV-2 infections [11,49,50]. Taken together, we show that ViralLink can highlight additional pathways through which SARS-CoV-2 could be affecting the lung epithelial cells, which cannot be identified by looking at the transcriptomic results in isolation.

Visualisation. Based on functional overrepresentation analysis, we created a function-specific network by subsetting the causal network. This visualisation was used to further explore the mechanisms of how specific signalling pathways are affecting the DEGs (S2A Fig and S4 Data). Specifically, we generated an innate-immunity associated subnetwork containing all upstream human signalling proteins associated with Reactome functions cytokine signalling in immune system, signaling by interleukins and MyD88-independent TLR4 cascade and all overrepresented functions of the DEGs (in place of the DEG nodes). These pathways contain 9/10 of the top betweenness centrality nodes (all except RHOA), evidencing the centrality and importance of the innate immune response to viral infection. Inspecting the TF layer of this immune subnetwork, we find a number of key TFs including STAT proteins (3 and 4), IRF proteins (1 and 5) and NFKB-related proteins (NFKB1, NFKBIA).

Cluster analysis. Finally, we evidenced the application of MCODE clustering analysis using the reconstructed SARS-CoV-2-infected NHBE cell causal network. We identified four clusters containing 15 or more nodes, making up 19% of the network (154/804) (S2B Fig and S2 Table and S1 Data). Assuringly, 9/10 of the top betweenness centrality nodes were included in these four clusters, further confirming the high connectivity and importance of these nodes in the causal network. Functional overrepresentation analysis of the cluster nodes highlighted a functional similarity between all four of the clusters (S2C and S2D Fig and S2 Data). Likely this is due to the high number of inter-cluster molecular interactions and because of the functional similarities between the top central nodes.

Collectively, we show that our systems biology workflow, ViralLink, reconstructs a functionally relevant intracellular signalling network affected by SARS-CoV-2 infection. Investigation of the networks through functional analysis, centrality measures and cluster analysis, combined with network visualisations, enables detailed study of the key proteins and pathways involved in signal transduction.