ANIMA: Association network integration for multiscale analysis

Contextual functional interpretation of -omics data derived from clinical samples is a classical and difficult problem in computational systems biology. The measurement of thousands of data points on single samples has become routine but relating ‘big data’ datasets to the complexities of human pathobiology is an area of ongoing research. Complicating this is the fact that many publicly available datasets use bulk transcriptomics data from complex tissues like blood. The most prevalent analytic approaches derive molecular ‘signatures’ of disease states or apply modular analysis frameworks to the data. Here we describe ANIMA (association network integration for multiscale analysis), a network-based data integration method using clinical phenotype and microarray data as inputs. ANIMA is implemented in R and Neo4j and runs in Docker containers. In short, the build algorithm iterates over one or more transcriptomics datasets to generate a large, multipartite association network by executing multiple independent analytic steps (differential expression, deconvolution, modular analysis based on co-expression, pathway analysis) and integrating the results. Once the network is built, it can be queried directly using Cypher (a graph query language), or by custom functions that communicate with the graph database via language-specific APIs. We developed a web application using Shiny, which provides fully interactive, multiscale views of the data. Using our approach, we show that we can reconstruct multiple features of disease states at various scales of organization, from transcript abundance patterns of individual genes through co-expression patterns of groups of genes to patterns of cellular behaviour in whole blood samples, both in single experiments as well in meta-analyses of multiple datasets.

report report report report report

Introduction
A frequent issue with bioinformatic analysis is the following scenario: a given dataset is analysed using various approaches in a linear workflow, in an attempt to extract biological features of interest from the data, often at different scales. For instance, a list of differentially abundant transcripts provides information on the system that differs from a list of co-expressed genes. While both such approaches are valid and provide independently useful information, these linear workflows do not not expose potentially informative relationships between different outputs; instead, multiple individual output items (plots, tables, etc.) are written to the output directory. However, multiple kinds of relationships, or associations, exist between entities of the same or different class in biological systems. For instance, two clinical variables might be correlated; the relationship is then expressed as correlation, its extent as the Pearson correlation coefficient, and the statistical significance of the relationship by a P-value, which may be corrected for multiple testing.
Such relationships can be discovered systematically in the analytic workflow. Typically, the result of such an analysis is written to disk in tabular form. However, a table of correlation statistics contains many non-significant entries, and therefore has a low signal to noise ratio, especially in large datasets. A more efficient approach would be to discard all results that fail to meet a pre-specified statistical cutoff and retain only the significant entries as a monopartite or bipartite network. Indeed, bipartite graphs are the dominant paradigm used to encode and represent information about relationships between entities of different classes. Bipartite graphs in systems biology and medicine have recently been reviewed 1 , illustrating numerous use cases and applications as well as describing their mathematical properties. Once the decision has been made to store analytic results and their relationships as bipartite graphs (or other types of network), the next question is where and how to store the data. Cytoscape 2 is extensively used for visualising network data, mainly in biology. Large networks can be searched and subsetted, and basic network analysis can be performed. However, using Cytoscape for a network data store rather than a visualisation tool becomes cumbersome. Recently, graph databases have emerged as a type of noSQL database. Data in a graph database is stored as nodes and edges (relationships between nodes). Both nodes and edges can be assigned multiple properties with values. One such database is Neo4j, which is freely available. Neo4j databases are queried using Cypher, an intuitive query language, and provide computational access to scripting languages via dedicated APIs. In our view, graph databases are the ideal data store for networktype data.
Reproducibility of research is critical for scientific progress. Bioinformatic analyses can be very complex, but usually the results obtained depend strongly on the methods used, software versions, and even operating system. Typical narrative descriptions of analytic methods provide insufficient information to guarantee reproducibility. Workflow tools like Taverna 3 exist but are more suited to performing tasks that link together functionality offered by services (such as webservices). We prefer scripted workflows, where running and re-running the same script on the same data is guaranteed to produce identical results. However, providing the code used in analysis and the raw data does not guarantee reproducibility, as the computational environment in which the analysis is run can also influence the outcome of computations. Typically, this becomes an issue when the functionality of a software package changes between versions. Therefore, in addition to data and source code, one has to provide the exact configuration and package versions of all software involved in the project. In bioinformatics, this can become daunting very quickly, and leads to the well-known problem of "dependency hell" where not only the software packages need the right version, but also their dependencies. A solution for this is to package the entire computational environment in one or more "software containers". The containerization platform Docker 4 is frequently used to fulfil this function and has enjoyed widespread adoption in reproducible research. This is exemplified by Nextflow 5 , a workflow management system using Docker.
A recent paper presents GeNNet 6 , which describes the rationale for scripted workflows and the use of graph databases in reproducible research. In this paper a scripted workflow in R 7 , use of Neo4j to store data and the use of Docker for reproducibility is described. Other recent efforts in the very wide field of -OMICS data integration, research reproducibility and data integration include the Omics Integrator package 8 which takes a variety of -OMIC data (such as transcriptomic and proteomic) as input and identifies possible underlying molecular pathways using network optimization algorithms, NDEx 9 , an online commons for sharing and searching biological networks.
Here we present ANIMA, Association Network Integration for Multiscale Analysis, a framework for producing and interrogating a multiscale association network, which allows summary and visualization of different, but simultaneously valid views of the state of the immune system under different conditions and at multiple scales. While ANIMA employs key strategies presented in the GeNNet paper, mainly "dockerisation", scripted workflows in R and storage of relationships in a graph database, it differs in detail in implementation as well as complexity.
In particular, we have defined novel functionality that allows the extraction of new information form the network structure in ANIMA, given an increased number of node types and increased complexity of the data model of the Neo4j database. In addition, GeNNet is implemented in a single container containing R and Neo4j, whereas ANIMA uses separate containers for each tool.

Amendments from Version 2
Minor changes were made for this version: 1. Support from the Crick Institute is now acknowledged in the "Grant information" section 2. The Introduction was revised to incorporate additional references to related work, and the logical flow has been improved so that it reads more clearly

REVISED
ANIMA is targeted both at bioinformaticians and computational biologists who wish to reproduce the results presented or analyse their own data, as well as biologists who would interact with the system through the provided graphical user interface.
In its current form, ANIMA is best suited for investigating human immunity. The human immune system can be regarded as a complex adaptive system 10 , even though it is integrated with the more complex system of the whole organism. A true systems view 11 of the immune system needs to account for the various aspects that characterize complex adaptive systems generally. This includes emergence, non-linearity, selforganization, noise, scaling, heterogeneity, a network architecture and preservation of context for individual observations 12 . Whole blood is a "window" into the immune system 13 , allowing a reasonably detailed assessment of the overall state of the immune system based on analysis of mRNA transcript abundance patterns from whole blood samples. Here we use the responses to three common infections (acute HIV infection, malaria and respiratory viral infections), as measured by transcript abundance patterns in whole blood, to demonstrate ANIMA functionality

Methods
Method overview ANIMA generates a multiscale association network (stored as a graph database) from multiple data types (expression data, clinical data and annotation data, e.g. biological pathways databases) by executing a comprehensive analytic workflow, enumerating bipartite graphs from the results, and merging all graphs into a single network. Figure 1 provides a conceptual overview of the multiscale analytics pipeline. As biological systems can be understood at multiple scales, our approach aims to integrate information across multiple scales which range from single genes to clinical phenotypes ( Figure 1A). This integration is summarised in the core data model for ANIMA ( Figure 1B) which represents relationships between various outputs of the analytic pipeline. Multiple bipartite graphs are generated from the outputs, and finally merged into a single data structure (Figure 1 C-E), which is then accessed to gain novel information about the system. The various steps as implemented in the ANIMA_build script are illustrated in Figure 1F.

Data ingestion and preparation
The first step in constructing an ANIMA database consists of accessing raw data, which consists of non-normalised Illumina BeadArray expression microarray data and clinical/phenotype data. A script (ANIMA_data.R) imports the data to the R workspace as well as the clinical data and creates a LumiBatch object of the imported data. Based on the experimental design, the data may be subsetted at this point. One or more LumiBatch objects are then saved to the project output folder for later re-use.

Scripted workflow execution
The second step in constructing an ANIMA database consists of iterating over all datasets included in the analysis (each saved as an RData file), sequentially performing thirteen analytic tasks (A1-A13 in Table 1 and Figure S7; see Supplementary Data for details). Array normalization, probe filtering and differential expression analysis are frequently performed on transcriptomics datasets, but the other analyses tend to be performed in isolation. To our knowledge, this is the first transcriptomics workflow to combine WGCNA and Chaussabel modular analyses and integrate these with cellular-level approaches. The chosen approaches were selected based on their perceived popularity in the systems immunology literature, and selected annotation sources are biased towards immunology and celltype specific data sets. This does not exclude the possibility of other approaches or annotation sets to be substituted in their place if the focus of investigation concerns other tissues or species of interest.

Enumeration of bipartite graphs
Following the transcriptional analytic tasks, the algorithm constructs twenty-nine bipartite networks from a multi-scale analytic pipeline ( Figure 1A, C), combining the output, and the relationships between different classes of output approaches (Table 1). Each network contains the associations between two distinct data types ( Figure 1B, D, Table 2, Table 3). The final association network is a result of graph union, merging all networks (Table 3) on shared node types while retaining all edges ( Figure 1E). This results in a data structure that exposes the relationships between modular gene expression and higher-level phenomena, while retaining key probe-and gene-level information. Three types of association (with their respective association indices) are utilised, resulting in three distinct edge types in the final data structure (See Figure 1B, 1E, Table 2 and supplementary methods; Supplementary File 1). Weighted gene co-expression network analysis 14 (WGCNA) is the core analytic method in ANIMA as this is used to discover biological processes in the system of interest. Supplementary Figure S7 indicates the points where the bipartite graphs are enumerated.
Network construction relies on three principles. Firstly, mathematical operations on data are performed, independent of prior knowledge (e.g. WGCNA networks). This aspect of the approach is completely unsupervised. The second process involves the testing of hypotheses, (e.g. determination of differential transcript abundance and differential co-expression of transcripts, requiring knowledge of phenotype/trait classes for the samples). Finally, the results of the first two processes are integrated in various ways with prior biological knowledge. The result is a collection of statistically robust analytic results and various associations between them and known biology, in the form of a multiple bipartite graphs.

ANIMA database
The twenty nine graphs share node types between them, as indicated in Supplemental Figure S7b, which describes the "data model" for the graph database. As each bipartite graph is enumerated, it is added to the Neo4j database instantiated at the start of the build process. Nodes and relationships are added using the "MERGE" command, ensuring that nodes are not added in duplicate. At the end of the build script run, a large graph has been generated, stored as a Neo4j graph database. This  is referred to as the ANIMA database, which makes the overall network structure as well as the individual nodes and edges accessible for further analysis (see supplementary methods; Supplementary File 1).

Accessing ANIMA
After network construction, information in the graph is accessible and utilized to expose new information not present in any of the individual steps. The key to making the ANIMA *a minimum of three and a maximum of ten results were included in the ANIMA database ** a minimum of four and a maximum of ten results were included in the ANIMA database *** results were reviewed algorithmically; in the case of no hits for enrichment at the specified cutoff, less stringent criteria were applied (See Supplementary information). In all cases, the corrected P-value was included as a property of the relationship, allowing filtering of the database entries.
X the PROBE-PROBE network is a monopartite correlation network and not a bipartite graph Abbreviations: BH, Benjamini-Hochberg multiple testing correction database useful lies in the use of functions and web applications (see supplementary methods; Supplementary File 1) that query this large multipartite graph and return visualization of relationships or tables of nodes and/or links (associations between nodes). This has been implemented in several R functions which underpin a Shiny web application called ANIMA REGO.
Constructing an ANIMA database from user data This paper and associated files provided as supplementary data allow the reconstruction of the particular ANIMA database presented here. However, it warrants emphasis that we are describing a method here, and that the provided datasets can be exchanged for a user's own data. The scripts at present only support Illumina microarray data, although future work aims to add support for other microarray platforms and RNAseq data.
In order to construct a new ANIMA database, the following is required in the source_data folder: 1. Non-normalised expression data (.txt format) 2. Clinical/phenotype data (csv format) 3. A file named "matrixPD" in csv format or similar that contains the names of the datasets to be analysed, as well as the classes of each of the variables included in the phenotype data, distinguishing numeric and categorical data types, and flagging certain variables as identifiers.
4. A file named "questions.R" which is an R script specifying a list which contains, for each dataset, multiple variable assignments required in the running of the analysis. Essentially, this file encapsulates the information required to execute the experimental design.

5.
A file named "setlist.R", another R script, providing high level project information, the url of the graph database, the various studies the datasets are drawn from, as well as the identifiers of data sets that contain matched samples, which enable additional analytic functionality.
With the required files in place, the scripts ANIMA_data. R and ANIMA_build.R are run in sequence, and the end state should be a new ANIMA database based on the user data. The supplementary information provides full replication instructions.

Using ANIMA
The novelty and value in ANIMA lies in its three broad approaches to data interpretation: Firstly, it allows detailed, multiscale investigation of a single dataset with a focused research question where two phenotype classes are compared (multiscale class comparison). Secondly, as datasets are stratified by default using two variables (typically one identifying two disease classes, and the other either a potential confounder (e.g. sex), or a second biologic variable of interest, (like co-infection with a second pathogen), we can examine the interaction of the second variable with the first using a factorial study design (factorial analysis). Finally, multiple datasets, and multiple conditions can be meaningfully compared and contrasted to identify similarities and differences (meta-analysis), both at cellular and modular level.

Application areas
In its current state, the ANIMA approach is optimised to address questions in immunology, but other application areas that may be tissue type specific (neuroscience, oncology) are also possible, but will require addition of annotation libraries and lists of cell-type specific genes.

Results
Example study: data sets We analysed three publicly available microarray datasets using the ANIMA pipeline and toolset (ArrayExpress/GEO identifiers: E-GEOD-29429, E-GEOD-34404/GSE34404, E-GEOD-68310/GSE68310). The first compares a cohort of subjects with acute HIV infection to healthy controls 25 , the second compares symptomatic malaria to asymptomatic controls in children from a malaria-endemic region 26 , and the last compares the host response in early symptomatic viral respiratory infections 27 . Where the datasets contained samples from multiple timepoints, we restricted the analysis to healthy controls and the first disease timepoint. The respiratory virus infection dataset contained samples from subjects with infections other than influenza or rhinovirus; we excluded those from this analysis. Figure S1 (Supplementary File 1) shows the experimental design for the factorial analysis in limma for each of the datasets, together with the numbers of samples in each of the individual groups. Each of the datasets also included clinical data, which was integrated in the analysis.

The ANIMA network
The result of the script that builds the ANIMA network is a large, mostly connected, multipartite graph. Figure S2 (Supplementary File 1) shows a subgraph of this network (for dataset HIVsetB, edge 5).
Searching the ANIMA network using network paths and filtering on node and edge properties In the most direct approach, the large ANIMA graph in the Neo4j database can be searched directly from a web browser using the Cypher Query Language (CQL).
Consider the following query: The query language combines commands and functions (e.g. "MATCH", "WHERE") with a typographic representation of nodes and relationships, where a node is indicated with round brackets ( ) and a relationship is indicated by dashes either side of a set of square brackets -[ ]-. Taken together, the query above conducts a search (or graph traversal) for a specific graph pattern matching MATCH WGCNA modules for the HIV positive vs control comparison in the HIVsetB dataset (n:wgcna {square:'HIVsetB',edge:5}) whose module eigengene (ME) is strongly correlated with any clinical variable (Pearson R > .6) (ph:pheno)-[r1]-; WHERE r1.weight > 0.6, and also returning the genes mapping to probes within those modules -[r2]-(p:PROBE) that have high levels of over-abundance in HIV-1 infected individuals AND p.logfc > 2.
In addition to the standard web browser interface for Neo4j, into which the above query can be directly entered and returned in the browser window (Figure 2A), we also provide a function in R (igraph_plotter) that returns the network found, and plots this within R ( Figure 2B), exports the network as node and edge lists for import in other software like Cytoscape ( Figure 2C), or returns the result as an igraph 23 object for further manipulation within R. (A file ANIMA_styles.xml is included in the common folder supplied with the source code; this is a Cytoscape stylesheet that reproduces the colouring shown in the Figure 2C) More sophisticated methods of using ANIMA are described next. These utilise both bespoke R functions as well dynamic interactivity provided in the Shiny web interface (see below).
Approach 1: Multiscale class comparison Accessing individual transcript abundance levels in multiple conditions. It is useful to view transcript abundance patterns of specified probesets, for instance to compare microarray data to RT-PCR validation data, or to investigate the behaviour of groups of biologically related genes in various conditions. We provide an interface for this in ANIMA. The user can submit a search string (in the form of a regular expression) containing gene names (HUGO Gene Nomenclature Committee (HGNC) symbols), and box-and whisker plots for the results are returned. In the original paper on acute HIV infection the authors discuss a gene set of six conserved genes that appear at multiple timepoints in an inferred regulatory network of viral set point 25 . We show the normalized expression data stratified by HIV status and sex in the two datasets included in the HIV analysis for these six genes (Figure 3). In the paper on acute viral respiratory infection, IFI27 and PI3 are identified to differ between acute influenza A and human rhinovirus infections. In influenza, IFI27 is upregulated and PI3 downregulated relative to human rhinovirus. The malaria study replicated prior knowledge of differential transcript abundance for C1QB, MMP9, C3AR1, IL18R and HMOX1; we show similar results for these transcripts. Supplementary Table 1 lists the results of differential expression analysis for the above transcripts in the three conditions, providing validation of data at individual transcript level.

Functional annotation of WGCNA modules. An important
question is what do WGCNA modules represent, given that these are groups of genes that co-vary across samples, but that can differ dramatically in size. Instead of searching only for evidence of co-regulation of expression by transcription factors, one should consider other causes of this co-variance. We propose the notion of WGCNA modules representing biological processes that may be regulated at different scales. For instance, a group of transcripts will co-vary across samples if these transcripts are expressed (predominantly) in a single cell type, and the proportions of this cell type varies between samples (Figure 4). Another group of transcripts may be expressed in multiple cell types, but represents a concerted transcriptional program executed in response to a specific stimulus, such as interferon-alpha stimulating the expression of a specific group of interferon-regulated genes (see below). Generally, we determine the function of WGCNA modules based on all statistically significant associations with pathways and cells.

Relationship of modules to clinical variables.
We performed Pearson correlation of WGCNA module eigengenes (ME) and clinical variables (described in supplementary methods). Figure 5 shows correlation of the pink ME with age and CD4 count in the HIVsetA data set. Clearly, the pink module is significantly associated with CD4 count, an association that is independent of age. Further investigation shows that this module is linked to interferon signalling, and probably expressed in a variety of cells. This agrees with our understanding of acute HIV infection, which is associated with a robust type-I interferon response and an acute drop in CD4 count 28 .
Investigating the structure of WGCNA modules. WGCNA modules are groups of co-expressed transcripts. Demonstrating the extent and direction of correlation of their constituent probes, and their relationships to biological pathways requires sophisticated visualization ( Figure 6). We showcase the example using the HIVsetB dataset and focus on two modules. The turquoise module shows the coordinated action of genes involved in cell division, relating to the cell proliferation in the lymphoid compartment. The yellow module is clearly related to interferon signalling; of interest here is that within this module there is a group of transcripts highly correlated with each other, all related to interferon signalling, suggesting that these transcripts may all be downstream of a single regulatory factor. Also clear from the figure is that this module is not limited to a single cell type, but rather to innate immune cells in general.

Relationships between module-based approaches.
Both WGCNA and the modular approach pioneered by Chaussabel, et al. 17 rely on clustering of similarity matrices to derive modules. A key difference in these methods in ANIMA is that the Chaussabel modules are pre-defined, whereas the WGCNA modules are derived from the transcriptional data under study. It is therefore interesting to discover the relationships between these two approaches. Figure 7A shows the bipartite network of WGCNA and Chaussabel modules derived from the HIVsetA dataset, as defined by the hypergeometric index. Figures 7B and 7C show the two projections of the bipartite graph. The hypergeometric index is not the only way associations between the two module types can be demonstrated; for instance, a more indirect association can be inferred when a WGCNA and Chaussabel module map to the same biological pathway. It is clear from the plots that only a subset of the list of Chaussabel modules associates with WGCNA modules in any given condition.
Deconvolution. An important question in analysing transcription data from complex tissues is whether differences in transcript abundance are attributable to transcriptional regulation in one or more cell types, or to changes in the composition of the overall leukocyte populations.    Figure S1; Supplementary File 1). Shown are WGCNA modules whose expression correlates with specific cell-type proportions (dark green, edges annotated with Pearson correlation coefficient R) and that are enriched for the genes specific to that cell type (medium green, suffixes xp_1-3 indicate the respective gene list on which the cell assignments were based, see Supplementary methods). The classes of cells are indicated in light green. The modules are annotated with coloured rings representing the difference in median eigengene values between cases and controls (diffME, see Supplementary methods); blue indicates modules which are under-expressed, and red indicates modules that are over-expressed in cases relative to controls. WGCNA module names are (arbitrarily) based on colours as per the convention of the WGCNA package, and modules were not renamed manually.

Virtual cells: an estimate of functional phenotypes of different immune cell types.
Given that the immune response is mediated by different types of cell, we attempted to re-create "virtual cells" based on the assumption the genes that co-vary in terms of transcript abundance across samples with that of celltype specific genes are expressed in that particular cell type. With these relationships, we generated virtual cells (probe coexpression matrices annotated with biological pathways and probe differential expression data). Figure S4 (Supplementary File 1) shows two cells (B-cells and neutrophils) in acute HIV infection; both characterized by interferon signalling. Given the "virtual cells" and their functions we can compare pathway-level transcript abundance in different cell types by creating a matrix of cell types and pathways, with each entry representing the "pathway activity" for a given cell and pathway combination. To illustrate this, we show replication of the finding of NK-cell activation in acute influenza infection (respInf dataset, Figure 8 A,B), and in addition provide more detail on which pathways are probably up-or downregulated in these and multiple other cells.
Approach 2: Factorial analysis Modules driven by sample class or sex. Our main interest in investigating transcriptomic datasets is to identify molecular and cellular processes that drive, or at least are associated with, specific phenotypic traits of the samples. Therefore, the experimental design for the differential abundance analysis and the probe filtering steps prior to WGCNA module discovery were both designed to highlight processes associated with two factors: disease class and sex; the former because this is the basis of the research question for the three studies, and the latter because sex has a distinct influence on immune system function 30 . We developed a novel approach to quantify the association of the WGCNA module expression with these two factors (See supplementary methods; Supplementary File 1). Figure 9 shows that, in the case of acute HIV infection, most modules are associated with the disease process (HIVsetB dataset), but that one module (purple) is strongly associated with sex. Figure S5 (Supplementary File 1) shows the module eigengenes for the yellow and purple modules, demonstrating differential class associations for WGCNA modules; this finding would have been missed had the data not been stratified by both HIV infection status and sex. Table S4 shows the module statistics for this dataset.
Additional modules were identified that associated with neither sample class nor sex. These represent biologic processes that manifest in heterogeneity of the sampled population. Figure S6 (Supplementary File 1) plots the study subjects in the HIV infection data set (HIVsetA) on two axes represented by two module eigengenes.
Approach 3: Meta-analysis of multiple datasets Module-level meta-analysis. Meta-analysis of multiple related expression datasets can lead to insights not available from analysis of any single datasets, and can highlight common where applicable. Note that the same coefficient is obtained for CD4 count as in panel A. Legends are shown for vertex type and diffME (a measure of differential co-expression (see Supplementary methods), i.e. the extent that the module eigengene median varies between two classes). Abbreviations: diffME, differential module eigengene.
patterns of transcript abundance across different conditions, or meaningful differences across highly common conditions. We implemented the approach pioneered 20 and refined 21 by Chaussabel, et al. to perform modular transcriptional repertoire analysis 31 on the six datasets. This approach is particularly suited to meta-analysis, as the composition of the modules is always identical. Figure 10 shows modular patterns for the six datasets in a clustered heatmap as well as the subset of modules with similar expression patterns across the six data sets, demonstrating universal patterns in the immune response to infection. Supplementary Table S3 (Supplementary File 2) lists module functions based on all significant enrichment associations for the modules in Figure 10B. Universal upregulation of interferon-related modules particularly stand out,

Meta-analysis at the cell-type level.
A second approach to meta-analysis was implemented using virtual cells based on WGCNA modules. Here we compared the pathway scores in a single or several cell types across multiple conditions. For instance, comparing the CD8 T-cell response in acute HIV, acute viral respiratory infection, and symptomatic malaria, we find that proliferation of activated CD8+ T cells characterizes acute HIV infection, and to a lesser extent symptomatic malaria infection. In contrast, there is a suppression of CD8 T cell activity in blood in other conditions, due in part to a reduction in CD8+ T-cell proportions in whole blood ( Figure 10C).

Runtime statistics and output
The build script produces a file called session_output.txt. This file captures various runtime statistics, as well as some textual output. This specific ANIMA database was created on a 2013 Mac Pro, 8 cores, 64 GB RAM, with 7 cores and 55 GB RAM allocated to Docker. The build phase took 34.3 hours to complete. Once the database is built, the various components of the Shiny web interface compute and render within seconds, depending on resources available.

Discussion
Systems immunology aims to understand the complex web of relationships between immune system components (cells, cytokines, effector molecules and other mediators) in immunemediated disease states. Much progress has been made in single-cell techniques that yield large amounts of information, e.g. single-cell RNAseq and mass cytometry. These approaches however are expensive. In contrast, RNAseq or microarray analysis of complex tissue samples (like blood) in principle contain information on the transcriptomic state of all cells present in the sample and is thus an unbiased approach. The main difficulty with this lies in the interpretation of the data, and in many cases a complexity-reducing analysis approach is employed, where the focus is placed on differentially expressed genes. Other approaches based on co-expression analysis often fail to explain the drivers of the co-expression patterns.
We demonstrate, using a novel method of aggregating information obtained from clinical and microarray data, an ability to reconstruct many aspects of the immune response, and to discuss this not in the language of probes, genes and signatures, but rather as coordinated biological processes and the cellular context for these processes, allowing the generation of hypotheses at multiple scales. Our multiscale class comparison approach can be used to validate findings from individual papers, for instance the intense NK-cell activation described in influenza 27 . In contrast to other approaches, e.g. PARADIGM 32 , which rank pathways in specified conditions, we utilise the multiscale association information to infer states of immune cells (virtual cells) and providing an interface to compare such cell states within and between datasets. Pathway and cell-activity ranking are provided as a summary of multiple cell states.
Using our factorial approach, we can begin to dissect inter individual heterogeneity in transcriptional patterns from transcriptional patterns that are in a causal relationship with defined factors. Application of the two meta-analysis approaches allows comparison of arbitrary datasets to detect similarities and differences at modular and cellular levels. An interesting and somewhat unexpected finding is that acute symptomatic malaria and acute respiratory viral illnesses are more similar to each other than to acute HIV infection, another viral illness. Despite these differences, we demonstrate that, at least for these three rather different infections, a broadly similar pattern of transcriptional module activity can be described.
In summary, ANIMA is both a robust implementation of various well-regarded analytic paradigms in microarray analysis, as well as a framework for integrating these various methods to expose relationships at multiple scales and render these computationally accessible.

Future work
ANIMA is an open-ended project. Future work may include a demonstration that the associations in the final network are robust to choice of package. This requires a technical comparison study incorporating multiple packages with the same goals. A priority of ongoing development is improving the ease with which users can import their own data, as well as adding support for more microarray platforms and RNAseq data.
Customisation of ANIMA at this stage requires modification of the source code to explicitly add additional functionality.
With modularization of code, we hope to make this process easier. The scope of ANIMA can be widened by including support for more tissue types and organisms. Integration of protein-protein interaction data using existing databases is planned. Finally, we envision a repository where completed ANIMA projects may be archived as well as hosted. Deffur and colleagues describes a transcriptomics data integration and analysis pipeline named ANIMA. ANIMA starts from raw microarray data sets, processes them using R and then applies WGCNA and another previously published module discovery method. The resulting networks (covering different relationships extracted from the data sets such as phenotype-differentially expressed genes, gene-gene co-expression, etc.) are stored using Neo4j graph database, allowing fast querying of relationships that satisfy certain criteria. In addition to the pipeline, the authors build a web application using Shiny. They also showcase the use of the pipeline by analyzing three publicly available data sets. The article explains the technical details of the pipeline clearly. Although ANIMA currently supports only Illumina BeadArray other platforms such as Affymetrix and sequencing data is planned to included in the future. Indeed, addition of other platforms and RNAseq data analysis capabilities would be a important improvement.

Software availability
The intro falls short on existing related work, several recent efforts aim to integrate, analyze and 1.

5.
The intro falls short on existing related work, several recent efforts aim to integrate, analyze and store omics data as well as provide a reproducible work flow such as Omics Integrator (Tuncbag et al I find using color names to represent modules in Figure 4 highly confusing, especially given that blue is also used to denote underexpression. They can be renamed and/or described better in the legend. Minor: The focus on immune system in the beginning of the intro constitutes a slight disconnection with respect to the abstract / rest of the text. I understand the focus due to the case studies but still reads isolated from the main theme of the paper, which is introducing ANIMA as a framework. The authors might as well redirect the focus to inference of cell state for immune response but then they would need to provide a more stringent validation of the obtained results in terms of the biological findings.
The use of subsections in the introduction feels rather unconventional for scientific writing and somewhat confusing given that the current approaches only mention one existing approach followed by the proposed approach.
Question 2 in problem definition fails to pose a general question, going to the technical details before having explained the graph under concern, the triples, P-value correction etc.
In fact, to me it seems that the authors describe the problem in the latter "rationale" sections, where they highlight the problem on the storage of research-related results and network data as well as the difficulty of reproducibility due to software / platform version changes.
"high order, high complexity and low dimensional state space vectors from high dimensional low order, low complexity data (i.e. non-normalised microarray data and clinical/sample phenotype" not clear what the authors mean by the order and complexity of the data (I see later in the figure 1 that it is biological complexity). Moreover, this sentence is very long and can use a revision.
Some references are not in right order (e.g., Chaussabel refered as 14 while being 17).

suffixes => suffices
Is the rationale for developing the new method (or application) clearly explained? Many thanks for the thorough review and comments.
The introduction has been reworked to address most of the comments made in the referee report. It now reads more clearly. Suggested references have been incorporated. With respect to comment (4), the legend of figure 3 addresses the question. The introduction now reduces the focus on the immune system, first discussing the ideas underlying the ANIMA framework, before showcasing ANIMA by example data derived from immune responses to infectious diseases.

Subsections in the introduction have been removed.
Vague allusions to the theory of complex systems have been either removed or made clearer.
While I agree that colour names for WGCNA modules can be confusing, I avoided relabelling the modules manually to remain true to the original output. Replacing colours with numbers would not offer a significant benefit, and coming up with function-related names would be to some extent subjective. Indeed, the annotation of the WGCNA modules suggests that assigning singular descriptions to them is not straightforward.
No competing interests were disclosed. Many thanks for reviewing the second version of the paper. Table 3 has been updated to indicate the "Chaussabel modules" are equivalent to the "baylor" datatype. In future versions of ANIMA we will change the datatype in the source code; as an aside it would be ideal to avoid eponyms for this, as in the case of WGCNA.
We did not include a table describing each individual datatype in the main text, but it is available in the Supplementary information as Table S5, and detailed descriptions are given on page 9-10 of the Supplementary file. Deffur et al. present ANIMA, a network-based method for integrating gene expression (microarray), annotation data, and clinical data. ANIMA utilizes a merged graph database to represent and allow exploration of associations/relationships between biological processes and cellular contexts at multiple scales. The authors provide software (R, Docker) and extensive documentation for building and interacting with ANIMA networks for several example infectious disease datasets; they also present a web interface (R/Shiny application) for interactive visualization of results. The example datasets are used to demonstrate three central approaches for interpreting data with ANIMA, highlighting selected results and features of the method and tools.
I would be interested in using this method for some of my own projects related to systems immunology and immuno-oncology, and the web tool seems like it would provide a nice tool for biologists who wish to explore a particular dataset. As a "framework," ANIMA would be even more appealing to computational biologists as a more modular/extensible system -but still offers a healthy breadth of analyses in its current state. While this manuscript is technically sound, there are issues with structure and organization that I would recommend addressing. In particular, the authors should (1) more clearly describe the rationale for their methods; and (2) present method details more prominently in the main body of the manuscript.

Major comments/critiques
Regarding rationale:

Minor/discretionary comments
For the Cypher query example, a corresponding "translation" of the syntax would be helpfuleither as something more SQL-like or a verbose, literal explanation of what the query is doing. There's a typo in Figure S10 ("F[i]lesystem"). In the "Discussion" section, the reference to "selection of one or a few cell types" as a limitation of single-cell RNAseq is increasingly less true with technologies like 10x -might want to revise. For differential transcript abundance prior to WGCNA, genes are "ranked" and the top 4000 are selected; similarly, the 500 are used in pathway enrichment. For Chaussabel modules, was any cutoff used to determine which genes were up-or down-regulated? Authors should clarify that for this particular iteration (due to the choices of cell types from Abbas et al. and the Chaussabel modules), the method is most applicable to datasets with an immune component or context.

Is the rationale for developing the new method (or application) clearly explained? Partly
Is the description of the method technically sound? 1.

Are sufficient details provided to allow replication of the method development and its use by others? Yes
If any results are presented, are all the source data underlying the results available to ensure full reproducibility? Yes Are the conclusions about the method and its performance adequately supported by the findings presented in the article? Partly No competing interests were disclosed.

Competing Interests:
Referee Expertise: Computational biology, systems biology, systems immunology, bioinformatics 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.
Author Response 20 May 2018 , University of Cape Town, South Africa

Armin Deffur
Many thanks for the thorough and comprehensive review. The following changes were made in order to address the major concerns: We discuss the key distinction of ANIMA compared to other methods as the ability to generate virtual cells from the expression data and annotation libraries; the specific example referenced (PARADIGM) outputs rankings of pathways We thank the reviewer for pointing out the Costa et al paper, which was published after our own development of ANIMA. While both approaches rely on scripted workflows, graph databases and reproducible computational environments, our approach is sufficiently novel in its ability to generate new information from the multipartite graph. A more detailed introduction is now provided, providing extended justification for the approaches used The Methods section was reworked by moving some material from the supplementary information to the main manuscript, and adding a new Table 3, which lists the individual bipartite graphs. Figure S7 is now cross-referenced to Table 1, and the various workflow components indicated in the Table and Figure. A section has been added for the requirements (hardware, data format, etc.) of using ANIMA for a user's own data, and the outline for the procedure is listed. A simpler version of Figure S9 has been added as a subfigure to Figure 1 in the main manuscript.
We feel that the Shiny web application is key to delivering the functionality of ANIMA from a 6.

7.
From a computational perspective, the motivation is clear. But from the perspective of a potential biologist why would they use this over something else.
What tools does this compare to and how is your method an improvement?
Is the rationale for developing the new method (or application) clearly explained? Partly

Is the description of the method technically sound? Partly
Are sufficient details provided to allow replication of the method development and its use by others? Yes If any results are presented, are all the source data underlying the results available to ensure full reproducibility? Yes Are the conclusions about the method and its performance adequately supported by the findings presented in the article? Partly No competing interests were disclosed.

Competing Interests:
Referee Expertise: Computational biology, systems biology, networks, statistics, infectious disease, immunology 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.
Author Response 20 May 2018 , University of Cape Town, South Africa Armin Deffur Many thanks for a thorough review. We made the following changes to address the stated concerns.
In response to points 1 and 2: We added adjusted P-values used as cutoffs in the analytic pipeline in Table 1, and have added a  new Table 3, providing detail on the bipartite graphs, based on the recommendation.
3. Figure 1 has been edited to clarify that both pathways and Chaussabel modules are gene sets forming prior information.
4. The term "validation study" has been replaced with "example study" 5. We emphasise in the Results section that the "virtual cells" and the cell matrices are new information we were able to generate based on the shared relationships found in ANIMA