Biogeography & environmental conditions shape bacteriophage-bacteria networks across the human microbiome

Viruses and bacteria are critical components of the human microbiome and play important roles in health and disease. Most previous work has relied on studying bacteria and viruses independently, thereby reducing them to two separate communities. Such approaches are unable to capture how these microbial communities interact, such as through processes that maintain community robustness or allow phage-host populations to co-evolve. We implemented a network-based analytical approach to describe phage-bacteria network diversity throughout the human body. We built these community networks using a machine learning algorithm to predict which phages could infect which bacteria in a given microbiome. Our algorithm was applied to paired viral and bacterial metagenomic sequence sets from three previously published human cohorts. We organized the predicted interactions into networks that allowed us to evaluate phage-bacteria connectedness across the human body. We observed evidence that gut and skin network structures were person-specific and not conserved among cohabitating family members. High-fat diets appeared to be associated with less connected networks. Network structure differed between skin sites, with those exposed to the external environment being less connected and likely more susceptible to network degradation by microbial extinction events. This study quantified and contrasted the diversity of virome-microbiome networks across the human body and illustrated how environmental factors may influence phage-bacteria interactive dynamics. This work provides a baseline for future studies to better understand system perturbations, such as disease states, through ecological networks. Author Summary The human microbiome, the collection of microbial communities that colonize the human body, is a crucial component to health and disease. Two major components to the human microbiome are the bacterial and viral communities. These communities have primarily been studied separately using metrics of community composition and diversity. These approaches have failed to capture the complex dynamics of interacting bacteria and phage communities, which frequently share genetic information and work together to maintain ecosystem homestatsis (e.g. kill-the-winner dynamics). Removal of bacteria or phage can disrupt or even collapse those ecosystems. Relationship-based network approaches allow us to capture this interaction information. Using this network-based approach with three independent human cohorts, we were able to present an initial understanding of how phage-bacteria networks differ throughout the human body, so as to provide a baseline for future studies of how and why microbiome networks differ in disease states.


Introduction
Viruses and bacteria are critical components of the human microbiome and play important roles in health and disease.Bacterial communities have been associated with disease states, including a range of skin conditions (1), acute and chronic wound healing conditions (2,3), and gastrointestinal diseases, such as inflammatory bowel disease (4,5), Clostridium difficile infections (6) and colorectal cancer (7,8).Altered human viromes (virus communities consisting primarily of bacteriophages) also have been associated with diseases and perturbations, including inflammatory bowel disease (5,9), periodontal disease (10), spread of antibiotic resistance (11), and others (12)(13)(14)(15)(16)(17).Viruses act in concert with their microbial hosts as a single ecological community (18).Viruses influence their living microbial host communities through processes including lysis, host gene expression modulation (19), influence on evolutionary processes such as horizontal gene transfer (22) or antagonistic co-evolution (26), and alteration of ecosystem processes and elemental stoichiometry (27).
Previous human microbiome work has focused on bacterial and viral communities, but have reduced them to two separate communities by studying them independently (5,9,10,(12)(13)(14)(15)(16)(17).This approach fails to capture the complex dynamics of interacting bacteria and phage communities, which frequently share genetic information and work together to maintain stable ecosystems.Removal of bacteria or phages can disrupt or even collapse those ecosystems (18,(28)(29)(30)(31)(32)(33)(34)(35)(36)(37).Relationship-based network approaches allow us to capture this interaction information.Studying such bacteria-phage interactions through community-wide networks built from inferred relationships could offer further insights into the drivers of human microbiome diversity across body sites and enable the study of human microbiome network dynamics overall.
In this study, we characterized human-associated bacterial and phage communities by their inferred relationships using three published paired virus and bacteria-dominated whole community metagenomic datasets (13,14,38,39).We leveraged machine learning and graph theory techniques to establish and explore the human bacteria-phage network diversity therein.This approach built upon previous large-scale phage-bacteria network analyses by inferring interactions from metagenomic datasets, rather than culture-dependent data (33), which is limited in the scale of possible experiments and analyses.
Our metagenomic interaction inference model improved upon previous models of phage-host predictions that have utilized a variety of techniques, such as linear models to predict bacteria-phage co-occurrence using taxonomic assignments (40), and nucleotide similarity models that were applied to both whole virus genomes (41) and related clusters of whole and partial virus genomes (42).Our approach uniquely included protein interaction data and was validated based on experimentally determined positive and negative interactions (i.e. who does and does not infect whom).Through this approach we were able to provide a basic understanding of the network dynamics associated with phage and bacterial communities on and in the human body.By building and utilizing a microbiome network, we found that different people, body sites, and anatomical locations not only support distinct microbiome membership and diversity (13,14,38,39,(43)(44)(45), but also support ecological communities with distinct communication structures and propensities toward community instability.Through an improved understanding of network structures across the human body, we empower future studies to investigate how these communities dynamics are influenced by disease states and the overall impact they may have on human health.

Cohort Curation and Sample Processing
We studied the differences in virus-bacteria interaction networks across healthy human bodies by leveraging previously published shotgun sequence datasets of purified viral metagenomes (viromes) paired with bacteria-dominated whole community metagenomes.Our study contained three datasets that explored the impact of diet on the healthy human gut virome (14), the impact of anatomical location on the healthy human skin virome (13), and the viromes of monozygotic twins and their mothers (38,39).We selected these datasets because their virome samples were subjected to virus-like particle (VLP) purification, which removed contaminating DNA from human cells, bacteria, etc.To this end, the publishing authors employed combinations of filtration, chloroform/DNase treatment, and cesium chloride gradients to eliminate organismal DNA (e.g.bacteria, human, fungi, etc) and thereby allow for direct assessment of both the extracellular and fully-assembled intracellular virome (Supplemental Figure S1 A-B) (14,39).Each research group reported quality control measures to ensure the purity of the virome sequence datasets, using both computational and molecular techniques (e.g.16S rRNA gene qPCR) (Table S1).These reports confirmed that the virome libraries consisted of highly purified virus genomic DNA.
The bacterial and viral sequences from these studies were quality filtered and assembled into contigs.We further grouped the related bacterial and phage contigs into operationally defined units based on their k-mer frequencies and co-abundance patterns, similar to previous reports (Supplemental Figure S2 -S3) (42).We referred to these operationally defined groups of related contigs as operational genomic units (OGUs).Each OGU represented a genomically similar sub-population of either bacteria or phages.Contig lengths within clusters ranged between 10 3 and 10 5.5 bp (Supplemental Figure S2 -S3).
While supplementing the previous quality control measures (Table S1) we found that, in light of the rigorous purification and quality control during sample collection and preparation, 77% (228 / 298 operational genomic units) still had some nucleotide similarity to a given bacterial reference genome (e-value < 10 −25 ).
As absence of bacterial contamination had been confirmed by sensitive molecular methods (Table S1), we interpreted this as evidence that the majority of the gut and skin bacteriophages were temperate and thereby shared elements with bacterial reference genomes that contained similar integrated phages when sequenced-a trend previously reported (14).Additionally, we identified two OGUs as being complete phages using the stringent Virsorter phage identification algorithm (class 1 confidence group) (47).

Evaluating the Model to Predict Phage-Bacteria Interactions
We predicted which phage OGUs infected which bacterial OGUs using a random forest model trained on experimentally validated infectious relationships from six previous publications (41,(48)(49)(50)(51)(52).Only bacteria and phages were used in the model.The training set contained 43 diverse bacterial species and 30 diverse phage strains, including both broad and specific ranges of infection (Supplemental Figure S4 A -B).While it is true that there are more known phages that infect bacteria, we were limited by the information confirming which phages do not infect certain bacteria and attempted to keep the numbers of positive and negative interactions similar.Phages with linear and circular genomes, as well as ssDNA and dsDNA genomes, were included in the analysis.Because we used DNA sequencing studies, RNA phages were not considered (Supplemental Figure S4 C-D).This training set included both positive relationships (a phage infects a bacterium) and negative relationships (a phage does not infect a bacterium).This allowed us to validate the false positive and false negative rates associated with our candidate models, thereby building upon previous work that only considered positive relationships (41).
Four phage and bacterial genomic features were used in a random forest model to predict infectious relationships between bacteria and phages: 1) genome nucleotide similarities, 2) gene amino acid sequence similarities, 3) bacterial Clustered Regularly Interspaced Short Palindromic Repeat (CRISPR) spacer sequences that target phages, and 4) similarity of protein families associated with experimentally identified protein-protein interactions (53).The resulting random forest model was assessed using nested cross validation, and the median area under its receiver operating characteristic (ROC) curve was 0.788, the median model sensitivity was 0.952, and median specificity was 0.615 (Figure 1 A).The most important predictor in the model was amino acid similarity between genes, followed by nucleotide similarity of whole genomes (Figure 1 B).Protein family interactions were moderately important to the model, and CRISPRs were largely uninformative, due to the minimal amount of identifiable CRISPRs in the dataset and their redundancy with the nucleotide similarity methods (Figure 1 B).Approximately one third of the training set relationships yielded no score and therefore were unable to be assigned an interaction prediction (Figure 1 C).
central each sample's network was on average.We accomplished this by utilizing two common centrality metrics: degree centrality and closeness centrality.Degree centrality, the simplest centrality metric, was defined as the number of connections each phage made with each bacterium.We supplemented measurements of degree centrality with measurements of closeness centrality.Closeness centrality is a metric of how close each phage or bacterium is to all of the other phages and bacteria in the network.A higher closeness centrality suggests that the effects of genetic information or altered abundance would be more impactful to all other microbes in the system.A network with higher average closeness centrality also indicates an overall greater degree of connections, which suggests a greater resilience against instability.This is because more highly connected networks are more stable as a result of pathway dependencies and the unlikelihood that a randomly removed bacteria or phage would cause major divisions across the network (30,58).We used this information to calculate the average connectedness per sample, which was corrected for the maximum potential degree of connectedness.Unfortunately our dataset was insufficiently powered to make strong conclusions toward this hypothesis, but this is an interesting observation that warrants further investigation.
Using our limited sample set, we observed that the gut microbiome network structures associated with high-fat diets appeared less connected than those of low-fat diets, although a greater sample size will be required to more properly evaluate this trend (Figure 2 A-B).Five subjects were available for use, all of which had matching bacteria and virome datasets and samples from 8-10 days following the initiation of their diets.
High-fat diets appeared to exhibit reduced degree centrality (Figure 2 A), suggesting bacteria in high-fat environments were targeted by fewer phages and that phage tropism was more restricted.High-fat diets also appeared to exhibit decreased closeness centrality (Figure 2 B), indicating that bacteria and phages were more distant from other bacteria and phages in the community.This would make genetic transfer and altered abundance of a given phage or bacterium less capable of impacting other bacteria and phages within the network.
In addition to diet, we found preliminary evidence that obesity influenced network structure.This was done using the three mother samples available from the twin sample set, all of which had matching bacteria and phage samples and confirmed BMI information.The obesity-associated network appeared to have a higher degree centrality (Figure 2 C), but less closeness centrality than the healthy-associated networks (Figure 2 D).These results suggested that the obesity-associated networks are less connected, having microbes further from all other microbes within the community.This again comes with the caveat that this is a preliminary observation with too small of a sample size to make more substantial claims.

Individuality of Microbial Networks
Skin and gut community membership and diversity are highly personal, with people remaining more similar to themselves than to other people over time (13,59,60).We therefore hypothesized that this personal conservation extended to microbiome network structure.We addressed this hypothesis by calculating the degree of dissimilarity between each subject's network, based on phage and bacteria abundance and centrality.We quantified phage and bacteria centrality within each sample graph using the weighted eigenvector centrality metric.This metric defines central phages as those that are highly abundant (A O as defined in the methods) and infect many distinct bacteria which themselves are abundant and infected by many other phages.Similarly, bacterial centrality was defined as those bacteria that were both abundant and connected to numerous phages that were themselves connected to many bacteria.We then calculated the similarity of community networks using the weighted eigenvector centrality of all nodes between all samples.Samples with similar network structures were interpreted as having similar capacities for maintaining stability and transmitting genetic material.
We used this network dissimilarity metric to test whether microbiome network structures were more similar within people than between people over time.We found that gut microbiome network structures clustered by person (ANOSIM p-value = 0.005, R = 0.958, Figure 3 A).Network dissimilarity within each person over the 8-10 day sampling period was less than the average dissimilarity between that person and others, although this difference was not statistically significant (p-value = 0.125, Figure 3 B).Four of the five available subjects were used because one of the subjects was not sampled at the initial time point.The lack of statistical confidence was likely due to the small sample size of this dataset.
Although there was evidence for gut network conservation among individuals, we found no evidence for conservation of gut network structures within families.The gut network structures were not more similar within families (twins and their mothers; intrafamily) compared to other families (other twins and mothers; inter-family) (p-value = 0.312, Figure 3 C).In addition to the gut, skin microbiome network structure was strongly conserved within individuals (p-value < 0.001, Figure 3 D).This distribution was similar when separated by anatomical sites.Most sites were statistically significantly more conserved within individuals (Supplemental Figure S6).

Human Skin Landscape
Extensive work has illustrated differences in diversity and composition of the healthy human skin microbiome between anatomical sites, including bacteria, virus, and fungal communities (13,44,59).These communities vary by degree of skin moisture, oil, and environmental exposure.As viruses are known to influence microbial diversity and community composition, we hypothesized that microbe-virus network structure would be specific to anatomical sites, as well.To test this, we evaluated the changes in network structure between anatomical sites within the skin dataset.
The average centrality of each sample was quantified using the weighted eigenvector centrality metric.
Intermittently moist skin sites (dynamic sites that fluctuate between being moist and dry) were significantly less connected than the more stable moist and sebaceous environments (p-value < 0.001, Figure 4 A).Also, skin sites that were occluded from the environment were much more highly connected than those that were constantly exposed to the environment or only intermittently occluded (p-value < 0.001, Figure 4 B).
To supplement this analysis, we compared the network signatures using the centrality dissimilarity approach described above.The dissimilarity between samples was a function of shared relationships, degree of centrality, and bacteria/phage abundance.When using this supplementary approach, we found that network structures significantly clustered by moisture, sebaceous, and intermittently moist status (Figure 4 C,E).
Occluded sites were significantly different from exposed and intermittently occluded sites, but there was no difference between exposed and intermittently occluded sites (Figure 4 D,F).These findings provide further support that skin microbiome network structure differs significantly between skin sites.

Discussion
Foundational work has provided a baseline understanding of the human microbiome by characterizing bacterial and viral diversity across the human body (13,14,(43)(44)(45)61).Here, we offer an initial understanding of how phage-bacteria networks differ throughout the human body, so as to provide a baseline for future studies of how and why microbiome networks differ in disease states.We developed and implemented a network-based analytical model to evaluate the basic properties of the human microbiome through bacteria and phage relationships, instead of membership or diversity alone.This enabled the application of network theory to provide a new perspective on complex ecological communities.We utilized metrics of connectivity to model the extent to which communities of bacteria and phages interact through mechanisms such as horizontal gene transfer, modulated bacterial gene expression, and alterations in abundance.
Just as gut microbiome and virome composition and diversity are conserved in individuals (13,43,44,60), gut and skin microbiome network structures were conserved within individuals over time.Gut network structure was not conserved among family members.These findings suggested that the community properties inferred from microbiome interaction network structures, such as stability (meaning a more highly connected network is more stable because a randomly removed bacteria or phage node is less likely to divide or disintegrate (30,58) the overall network), the potential for horizontal gene transfer between members, and co-evolution of populations, were person-specific.These properties may be impacted by personal factors ranging from the body's immune system to external environmental conditions, such as climate and diet.
We observed evidence supporting the ability of environmental conditions to shape gut and skin microbiome interaction network structure by observing that diet and skin location were associated with altered network structures.We found evidence that diet was sufficient to alter gut microbiome network connectivity, although this needs to be interpreted as a case observation, due to the small sample size.Although our sample size was small, our findings provided some preliminary evidence that high-fat diets were less connected than low-fat diets and that high-fat diets therefore may lead to less stable communities with a decreased ability for microbes to directly influence one another.We supported this finding with the observation that obesity may have been associated with decreased network connectivity.Together these findings suggest the food we eat may not only impact which microbes colonize our guts, but may also impact their interactions with infecting phages.Further work will be required to characterize these relationships with a larger cohort.
In addition to diet, the skin environment also influenced the microbiome interaction network structure.
Network structure differed between environmentally exposed and occluded skin sites.The sites under greater environmental fluctuation and exposure (the exposed and intermittently exposed sites) were less connected and therefore were predicted to have a higher propensity for instability.Likewise, intermittently moist sites demonstrated less connectedness than the more stable moist and sebaceous sites.Together these data suggested that body sites under greater degrees of fluctuation harbored less connected, potentially less stable microbiomes.This points to a link between microbiome and environmental stability and warrants further investigation.
While these findings take us an important step closer to understanding the microbiome through interspecies relationships, there are caveats to and considerations regarding the approach.First, as with most classification models, the infection classification model developed and applied is only as good as its training set -in this case, the collection of experimentally-verified positive and negative infection data, where genomes of all members are fully sequenced.Large-scale experimental screens for phage and bacteria infectious interactions that report high-confidence negative interactions (i.e., no infection) are desperately needed, as they would provide more robust model training and improved model performance.Furthermore, just as we have improved on previous modeling efforts, we expect that new and creative scoring metrics will be integrated into this model to improve future performance.
Second, although our analyses utilized the best datasets currently available for our study, this work was done retrospectively and relied on existing data up to seven years old.These archived datasets were limited by the technology and costs of the time.For example, the diet and twin studies, relied on multiple displacement amplification (MDA) in their library preparations-an approach used to overcome the large nucleic acids requirements typical of older sequencing library generation protocols.It is now known that MDA results in biases in microbial community composition (62), as well as toward ssDNA viral genomes (63,64), thus rendering the resulting microbial and viral metagenomes largely non-quantitative.Future work that employs larger sequence datasets and that avoids the use of bias-inducing amplification steps will build on and validate our findings, as well as inform the design and interpretation of further studies.
Finally, the networks in this study were built using operational genomic units (OGUs), which represented groups of highly similar bacteria or phage genomes or clustered genome fragments.Similar clustering definition and validation methods, both computational and experimental, have been implemented in other metagenomic sequencing studies, as well (42,(65)(66)(67).These approaches could offer yet another level of sophistication to our network-based analyses.While this operationally defined clustering approach allows us to study whole community networks, our ability to make conclusions about interactions among specific phage or bacterial species or populations is inherently limited, compared to more focused, culture-based studies such as the work by Malki et al (68).Future work must address this limitation, e.g., through improved binning methods and deeper metagenomic shotgun sequencing, but most importantly through an improved conceptual framing of what defines ecologically and evolutionarily cohesive units for both phage and bacteria (69).Defining operational genomic units and their taxonomic underpinnings (e.g., whether OGU clusters represent genera or species) is an active area of work critical to the utility of this approach.As a first step, phylogenomic analyses have been performed to cluster cyanophage isolate genomes into informative groups using shared gene content, average nucleotide identity of shared genes, and pairwise differences between genomes (70).Such population-genetic assessment of phage evolution, coupled with the ecological implications of genome heterogeneity, will inform how to define nodes in future iterations of the ecological network developed here.Even though we are hesitant to speculate on phage host ranges at low taxonomic levels in our dataset, the data does aggree with previous reports of instances of broad phage host range (68,71).
Together our work takes an initial step towards defining bacteria-virus interaction profiles as a characteristic of human-associated microbial communities.This approach revealed the impacts that different human environments (e.g., the skin and gut) can have on microbiome connectivity.By focusing on relationships between bacterial and viral communities, they are studied as the interacting cohorts they are, rather than as independent entities.While our developed bacteria-phage interaction framework is a novel conceptual advance, the microbiome also consists of archaea and small eukaryotes, including fungi and Demodex mites (1, 72)-all of which can interact with human immune cells and other non-microbial community members (73).
Future work will build from our approach and include these additional community members and their diverse interactions and relationships (e.g., beyond phage-bacteria).This will result in a more robust network and a more holistic understanding of the evolutionary and ecological processes that drive the assembly and function of the human-associated microbiome.

Materials & Methods
trimmed using the Fastx toolkit (v0.0.14) to exclude bases with quality scores below 33 and shorter than 75 bp (74).Paired end reads were filtered to exclude sequences missing their corresponding pair using the get_trimmed_pairs.pyscript available in the source code.

Contig Assembly
Contigs were assembled using the Megahit assembly program (v1.0.6) (75).A minimum contig length of 1 kb was used.Iterative k-mer stepping began at a minimum length of 21 and progressed by 20 until 101.All other default parameters were used.

Contig Abundance Calculations
Contigs were concatenated into two master files prior to alignment, one for bacterial contigs and one for phage contigs.Sample sequences were aligned to phage or bacterial contigs using the Bowtie2 global aligner (v2.2.1) (76).We defined a mismatch threshold of 1 bp and seed length of 25 bp.Sequence abundance was calculated from the Bowtie2 output using the calculate_abundance_from_sam.plscript available in the source code.

Operational Genomic Unit Binning
Contigs often represent large fragments of genomes.In order to reduce redundancy and the resulting artificially inflated genomic richness within our dataset, it was important to bin contigs into operational units based on their similarity.This approach is conceptually similar to the clustering of related 16S rRNA sequences into operational taxonomic units (OTUs), although here we are clustering contigs into operational genomic units (OGUs) (61).
Contigs were clustered using the CONCOCT algorithm (v0.4.0) (77).Because of our large dataset and limits in computational efficiency, we randomly subsampled the dataset to include 25% of all samples, and used these to inform contig abundance within the CONCOCT algorithm.CONCOCT was used with a maximum of 500 clusters, a k-mer length of four, a length threshold of 1 kb, 25 iterations, and exclusion of the total coverage variable.OGU abundance (A O ) was obtained as the sum of the abundance of each contig (A j ) associated with that OGU.The abundance values were length corrected such that: Where L is the length of each contig j within the OGU.

Operational Genomic Unit Identification
To confirm a lack of phage sequences in the bacterial OGU dataset, we performed blast nucleotide alignment of the bacterial OGU representative sequences using an e-value < 10 −25 , which was stricter than the 10 −10 threshold used in the random forest model below, against all of the phage reference genomes available in the EMBL database.We used a stricter threshold because we know there are genomic similarities between bacteria and phage OGUs from the interactive model, but we were interested in contigs with high enough similarity to references that they may indeed be from phages.We also performed the converse analysis of aligning phage OGU representative sequences to EMBL bacterial reference genomes.Finally, we ran both the phage and bacteria OGU representative sequences through the Virsorter program (1.0.3) to identify phages (all default parameters were used), using only those in the high confidence identification category "class 1" (47).

Open Reading Frame Prediction
Open reading frames (ORFs) were identified using the Prodigal program (V2.6.2) with the meta mode parameter and default settings (78).
As a result of gene transfer or phage genome integration during infection, phages may share genes with their bacterial hosts, providing us with evidence of phage-host pairing.We identified shared genes between bacterial and phage genomes by assessing amino acid similarity between the genes using the Diamond protein alignment algorithm (v0.7.11.60) (82).The mean alignment bitscores for each genome pair were recorded for use in our classification model.

Protein -Protein Interactions
The final method used for predicting infectious interactions between bacteria and phages was the detection of pairs of genes whose proteins are known to interact.We assigned bacterial and phage genes to protein families by aligning them to the Pfam database using the Diamond protein alignment algorithm.We then identified which pairs of proteins were predicted to interact using the Pfam interaction information within the Intact database (53).The mean bitscores of the matches between each pair were recorded for use in the classification model.

Interaction Network Construction
The bacteria and phage operational genomic units (OGUs) were scored using the same approach as outlined above.The infectious pairings between bacteria and phage OGUs were classified using the random forest model described above.The predicted infectious pairings and all associated metadata were used to populate a graph database using Neo4j graph database software (v2.3.1)(83).This network was used for downstream community analysis.

Centrality Analysis
We quantified the centrality of graph vertices using three different metrics, each of which provided different information graph structure.When calculating these values, let G(V, E) be an undirected, unweighted graph with |V | = n nodes and |E| = m edges.Also, let A be its corresponding adjacency matrix with entries a ij = 1 if nodes V i and V j are connected via an edge, and a ij = 0 otherwise.
Briefly, the closeness centrality of node V i is calculated taking the inverse of the average length of the shortest paths (d) between nodes V i and all the other nodes V j .Mathematically, the closeness centrality of node V i is given as: The distance between nodes (d) was calculated as the shortest number of edges required to be traversed to move from one node to another.
Intuitively, the degree centrality of node V i is defined as the number of edges that are incident to that node: where a ij is the ij th entry in the adjacency matrix A.
The eigenvector centrality of node V i is defined as the i th value in the first eigenvector of the associated adjacency matrix A. Conceptually, this function results in a centrality value that reflects the connections of the vertex, as well as the centrality of its neighboring vertices.
The centralization metric was used to assess the average centrality of each sample graph G. Centralization was calculated by taking the sum of each vertex V i 's centrality from the graph maximum centrality C w , such that: The values were corrected for uneven graph sizes by dividing the centralization score by the maximum theoretical centralization (T) for a graph with the same number of vertices.

Network Relationship Dissimilarity
We assessed similarity between graphs by evaluating the shared centrality of their vertices, as has been done previously.More specifically, we calculated the dissimilarity between graphs G i and G j using the Bray-Curtis dissimilarity metric and eigenvector centrality values such that: Where C ij is the sum of the lesser centrality values for those vertices shared between graphs, and C i and C j are the total number of vertices found in each graph.This allows us to calculate the dissimilarity between graphs based on the shared centrality values between the two graphs.

Statistics and Comparisons
Differences in intrapersonal and interpersonal network structure diversity, based on multivariate data, were calculated using an analysis of similarity (ANOSIM).Statistical significance of univariate Eigenvector centrality differences were calculated using a paired Wilcoxon test.
Statistical significance of differences in univariate eigenvector centrality measurements of skin virome-microbiome networks were calculated using a pairwise Wilcoxon test, corrected for multiple hypothesis tests using the Holm correction method.Multivariate eigenvector centrality was measured as the mean differences between cluster centroids, with statistical significance measured using an ANOVA and post hoc Tukey test.Eigen Centrality Differences in Mean Levels of Group

Figure 1 :
Figure 1: Summary of Multi-Study Network Model.(A) Median ROC curve (dark red) used to create the microbiome-virome infection prediction model, based on nested cross validation over 25 random iterations.The maximum and minimum performance are shown in light red.(B) Importance scores associated with the metrics used in the random forest model to predict relationships between bacteria and phages.The importance score is defined as the mean decrease in accuracy of the model when a feature (e.g.Pfam) isexcluded.Features include the local gene alignments between bacteria and phage genes (denoted blastx; the blastx algorithm in Diamond aligner), local genome nucleotide alignments between bacteria and phage OGUs, presence of experimentally validated protein family domains (Pfams) between phage and bacteria OGUs, and CRISPR targeting of bacteria toward phages (CRISPR).(C) Proportions of samples included (gray) and excluded (red) in the model.Samples were excluded from the model because they did not yield any scores.Those interactions without scores were automatically classified as not having interactions.(D) Bipartite visualization of the resulting phage-bacteria network.Phage OGUs are presented in orange,bacteria OGUs in red, and their interaction edges are represented as connecting lines.This network includes information from all three published studies.(E) Network diameter (measure of graph size; the greatest number of traversed vertices required between two vertices), (F) number of vertices, and (G) number of edges (relationships) for the total network (orange) and the individual study sub-networks (diet study = red, skin study = yellow, twin study = green).

Figure 2 :
Figure 2: Impact of Diet and Obesity on Gut Network Structure.(A) Quantification of average degree centrality (number of edges per node) and (B) closeness centrality (average distance from each node to every other node) of gut microbiome networks of subjects limited to exclusively high-fat or low-fat diets.Each point represents the centrality from a human subject stool sample that was collected 8-10 days following the beginning of their defined diet.There are five samples here, compared to the four in figure 3, because one of the was only sampled post-diet, providing us data for this analysis but not allowing us to compare to a baseline for figure 3. (C) Quantification of average degree centrality and (D) closeness centrality between obese and healthy adult women from the Twin gut study.Each point represents a stool sample taken from one of the three adult woman confirmed as obese or healthy and with matching virus and bacteria data.

Figure 3 :
Figure 3: Intrapersonal vs Interpersonal Network Dissimilarity Across Different Human Systems.(A) NMDS ordination illustrating network dissimilarity between subjects over time.Each sample is colored by subject, with each colored sample pair collected 8-10 days apart.Dissimilarity was calculated using the Bray-Curtis metric based on abundance weighted eigenvector centrality signatures, with a greater distance representing greater dissimilarity in bacteria and phage centrality and abundance.Only four subjects were included here, compared to the five used in figure 2, because one of the subjects was missing the inital sampling time point and therefore lacked temporal sampling.(B)Quantification of gut network dissimilarity within the same subject over time (intrapersonal) and the mean dissimilarity between the subject of interest and all other subjects (interpersonal).The p-value is provided near the bottom of the figure.(C) Quantification of gut network dissimilarity within subjects from the same family (intrafamily) and the mean dissimilarity between subjects within a family and those of other families (interfamily).Each point represents the inter-family and intra-family dissimilarity of a twin or mother that was sampled over time.(D) Quantification of skin network dissimilarity within the same subject and anatomical location over time (intrapersonal) and the mean dissimilarity between the subject of interest and all other subjects at the same time and the same anatomical location (interpersonal).All p-values were calculated using a paired Wilcoxon test.

Figure 4 :
Figure4: Impact of Skin Micro-Environment on Microbiome Network Structure.(A) Notched box-plot depicting differences in average eigenvector centrality between moist, intermittently moist, and sebaceous skin sites and (B) occluded, intermittently occluded, and exposed sites.Notched box-plots were created using ggplot2 and show the median (center line), the inter-quartile range (IQR; upper and lower boxes), the highest and lowest value within 1.5 * IQR (whiskers), outliers (dots), and the notch which provides an approximate 95% confidence interval as defined by 1.58 * IQR / sqrt(n).Sample sizes for each group were: Moist = 81, Sebaceous = 56, IntMoist = 56, Occluded = 106, Exposed = 61, IntOccluded = 26.(C) NMDS ordination depicting the differences in skin microbiome network structure between skin moisture levels and (D) occlusion.Samples are colored by their environment and their dissimilarity to other samples was calculated as described in figure3.(E) The statistical differences of networks between moisture and (F) occlusion status were quantified with an anova and post hoc Tukey test.Cluster centroids are represented by dots and the extended lines represent the associated 95% confidence intervals.Significant comparisons (p-value < 0.05) are colored in red, and non-significant comparisons are gray.

Figure S1 :
Figure S1: Sequencing Depth Summary.Number of sequences that aligned to (A) Phage and (B) Bacteria operational genomic units per sample and colored by study.

Figure S2 :
Figure S2: Contig Summary Statistics.Scatter plot heat map with each hexagon representing the abundance of contigs.Contigs are organized by length on the x-axis and the number of aligned sequences on the y-axis.

Figure S3 :
FigureS3: Operational Genomic Unit Summary Statistics.Scatter plot with operational genomic unit clusters organized by average contig length within the cluster on the x-axis and the number of contigs in the cluster on the y-axis.Operational genomic units of (A) bacteriophages and (B) bacteria are shown.

Figure S4 :
Figure S4: Summary information of validation dataset used in the interaction predictive model.A) Categorical heat-map highlighting the experimentally validated positive and negative interactions.Only bacteria species are shown, which represent multiple reference strains.Phages are labeled on the x-axis and bacteria are labeled on the y-axis.B) Quantification of bacterial host strains known to exist for each phage.C) Genome strandedness and D) linearity of the phage reference genomes used for the dataset.

Figure S5 :
FigureS5: Structure of the interactive network.Metadata relationships to samples (Phage Sample ID and Bacteria Sample ID) included the associated time point, the study, the subject the sample was taken from, and the associated disease.Infectious interactions were recorded between phage and bacteria operational genomic units (OGUs).Sequence count abundance for each OGU within each sample was also recorded.

Figure S6 :
Figure S6: Intrapersonal vs Interpersonal Dissimilarity of the Skin.Quantification of skin network dissimilarity within the same subject and anatomical location over time (intrapersonal) and the mean dissimilarity between the subject of interest and all other subjects at the same time and the same anatomical location (interpersonal), separated by each anatomical site (forehead [Fh], palm [Pa], toe web [Tw], umbilicus [Um], antecubital fossa [Ac], axilla [Ax], and retroauricular crease [Ra]).P-value was calculated using a paired Wilcoxon test.