Evaluation criteria

For the following analyses, we assessed the quality of the input networks and network integration method outputs using three evaluation criteria: (1) gene coannotation prediction; (2) gene module detection and (3) supervised gene function prediction. First, we used an established precision-recall evaluation strategy^22,23 to determine how well gene– gene relationships produced by the given method overlapped with gene pairs coannotated to the same term in a particular functional standard. Second, we evaluated the capacity of each method to produce biological modules by comparing clusters computed from the output of each method to known modules such as protein complexes, pathways and biological processes. These two evaluations measure the intrinsic quality of the outputs generated by the integration methods, that is without training any additional models on top of the outputs. Finally, the supervised gene function prediction evaluation determines how discriminative the method outputs are for predicting known gene functions. Here, a portion of the genes and corresponding labels (known functional classes such as protein complex membership) were held out and used to evaluate the accuracy of a support vector machine (SVM) classifier²⁴, which is trained on the remaining gene features, output from the given integration method, to predict the held-out labels⁷. This constitutes an extrinsic evaluation, indicating how effectively the method outputs can be used in conjunction with an additional classification model.

In the following experiments, to ensure a fair choice of hyperparameters across BIONIC and the integration methods we compared to, we performed a hyperparameter optimization step using an independent set of Schizosaccharomyces pombe networks as inputs^25–27 and a set of Gene Ontology (GO) curated pombe protein complexes²⁸ for evaluation. The best performing hyperparameters for each approach were used (Methods).

Evaluation of BIONIC features and input networks

We first used the unsupervised BIONIC to integrate three diverse yeast networks: a comprehensive network of correlated genetic interaction profiles (4,529 genes, 33,056 interactions)²², a coexpression network derived from transcript profiles of yeast strains carrying deletions of transcription factors (1,101 genes, 14,826 interactions)²¹ and a protein–protein interaction network obtained from an affinity-purification mass-spectrometry assay (2,674 genes, 7,075 interactions)²⁰, which combine for a total of 5,232 unique genes and 53,351 unique interactions (Fig. 2 and Supplementary Data File 2). Compared to the input networks, BIONIC integrated features have equivalent or superior performance on all evaluation criteria over three different functional benchmarks: IntAct protein complexes²⁹, Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways³⁰ and GO Biological Processes²⁸ (Fig. 2a and Supplementary Data File 3). The coannotation and module detection benchmarks contain between 1,786 and 4,170 genes overlapping the integration results. The module detection benchmarks define between 107 and 1,809 modules. The IntAct, KEGG and GO Biological Process gene function prediction benchmarks cover 567, 1,770 and 1,211 genes overlapping the integration results, and 48, 53 and 63 functional classes, respectively (Supplementary Data File 2). As an additional test, BIONIC produces high-quality features that accurately predict a diverse set of yeast biological process annotations per gene²² (Fig. 2b). Some categories in this last test do better than others. These performance patterns were mirrored in the individual input networks, indicating that this is the result of data quality, rather than method bias.

Open in a separate window

Fig. 2 |

Comparison of BIONIC integration to three input networks.

a, Functional evaluations for three yeast networks, and unsupervised BIONIC integration. Data are presented as mean values. Error bars indicate the 95% confidence interval for $n=10$ independent samples. Number of captured modules are indicated above the module detection bars as determined by a 0.5 overlap (Jaccard) score cutoff. b, Evaluations over high-level functional categories, split by category. Numbers above columns indicate gene overlap with integration results and the average performance of each method is reported (right of each row). c, Top row: comparison of overlap scores between known complexes and predicted modules. Each point is a protein complex. The axes indicate the overlap (Jaccard) score, where 0 indicates no members of the complex were captured, and 1.0 indicates the complex was captured perfectly. The diagonal indicates equivalent performance. Points above the diagonal are complexes where BIONIC outperforms the given network, and points below the diagonal are complexes where BIONIC underperforms. The arrows indicate the LSM2–7 complex, shown in d. A Venn diagram describes the overlap of captured complexes (score of 0.5 or higher) between the input networks and BIONIC integration. Numbers in brackets denote the total captured complexes for each method. Bottom row: the distribution of overlap scores between predicted and known complexes. The dashed line indicates the mean. d, Functional relationships between predicted LSM2–7 complex members and genes in the local neighborhood, as given by the three input networks and BIONIC integration. The predicted cluster best matching the LSM2–7 complex in each network is circled. The overlap score of the predicted module with the LSM2–7 complex is shown. Edges correspond to protein–protein interactions in PPI²⁰, Pearson correlation between gene profiles in coexpression²¹ and genetic interaction²² networks, and cosine similarity between gene features in the BIONIC integration. The complete LSM2–7 complex is depicted on the right. Edge weight corresponds to the strength of the functional relationship. PPI, protein–protein interaction; COEX, coexpression; GI, genetic interaction; BP, biological process.

We observed that features obtained through BIONIC network integration often outperformed the individual input networks at capturing functional modules (Fig. 2a) and captured more modules (Fig. 2c and Supplementary Data File 4), demonstrating the use of the combined features over individual networks for downstream applications such as module detection. Here we treated the network adjacency profiles (rows in the adjacency matrix) as gene features. We then examined how effectively the input networks and integrated BIONIC features captured known protein complexes, by matching each individual known complex to its best-matching predicted module and quantifying the overlap (Fig. 2c). We then compared the overlap scores from each network to the BIONIC overlap scores to identify complexes where BIONIC performs either better or worse than the input networks. Of 344 protein complexes tested, BIONIC strictly improved 196, 309 and 222 complex predictions and strictly worsened 82, 17 and 98 complex predictions compared to the input protein–protein interaction, coexpression and genetic interaction networks, respectively. The distributions of complex overlap scores for each dataset indicate that BIONIC predicts protein complexes more accurately than the input networks on average. Indeed, if we use an overlap score of 0.5 or greater to indicate a successfully captured complex, the integrated BIONIC features, containing information from three networks, capture 121 complexes, compared to 88, 3 and 74 complexes for the individual protein–protein interaction, coexpression and genetic interaction networks, respectively (Fig. 2c). We also repeated this module analysis while optimizing the clustering parameters on a per-module basis, an approach that tests how well each network and BIONIC perform at capturing modules under optimal clustering conditions for each module. Here too, the integrated BIONIC features capture more modules and with a greater average overlap score than the individual input networks (Extended Data Figs. 4 and ​and55 and Supplementary Data File 5).

To better understand how BIONIC is able to improve functional gene module detection compared to the input networks, we examined the LSM2–7 complex, which was identified in our module detection evaluation (Fig. 2a) as an example to show how BIONIC effectively combines gene–gene relationships across different networks and recapitulates known biology. The LSM2–7 complex localizes to the yeast nucleoli and is involved in the biogenesis or function of the small nucleolar RNA SNR5 (ref. ³¹). LSM2–7 is made up of the protein products of six genes: LSM2, LSM3, LSM4, LSM5, LSM6 and LSM7. We found that the cluster that best matched the LSM2–7 complex in each input network only captures a subset of the full complex (Supplementary Data File 4). The BIONIC module, however, contains five out of six members of the LSM2–7 complex, along with two additional members: LSM1 and PAT1, which are functionally associated with the LSM2–7 complex³². The missing member, LSM5, is in the local neighborhood of the cluster in the BIONIC feature space. We examined the best-matching clusters and their local neighborhood, consisting of genes that show a direct interaction with predicted members of the LSM2–7 complex, in the input networks, and in a profile similarity network obtained from the integrated BIONIC features of these networks (Fig. 2d). We found that both the PPI and genetic interaction networks captured two members of the LSM2–7 complex, with two additional members in the local neighborhood. The coexpression network only identified one complex member, and the local neighborhood of the best-matching module did not contain any additional known complex members. Finally, BIONIC used the interaction information across input networks to better identify the LSM2–7 module, with the addition of two functionally related proteins. This analysis demonstrates the use of BIONIC for identifying meaningful biological modules by effectively combining information across input networks. Indeed, when we optimized the module detection procedure to specifically resolve the LSM2–7 complex, we found that BIONIC was able to capture the complex with a higher overlap score (0.83) than any of the input networks (0.33, 0.17 and 0.50 for the PPI, coexpression and genetic interactions networks, respectively), and it outperformed other integration methods (0.43, 0.22, 0.44, 0.60 and 0.68 for the Union, iCell, deepNF, Mashup and multi-node2vec methods, respectively) (Supplementary Data File 5).

We also performed an analysis to examine how BIONIC features encode the input networks and found that the input networks are generally encoded uniformly across feature dimensions (Supplementary Note 1 and Extended Data Fig. 6).

Evaluation of BIONIC and established unsupervised integration methods

We compared network integration results from the unsupervised BIONIC (Fig. 2) to those derived from several different established integration approaches: a naive union of networks (Union), a nonnegative matrix tri-factorization approach (iCell)², a deep learning multi-modal autoencoder (deepNF)¹¹, a low-dimensional diffusion state approximation approach (Mashup)⁷ and a multi-network extension of the node2vec (ref. ¹³) model (multi-node2vec)³³ (Fig. 3). These unsupervised integration methods cover a wide range of methodologies and the main possible output types (networks for Union and iCell, features for deepNF, Mashup and multi-node2vec). BIONIC performs as well as, or better than the tested integration methods across all evaluation types and benchmarks (Fig. 3a). We also evaluated BIONIC and the other integration approaches on a per-biological process basis (Fig. 3b). Here we found BIONIC generally outperforms the established integration approaches on each biological process, with the exception of several biological processes when compared to deepNF. Averaging over the performance for each biological process, we found BIONIC performs on par with deepNF (average precision of 0.53 for BIONIC compared to 0.52 for deepNF). DeepNF performs competitively on the per-biological process evaluations (Fig. 3b), but it underperforms on the global performance evaluations (Fig. 3a). The per-biological process evaluations assess how well a method predicts large-scale biological process coannotation, whereas the global performance evaluations measure how well a method predicts smaller-scale functional modules (that is, protein complexes). This discrepancy in performance indicates deepNF is able to capture broad-scale functional organization, but it fails to resolve smaller functional modules. BIONIC performs well on both of these evaluations, however, indicating it can learn gene features that resolve both broad and detailed functional organization. To ensure these results are consistent under different input networks, we integrated a set of yeast-two-hybrid networks and found similar performance patterns (Supplementary Note 2 and Extended Data Fig. 7).

Open in a separate window

Fig. 3 |

Comparison of BIONIC to existing integration approaches.

a, Coannotation prediction, module detection and gene function prediction evaluations for three yeast networks integrated by the tested unsupervised network integration methods. The input networks and evaluation standards are the same as in Fig. 2. Data are presented as mean values. Error bars indicate the 95% confidence interval for $n=10$ independent samples. Numbers above the module detection bars indicate the number of captured modules, as determined by a 0.5 overlap (Jaccard) score cutoff. b, Evaluation of integrated features using high-level functional categories, split by category. Numbers above columns indicate gene overlap with integration results and the average performance of each method across categories is reported (right of each row). PPI, protein– protein interaction; BP, biological process.

Evaluation of BIONIC in a semisupervised setting

We also tested how BIONIC performs in a semisupervised setting (Fig. 4). Here, we compared BIONIC trained with no labeled data (unsupervised), BIONIC trained with a held-out set of functional labels given by IntAct, KEGG and GO (semisupervised), and a supervised integration algorithm using the same labels (GeneMANIA⁵). For each of these methods, we integrated the yeast protein–protein interaction, coexpression and genetic interaction networks from the Fig. 2 analysis. In each benchmark, 20% of genes (IntAct, KEGG, GO) were randomly held out and used as a test set, while the remaining 80% of genes were used for training. The unsupervised BIONIC did not use any gene label information for training, but it was evaluated using the same test set as the supervised methods to ensure a consistent performance comparison. To control for variability in the train-test set partitioning, this procedure was repeated ten times and the average performance across test sets was reported (Methods). We found that adding labeled data can substantially improve the features BIONIC learns and these features also outperform the integration results produced by the supervised GeneMANIA method. We also found that even without labeled data, BIONIC performs as well as, or exceeds, GeneMANIA performance. Notably, the performance of the unsupervised and semisupervised BIONIC is similar for gene function prediction. This indicates unsupervised BIONIC features are already sufficiently discriminative for classifiers to perform well. Thus, BIONIC can be used effectively in both an unsupervised and semisupervised setting, which demonstrates its versatility as a biological network integration algorithm. We also analyzed the performance of BIONIC in a scenario where labels are noisy (Supplementary Note 3 and Extended Data Fig. 8).

Open in a separate window

Fig. 4 |

Supervised performance of BIONIC compared with an existing supervised integration approach.

Performance comparison between a supervised network integration algorithm trained with labeled data (GeneMANIA), BIONIC trained without any labeled data (unsupervised) and BIONIC trained with labeled data (semisupervised). Bars indicate the average performance over ten trials of random train-test splits for the given benchmark (Methods). Data are presented as mean values. Error bars indicate the 95% confidence interval. $n=10$ independent samples for the coannotation prediction and gene function prediction evaluations, and $n=100$ for the module detection evaluation. BP, biological process.

Scalability of BIONIC and established integration approaches

An effective integration algorithm should be able to scale to many network inputs and networks with many nodes. To test network input scalability, we randomly sampled progressively larger sets of yeast gene coexpression networks (Fig. 5a and Supplementary Data File 2) and assessed the performance of the resulting integrations of these sets. We similarly tested node scalability by randomly subsampling progressively larger gene sets of four human protein–protein interaction networks^34–37 (Fig. 5b and Supplementary Data File 2). BIONIC can integrate numerous networks (Fig. 5a), and networks with many nodes (Fig. 5b), outperforming all other methods assessed for progressively more and larger networks. To achieve this scalability, BIONIC takes advantage of the versatile nature of deep learning technology by learning features for small batches of genes and networks at a time, reducing the computational resources required for any specific training step. To learn gene features over large networks, BIONIC learns features for random subsets of genes at each training step and randomly subsamples the local neighborhoods of these genes to perform the graph convolution (Methods), maintaining a small overall computational footprint. This subsampling allows BIONIC to integrate networks with many genes, whereas methods such as Mashup can only do so with an approximate algorithm that reduces integration performance (Supplementary Fig. 1). To integrate many networks, BIONIC uses a network-wise sampling approach, where a random subset of networks is integrated at a time during each training step. This reduces the number of parameter updates required at once, since only GCNs corresponding to the subsampled networks are updated in a given training step. We also tested the extent of BIONIC scalability in terms of computational resources (Supplementary Note 4 and Extended Data Fig. 9).

Open in a separate window

Fig. 5 |

Network quantity and network size performance comparison across integration methods.

a, Performance comparison of unsupervised integration methods across different numbers of randomly sampled yeast coexpression input networks on KEGG pathways gene coannotations. b, Performance comparison of unsupervised integration methods across four human protein–protein interaction networks for a range of subsampled nodes (genes) on CORUM complexes protein coannotations. In these experiments, the Mashup method failed to scale to seven or more networks (a) and 4,000 or more nodes (b), as indicated by the absence of bars in those cases (Methods). Data are presented as mean values. Error bars indicate the 95% confidence interval for $n=10$ independent samples. multi-n2v, multi-node2vec.

Implementation details

BIONIC was implemented using PyTorch⁶⁰, a popular Python-based deep learning framework, and relies on functions and classes from the PyTorch Geometric library⁶¹. It uses the Adam⁶² optimizer to train and update its weights. To be scalable in the number of networks, BIONIC uses an optional network batching approach where subsets of networks are sampled and integrated at each training step. The sampling procedure is designed so that each network is integrated exactly once per training step. Network batching yields a constant memory footprint at the expense of increased runtime with no empirical degradation of feature quality. This feature is provided for additional scalability over what is demonstrated in this work and was not used in any of our reported experiments. Additionally, BIONIC is scalable in the number of network nodes. It uses a node sampling approach (equivalent to mini-batch training, where nodes are samples) to learn features for subsets of nodes in a network, and a neighborhood sampling procedure to subsample node neighborhoods. Node sampling ensures only part of a network needs to be retained in memory at a time while neighborhood sampling reduces the effective higher-order neighborhood size in sequential GAT passes, again reducing the number of nodes required to be retained in memory at any given time, further reducing BIONIC’s memory footprint.

For very large networks where the initial node feature matrix (that is, the identity matrix) cannot fit into memory due to limitations with PyTorch matrix operations, BIONIC incorporates a singular value decomposition-based approximation. First, the union of networks is computed by creating a network that contains the nodes and edges of all input networks. If an edge occurs in multiple networks, the maximum weight is used. A low-dimensional singular value decomposition approximation of the normalized Laplacian matrix of the union network is computed and used as the initial node features for each network. Finally, BIONIC uses sparse representations of network adjacency matrices (except for the input node feature matrix, see above), further reducing memory footprint. All BIONIC integration experiments in this paper were run on an NVIDIA Titan Xp graphical processing unit with 12 GB of video RAM, no more than 16 GB of system RAM and a single 2.4 GHz Intel Xeon CPU.

Network preprocessing

The yeast protein–protein interaction network²⁰ and human protein–protein interaction networks^34–37 were obtained from BioGRID⁶³, genetic interaction profiles²² were obtained directly from the published supplementary data of Costanzo et al.²², and gene expression profiles²¹ were obtained from the SPELL database⁶⁴. These networks were chosen since they had the most functional information compared to other networks in their class (that is, protein–protein interaction networks, coexpression networks and genetic interaction networks). To create a network from the genetic interaction profiles, genes with multiple alleles were collapsed into a single profile by taking the maximum profile values across allele profiles. Pairwise Pearson correlation between the profiles was then calculated, and gene pairs with a correlation magnitude greater than or equal to 0.2 were retained as edges, as established²². For the gene expression profiles, networks were constructed by retaining gene pairs with a profile Pearson correlation magnitude in the 99.5th percentile. Coexpression and genetic interaction networks had their edge weights normalized to the range [0, 1].

Obtaining integrated results

The naive union of networks benchmark was created by taking the union of node sets and edge sets across input networks. For edges common to more than one network, the maximum weight was used. For all other methods, automated hyperparameter optimization was performed to ensure hyperparameters were chosen consistently and fairly. Here, a S. pombe genetic interaction network²⁷, coexpression network²⁶ and protein–protein interaction network²⁵ were used as inputs to the integration methods. To perform one iteration of the hyperparameter optimization, a random combination of hyperparameters was uniformly sampled over a range of reasonable values for each method and used to integrate the three pombe networks. The integration results were then evaluated using a pombe protein complex standard (obtained from https://www.pombase.org/data/annotations/Gene_ontology/GO_complexes/Complex_annotation.tsv). The evaluations consisted of a coannotation prediction, module detection and gene function prediction assessment (Evaluation methods). This procedure was repeated for 50 combinations of hyperparameters, for each method. For methods that produced features (deepNF, Mashup, multi-node2vec and BIONIC), a feature dimension of 512 was used to ensure results were comparable across methods. For methods that required a batch size parameter (deepNF and BIONIC), the batch size was set to 2,048 to ensure reasonable computation times. Hyperparameter combinations were then ranked for each method across the three evaluation types and the hyperparameter combination corresponding to the highest average rank across evaluation types was chosen. The hyperparameter optimization results are found in Supplementary Data File 1. Note that the union method was not included in the hyperparameter optimization because it has no hyperparameters. Additionally, the Mashup method used 44 hyperparameter combinations rather than 50, as six hyperparameter combinations exhausted the available memory resources and did not complete.

All integration results reported were obtained by integrating networks using the set of hyperparameters identified in the hyperparameter optimization procedure. BIONIC features used in the Figs. 2, ​,33 and ​and66 analyses are found in Supplementary Data File 8. Coannotation prediction, module detection and gene function prediction standards used in Figs. 2–5 are found in Supplementary Data File 9.

Benchmark construction

Functional benchmarks were derived from GO Biological Process ontology annotations, KEGG pathways and IntAct complexes for yeast, and CORUM complexes for human (Supplementary Data File 3). Analyses were performed using positive and negative gene pairs, clusters or functional labels obtained from the standards as follows: the GO Biological Process benchmark was produced by filtering Inferred from Electronic Annotation (IEA) annotations, as they are known to be lower quality, removing genes with dubious open reading frames and filtering terms with more than 30 annotations (to prevent large terms, such as those related to ribosome biogenesis, from dominating the analysis⁶⁵). We found the performance evaluations to be robust to this threshold (Supplementary Fig. 6). For the coannotation benchmark, all gene pairs sharing at least one annotation were retained as positive pairs, while all gene pairs not sharing an annotation were considered to be negative pairs. KEGG, IntAct and CORUM benchmarks were produced analogously, without size filtering.

For the module detection benchmark, clusters were defined as the set of genes annotated to a particular term, for each standard. Modules of size 1 (singletons) were removed from the resulting module sets as they are uninformative. For the per-module analyses in Fig. 2c, Extended Data Figs. 4 and ​and55 and Supplementary Data Files 4 and 5, we also removed any modules of size 2 since these modules had highly variable Jaccard scores.

The supervised standards were obtained by treating each gene annotation as a class label, leading to genes with multiple functional classes (that is, a multi-label classification problem). The standards were filtered to only include classes with 20 or more members for GO Biological Process and KEGG, or ten members for IntAct. This was done to remove classes with very few data points, ensuring more robust evaluations.

The granular function standard in Figs. 2b, ​,3b3b and ​and66 was obtained from the Costanzo et al.²² supplementary materials. Any functional category with fewer than 20 gene members was removed from the analysis to ensure only categories with robust evaluations were reported.

Evaluation methods

We used a precision-recall-based coannotation framework to evaluate individual networks and integrated results. We used precision-recall instead of receiving operator curve because of the substantial imbalance of positives and negatives in the pairwise benchmarks for which the receiving operator curve would overestimate performance. Here, we computed the pairwise cosine similarities between gene profiles in each network or integration result. Due to the high-dimensionality of the datasets, cosine similarity is a more appropriate measure than Euclidean distance since the contrast between data points is reduced in high-dimensional spaces under Euclidean distance⁶⁶. Precision-recall operator points were computed by varying a similarity threshold, above which gene or protein pairs are considered positives and below which pairs are considered negative. Each set of positive and negative pairs was compared to the given benchmark to compute precision and recall values. To summarize the precision-recall curve into a single metric, we computed average precision (AP) given by:

where $n$ is the number of operator points (that is, similarity thresholds) and $P_i$ and $R_i$ are the precision and recall values at operator point $i$, respectively. This gives the average of precision values weighted by their corresponding improvements in recall. We chose this measure over the closely related area under the precision-recall curve measure since this interpolates between operator points and tends to overestimate actual performance⁶⁷.

The module detection evaluation was performed by clustering the integrated results from each method and comparing the coherency of resulting clusters with the module-based benchmarks. Since the benchmarks contain overlapping modules (that is, one gene can be present in more than one module) that prevent the use of many common clustering evaluation metrics (since these metrics assume unique assignment of gene to cluster), the module sets are subsampled during the evaluation to ensure there are no overlapping modules (the original module sets are used as-is for the per-module-optimized experiments in Extended Data Fig. 5 and Supplementary Data File 5). Next, the integrated results are hierarchically clustered with a range of distance metrics (Euclidean and cosine), linkage methods (single, average and complete) and thresholds to optimize benchmark comparisons over these clustering parameters (this is done for all methods that are compared). The resulting benchmark-optimized cluster sets are compared to the benchmark module sets by computing adjusted mutual information: an information theoretic comparison measure that is adjusted to normalize against the expected score from random clustering. The highest adjusted mutual information score for each integration approach is reported, ensuring the optimal cluster set for each dataset across clustering parameters is used for the comparison and that our results are not dependent on clustering parameters. Finally, this procedure is repeated ten times to control for differences in scores due to the cluster sampling procedure. The sets of clustering parameter-optimized BIONIC clusters obtained from the Fig. 2 integration for each standard are in Supplementary Data File 4.

To perform the supervised gene function prediction evaluation, ten trials of five-fold cross validation were performed using SVM classifiers each using a radial basis function kernel²⁴. The classifiers were trained on a set of gene features obtained from the given integration method with corresponding labels given by the IntAct, KEGG and GO Biological Process supervised benchmarks in a one-versus-all fashion (since each individual gene has multiple labels). Each classifier’s regularization and gamma parameters were tuned in the validation step. For each trial, the classifier results were evaluated on a randomized held-out set consisting of 10% of the gene features not seen during training or validation and the resulting classification accuracy was reported. We repeated this entire procedure for a random forest⁶⁸ and a gradient boosted trees⁶⁹ classifier and found BIONIC also outperforms the compared integration methods, indicating the SVM classifier is not biased toward improving BIONIC performance (Supplementary Fig. 7).

The granular functional evaluations in Figs. 2b and ​and3b3b were generated by computing the average precision (as mentioned in the precision-recall evaluation framework description) for the gene subsets annotated to the given functional categories.

To perform the module comparison analysis in Fig. 2c, we additionally applied the module detection analysis performed in Fig. 2a to the input networks. Here, the interaction profiles of the networks were treated as gene features and the clustering parameters were optimized to best match the IntAct complexes standard. We compared the resulting module sets from the input networks and BIONIC features to known protein complexes given by the IntAct standard. For each complex in the standard, we reported the best-matching predicted module in each dataset as determined by the overlap (Jaccard) score between the module and the known complex (Supplementary Data File 4). To generate the Venn diagram, we defined a complex to have been captured in the dataset if it had an overlap score of 0.5 or greater with a predicted module.

To perform the LSM2–7 module analysis in Fig. 2d, we analyzed the predicted module in each dataset that had the highest overlap score with the LSM2–7 complex. We created a network from the BIONIC features by computing the cosine similarity between all pairs of genes and setting all similarities below 0.5 to zero. The resulting nonzero values were then treated as weighted edges to form a network. We extracted a subnetwork from each of the protein–protein interaction, coexpression, genetic interaction and newly created BIONIC networks, consisting of the best scoring predicted module and the genes showing direct interactions with those in the predicted module. We laid out these networks using the spring-embedded layout algorithm in Cytoscape⁷⁰. The edges in the protein–protein interaction network correspond to direct, physical interactions, and the edges in the coexpression and genetic interaction networks correspond to the pairwise Pearson correlation of the gene profiles, as described above.

To perform the semisupervised network integration experiment in Fig. 4, we first generated randomized train and test sets. Here, 20% of genes were randomly held out in each gene function benchmark (IntAct, KEGG and GO Biological Process) separately and retained for downstream evaluations. These benchmarks consist of functional labels for a set of yeast genes (protein complex membership in IntAct, pathway membership in KEGG and biological process annotation in GO Biological Process) and are the same benchmarks used in the gene function prediction evaluation (Figs. 2a and ​and3a).3a). The remaining 80% of genes were used for training GeneMANIA and BIONIC. To generate test sets for the coannotation prediction benchmarks, we removed any coannotations where both genes were present in the training set. To generate test sets for the module detection benchmarks, we removed any modules consisting entirely of genes in the training set. We then integrated the three yeast networks from the Figs. 2 and ​and33 analysis (a protein–protein interaction²⁰, gene coexpression²¹ and genetic interaction network²²) using the supervised GeneMANIA, BIONIC without using any labeled data (unsupervised) and a semisupervised mode of BIONIC that uses the labeled data (semisupervised). Each integration result was then evaluated using the held-out test data. For the coannotation prediction and module detection evaluations, the integrated features from BIONIC (both unsupervised and semisupervised) and the integrated network from GeneMANIA were evaluated. Both GeneMANIA and the semisupervised BIONIC generate gene label predictions directly, without the need for an additional classifier such as in the Fig. 2a gene function prediction evaluation. However, the unsupervised BIONIC does not generate gene label predictions (since it is given no labeled information to begin with). To ensure a consistent comparison with GeneMANIA and the semisupervised BIONIC, we trained a classification head on top of the unsupervised BIONIC. The classification head architecture is identical to the semisupervised BIONIC classification head, however, in the unsupervised case we only allow gradients from the classification loss objective to backpropagate to the classification head, not the rest of the model. This ensures a comparable classification model can be trained on top of the unsupervised BIONIC model, without the labeled data affecting the model weights such as in the semisupervised case. GeneMANIA does not generate multi-label predictions, and so we used GeneMANIA to generate label predictions for each class individually and then performed Platt scaling to convert these binary class predictions to multi-label predictions^7,71. The gene function prediction evaluations were then performed by comparing the gene label predictions from the integration methods, to the held-out test labels. This entire procedure, starting with the train-test set partitioning, to the final evaluations, was repeated a total of ten times to control for performance variability due to the partitioning procedure.

Network scaling experiment

To perform the network scaling experiment, we uniformly sampled subsets of the yeast coexpression networks (Supplementary Data File 2). We performed ten integration trials for each network quantity, and these trials were paired (that is, each method integrated the same randomly sampled sets of networks). The average precision scores of the resulting integrations with respect to the KEGG pathways coannotation standard (Supplementary Data Files 3) were then reported. The Mashup method did not scale to the seven-network input size or beyond on a machine with 64 GB of RAM.

Node scaling experiment

The node scaling experiment was performed by uniformly subsampling the nodes of four large human protein–protein interaction networks^34–37 (Supplementary Data File 2) for a range of node quantities and integrating these subsampled networks. Ten trials of subsampling were performed for each number of nodes (paired, as above) and the average precision scores with respect to the CORUM complexes coannotation standard (Supplementary Data File 3) were reported. The Mashup method did not scale to 4,000 nodes or beyond on a machine with 64 GB of RAM.

Gene chemical sensitivity predictions

Chemical–genetic profiles against a diagnostic set of 310 nonessential yeast gene deletion mutants were obtained from a previous study¹⁸. The genes were chosen using the COMPRESS-GI algorithm, which selected a set of 157 genes capturing most of the functional information within genome-wide genetic interaction data⁷², along with 153 genes that were manually selected to complement the set. Haploid deletion mutants for the gene set were constructed in a genetic background that conferred drug hypersensitivity ($pdr\mathit{1}\Delta pdr\mathit{3}\Delta snq\mathit{2}\Delta$) using synthetic genetic analysis technology, and each mutant strain was barcoded with a unique 20 bp DNA identifier adjacent to a common priming site. The mutant collection was grown and stored as a pooled library in yeast extract peptone-glycerol (15% v/v). A set of approximately 10,000 compounds from the RIKEN Natural Product Depository (NPDepo) were interrogated. Screens were done in 96-well format, where a single well contained the entire pool of 310 mutants at a density of 4.65 × 10⁵ cells per ml and 196 μl of YPGal media (1% yeast extract, 2% peptone, 2% galactose). Each well was treated with 2 μl of compound (1 mg ml⁻¹ stock dissolved in dimethylsulfoxide (DMSO)). After 48 h of growth in 30 °C, genomic DNA was extracted from each compound-treated pool with an automated high-throughput nucleic acid purification robot (QIAcube HT, Qiagen). Mutant-specific barcodes and well-specific index tags were PCR-amplified using multiplex primers and a communal U2 primer. PCR products were pooled in 768-plex and gel-purified from 2% agarose gels using a Geneclean III kit. Amplicons were quantified using a Kapa quantitative PCR kit and were sequenced with an Illumina Hiseq 2500 machine at the RIKEN Center for Life Science Technologies. Sequencing data was processed using the BEAN-counter software⁴¹, which generated chemical–genetic interaction $z$-scores normalized against DMSO-only (1% DMSO) treated samples. False discovery rates were estimated for biological process²² predictions, for each compound, and those compounds with a false discovery rate of <25% were retained, resulting in a set of 1,522 compounds and 289 genes (high-confidence set)¹⁸. Next, interquartile range (IQR) scores were calculated from the chemical–genetic scores as follows:

Here, ${\text{CGs}}_i$ is the chemical–geneticscore for the $i$th replicate, $\overline{\text{CGs}}$ is the median of all chemical–genetic scores, $Q{3}_{\text{CGs}}$ is the 75th percentile of chemical–genetic scores and $Q{1}_{\text{CGs}}$ is the 25th percentile of chemical–genetic scores. Tukey’s test⁷³ was used to determine outliers based on the IQR of the distribution of IQR scores in the screen. Genes with at least one replicate that had a negative (sensitive) chemical–genetic score more than three times the IQR of the compound profile (that is, ‘outlier’ genes) were retained.

To predict chemical–genetic interactions using BIONIC, we first selected a set of 50 compounds to generate predictions on and experimentally validate. For each diagnostic pool compound, we filtered out any genes not present in the integrated BIONIC features (the same features used for the Figs. 2 and ​and33 analyses, referred to as PEG features). Any compounds with fewer than two outlier sensitive genes were then removed. For each of the remaining compounds, we randomly split the sensitive genes into train and test sets. Next, for a given compound, we computed BIONIC predictions for the test set genes. We did this by averaging the corresponding BIONIC PEG features for each gene in the training set under a cosine distance metric to get a representative feature vector in gene feature space for the given compound. The BIONIC predictions for the compound were then obtained by identifying the top 20 nearest genes to this feature vector (excluding genes in the training set).

To obtain a score for the BIONIC predictions, an ordered Fisher’s exact test was performed between the test set genes and the BIONIC predictions as follows:

where

$p$ corresponds to the minimum $P$ value obtained for progressively larger subsets of BIONIC’s 20 predictions, starting from the top prediction to the full set of 20 predictions. $n$ is the number of total predictions made by BIONIC and $k$ is the number of those predictions that are correct. $k_i$ corresponds to the number of correct predictions for the first $i$ genes in the BIONIC predictions. $f$ is the probability mass function of the hypergeometric distribution. Here, $K$ corresponds to the number of genes found to be sensitive to the given compound. $N$ is the total number of yeast essential genes in the analysis, specifically, essential genes for which TS mutants could be made and are also present in the BIONIC features (847 total genes). We chose the ordered Fisher’s exact test over the commonly used unordered version because BIONIC produces a ranked list of predictions. Taking into account the ordering of BIONIC predictions is a fairer assessment, since, for example, a compound may only have a small number of sensitive genes (fewer than 20). In this case, BIONIC’s top predictions may include these essential genes, however, an unordered Fisher’s exact test would not consider this ranking and treat the full set of 20 predictions as equivalent, whereas the ordered test would consider the ranking.

The above process was repeated five times for new randomly sampled train and test gene splits, or up to the maximum number of train-test splits possible for compounds with fewer than five sensitive genes. Final $P$ values were obtained for each compound by averaging over the $P$ values from each trial. Compounds were ranked by most significant $P$ values and the top 50 compounds were selected for further screening. Sensitive essential gene predictions for a given compound were then generated by using the full set of sensitive diagnostic pool genes as the training set, computing a representative compound feature vector by averaging the training set BIONIC gene features, and identifying the top 20 nearest essential genes to this compound feature vector.

The BIONIC gene chemical sensitivity predictions were benchmarked against experimental data obtained from chemical–genetic screens using a collection of TS mutants for essential genes. We previously constructed a drug-hypersensitive, barcoded set of TS mutants for 1,181 TS alleles spanning 837 essential genes⁴⁰. Similar to the diagnostic set of nonessential genes, this collection also contained the $\text{pdr}1\text{Δpdr}3\text{Δsnq}2\text{Δ}$ triple deletion and a 20 bp barcode was inserted next to a common priming site upstream of a natMX cassette integrated at the $\text{pdr}3\text{Δ}$ locus. We conducted chemical–genetic screens against the 50 compounds initially selected for BIONIC analysis using the same method that was used to generate the diagnostic set profiles, except that the TS mutant pools were incubated at 25 °C instead of 30 °C for 48 h. We calculated chemical–genetic interaction $z$-scores and removed nonspecific technical effects using BEAN-counter software⁴¹. IQR scores were calculated as described above. Negative (sensitive) interactions that were more than four times the IQR (classified as ‘far outliers’) were used to validate the gene chemical sensitivities predicted by BIONIC.

The significance of BIONIC sensitive essential gene predictions for each compound was determined by using an ordered Fisher’s exact test, as detailed above. The Benjamini–Hochberg procedure⁷⁴ was applied to the resulting $P$ values at a false discovery rate of 5%.

To generate biological process²² predictions as reported in Fig. 6a,​,b,b, a Fisher’s exact test was performed between the full set of 20 BIONIC gene predictions and biological process gene annotations. We used the same annotations as in Figs. 2b and ​and3b3b (ref. ²²). If the BIONIC sensitive gene predictions were enriched for one or more bioprocesses, and these bioprocesses overlapped with the annotated bioprocess of the true sensitive genes, we considered this a correct bioprocess prediction. To generate the random benchmark in Fig. 6b, the gene labels of the BIONIC integrated features were randomly permuted and new essential sensitive gene predictions for the 50 selected compounds were generated in the same manner as the original BIONIC predictions (detailed above). This process was repeated for 1,000 random gene label permutations to generate the benchmark distributions. The circle plot in Fig. 6d was produced by first hierarchically clustering the integrated BIONIC gene features, subsetted to essential genes annotated to the glycosylation, protein folding/targeting, cell wall biosynthesis bioprocess. Two clustering thresholds were chosen to generate clusters, broadly indicating the hierarchical organization of the BIONIC gene features. The first, most granular clustering threshold was adaptively chosen to generate clusters best-matching known protein complexes, as defined by the IntAct complexes standard²⁹. For each protein complex in the standard, the clustering threshold was optimized to produce the cluster best matching this protein complex. For clusters not matching known complexes, the largest complex optimized threshold was used. The second, higher clustering threshold was set to a cophenetic distance of 0.9.

The BIONIC essential gene sensitivity predictions can be found in Supplementary Data File 7.

Quantification of $\text{β}$-1,6-glucan levels

Wild-type ($his\mathit{3}\Delta$ in the BY4741 background) and the $kre\mathit{6}\Delta$ strain (YOC5627) of $S$. cerevisiae were grown in yeast extract peptone at 25 °C with shaking at 200 r.p.m. to 1 × 10⁷ cells per ml. Wild-type cells were treated with 5 μg ml⁻¹ jervine (J0009; Tokyo Chemical Industry, Tokyo, Japan) for 4 h. We used jervine since it is chemically similar to NP329 and is more commercially available than NP329. The samples were centrifuged at 15,000g for 3 min, and the supernatant was discarded. The pellet was washed and suspended in PBS, adjusted to 1 × 10⁶ cells per ml and autoclaved for 20 min. After centrifugation at 15,000g for 1 min, the supernatant was stored on ice (sample A) and the pellet was further extracted. The $\text{β}$-1,6-glucan was extracted from the pellet using a slightly modified version of the protocol of Kitamura et al.⁷⁵. First, 500 ml of 10% TCA was added to the culture, which was incubated on ice for 10 min. After centrifugation at 15,000g for 3 min, the samples were washed twice with deionized water. The pellet was suspended in 500 μl of 1 N NaOH and incubated at 75 °C for 1 h. The solution was mixed with 500 μl of 1 M HCl and Tris buffer (10 mM Tris-HCl, pH 7). After centrifugation at 15,000g for 1 min, the supernatant was stored on ice (sample B). The total amounts of $\text{β}$-1,6-glucan in samples A and B were measured according to the method of Yamanaka et al.⁷⁶.

We thank B. Andrews, M. Costanzo and C. Myers for their insightful comments. We also thank M. Fey for adding important features to the PyTorch Geometric library for us. This work was supported by NRNB (US National Institutes of Health, National Center for Research Resources grant number P41 GM103504). Funding for continued development and maintenance of Cytoscape is provided by the US National Human Genome Research Institute under award number HG009979. This work was also supported by the Canadian Institutes of Health Research Foundation grant number FDN-143264, US National Institutes of Health grant number R01HG005853 and joint funding by Genome Canada (OGI-163) and the Ministry of Economic Development, Job Creation and Trade, under the program Bioinformatics and Computational Biology. This work was supported by the National Research Council of Canada through the AI for Design program. This work was supported by CIFAR AI Chair programs. This work was also supported by JSPS KAKENHI grant numbers JP15H04483 (C.B. and Y.O.), JP17H06411 (C.B. and Y.Y.), JP18K14351 (K.I.-N.), JP19H03205 (Y.O.), JP20K07487 (D.Y.) and a RIKEN Foreign Postdoctoral Fellowship (S.C.L.). This research was enabled in part by support provided by SciNet and the Digital Research Alliance of Canada.