Blast (Altschul et al., 1990) was used to compare the sequence similarity of all genome genes in one genus. In order to quickly and strictly find the similar genes, we directly adopted the reciprocal best hits comparison and controlled the identities and cover ratios to 65%. Taking Flavobacterium as an example, a total of 16,830 orthologs were obtained using OrthoFinder 2.0 (Emms and Kelly, 2019), and 51,247 nodes were obtained using our method, of which 50,998 nodes contained only one orthologs (99.5%). Therefore, our method is stricter than directly using orthologs. Node2vec algorithm (Grover and Leskovec, 2016) was used to generate graph embeddings for each node.

Positional relationship pairs between two genes from each genome were constructed using the data of coding DNA sequence (CDS) (Figure 2). Through the correspondence between genes and nodes, all positional pairs were connected into a single gene synteny network, in which there could be more than one connection between two nodes. The adjacency matrix was constructed in accordance with the number of connections between nodes.

Given an undirected graph with node feature matrix X and adjacency matrix A, the graph convolution operation (Kipf et al., 2016) is defined as: H =   σ ( D 1 2 A ^ D − 1 2 X W ) (3) A ^ = A + I ,       D i i =   ∑ j A i j ^ (4) where I is the identity matrix, W is the matrix of trainable weights in the neural network, X is the feature matrix before the update, H is the feature matrix after the update, and σ is the activation function (ReLU). The graph convolution operation iteratively calculates the weighted average of the node attributes of the neighbors of the current node to obtain the new feature matrix of the node. In this framework, the features of unlabeled nodes (nodes without assigned functional pathways) and the features of nearby labeled nodes (nodes with assigned functional pathways) are mixed to be propagated through the synteny network diagram. If two nodes have the same neighbor structure and neighbor features, their embedded feature matrix H will be exactly the same.

Python’s PyTorch Geometric Module was used to implement PPA-GCN. Multiple graph convolutional layers can be stacked to enable learning on a larger domain structure. After testing, a three-layer stack was found to perform the best. The two-class cross entropy was used as the loss function because of the multilabel nature of the problem.

The original input was fed into the framework, and 50 epochs of random sampling verification training were performed with the test set. The nodes with an average cross-validation accuracy rate of less than 30% are removed from the training set, and nodes and labels with a assignment stability of 90% in the test set (that is, the same label is assigned more than 45 times) are added to the training set. After many iterations, when the number of nodes in the training set reached more than 90% of the total number of nodes in the gene synteny network, the training was considered to have reached a stable state, and the final assignment results were output.

The degree is defined as the number of all edge connections of a node in a graph, describing the first-order connection degree of the node. The degree distribution is an overall description of the nodes in a network, that is, the probability distribution or statistical distribution of the node degrees.

The clustering coefficient is used to describe the degree of clumping among the vertices of a network. Specifically, it is the degree of interconnection among the adjacent nodes of a node, describing the second-order connection degree of the node. For node i with degree k _i , the local clustering coefficient is defined as: C i =   2 L i k i ( k i − 1 )   (5) where L _i is the number of connections among the k _i neighbors of node i. The overall aggregation coefficient of the network is characterized as the average value of the aggregation coefficients of all nodes.

PRA is the accuracy at the node level and is defined as the proportion of genes with at least one label assigned correctly.

The TLPR is the accuracy at the label level and is defined as the number of correctly assigned labels divided by the total number of labels.

When a label is predicted for a node, we assign weights in accordance with the assignment probability, sum the WPRA of each label of a node to obtain the WPRA of that node, and divide by the total number of nodes to obtain the overall WPRA: w p r e d i c t i o n =   1 N ∑ k ∈ T i 2 ( I + 1 − k ) I   ( I + 1 ) (6) where N is the total number of nodes, I is the number of labels for node i, T _i is the order of the correct label probabilities assigned for node i (from large to small), and k is the k-th ranked probability label that was assigned correctly.

The kappa coefficient is often used for testing consistency, that is, whether the assignment effect of the model is consistent with the actual classification effect. Its value is between -1 and 1. When the value is greater than 0.6, it is considered substantial, and when it is greater than 0.8, it is considered almost perfect. The calculation of the kappa coefficient is based on the confusion matrix: k a p p a =   p 0 − p e 1 − p e (7) p 0 =   ∑ i M i i ∑ i j M i j ,     p e =   ∑ i M i . M . i ( ∑ i j M i j ) 2 (8) where M is the confusion matrix of the assignment results.

The Hamming distance is measure of the distance between the assigned and real labels, with a value between 0 and 1. A distance of 0 means that the assigned results are exactly the same as the real results, and a distance of 1 means that the model’s results are completely opposite to the desired results. This indicator is calculated as the number of erroneously assigned labels divided by the total number of labels.

This coefficient is an indicator for comparing the similarity of two finite sets, defined as the size of the intersection of two label sets (the true label set and the assigned label set) divided by the size of the union. When this coefficient is 1, the assigned results are completely consistent with the actual situation; when the coefficient is 0, the assigned results are completely inconsistent with the actual situation.

We tested PPA-GCN with 5-fold cross-validation on three data sets. PPA-GCN achieved prediction rates of assignment (PRAs) of 84.8% (Flavobacterium), 77.0% (Pseudomonas) and 71.0% (Staphylococcus) on the three prokaryotic bacterial genera (Figure 3). According to the evaluation index results (Table 1), PPA-GCN is well adapted to all three genera.

In addition, we compared PPA-GCN with five other machine learning methods. deepNF (Gligorijevic et al., 2018), Mashup (Cho et al., 2016) and Pseudo2GO (Fan and Zhang, 2020) are three deep learning methods that use graph information for function prediction. Support vector machines (SVM) and deep neural networks (DNN) are two machine learning models that are not based on graph information. Using the Staphylococcus genome as the test data set, all methods use the same features in PPA-GCN as input, and use 5-fold cross-validation to test performance. The results (Table 2) show that, PPA-GCN achieves the best performance among all indicators.

In order to evaluate the adaptability of PPA-GCN to new data, the genes of the new genome were classified into network nodes. The test set node of the newly assigned functional path label in the network was used as the evaluation object, and the difference between the assigned output label and the real label is directly compared. The results are shown in Table 3, which proves that the results of PPA-GCN is reliable.

Because functional pathway assignment for bacterial genomes is still in the development stage, there will inevitably be some false pathway labels on the bacterial genes. Hence, we needed to test the fault tolerance of PPA-GCN. All assigned labels were assumed to be correct. In each epoch of training, some unlabeled nodes were given random labels to also participate in the training process. Two sets of experiments were conducted. In one, a certain percentage (5–20%) of incorrectly labeled samples were added in each epoch independently, and in the other, incorrectly labeled samples were added accumulatively. The PRA without the addition of incorrect labels was taken as the standard, and the PRA after the addition of incorrect labels was divided by the standard PRA to serve as the performance indicator. The results (Figure 4, Supplementary Table S2) show that PPA-GCN can still maintain more than 75% performance with the addition of incorrect labels at a rate of up to 100% (that is, the incorrectly labeled samples compose up to 50% of the training set). Because the distribution of wrong labels is random, and the distribution of correct labels is ordered, the influence of correct labels on the training results is greater than that of wrong labels, which enhances the fault tolerance of the framework. With an increasing proportion of incorrect labels, the efficiency of the framework did not drop sharply. This result shows that PPA-GCN has strong fault tolerance.

A graph neural network can achieve excellent prediction accuracy, but it is difficult to give practical meaning to features. To evaluate the importance of the selected features, the PRAs before and after feature removal were compared (Figure 5, Supplementary Table S3). There are three important features in the PPA-GCN input: the node scale, the adjacency probability matrix and the gene synteny network.

The node scale is defined as the number of genes grouped into one node. The node scale was selected as an input feature because it can reflect the characteristics of a group of genomes. Staphylococcus has a closed pan-genome with an average node scale of 132.3, that is, an average of approximately 132 genes grouped into one node. Flavobacterium and Pseudomonas have open pan-genomes with average node scales of only 4.8 and 6.9, respectively. The node scale was one of the major observed differences between the labeled (training set) and unlabeled (test set) node sets in the gene synteny network. PPA-GCN showed no significant difference in performance when the node scale information was removed from the input (Figure 5A). The node scale has no effect on framework training, and this is beneficial for the applicability of the framework to unlabeled nodes.

The locations of genes in genomes are often specific, and the gene synteny network extracted from the same genus could reflect the intrinsic properties of the genus. The adjacency probability matrix is defined as the probability that two specific nodes can achieve a certain number of connections in a specific genome synteny network. Adding the adjacency probability matrix to the input was found to greatly improve the performance of the framework (Figure 5B). The adjacency probability matrix provides PPA-GCN with an information dissemination pattern for a specific bacterial genus in the gene synteny network.

Since the adjacency probability matrix can be used to extract synteny information patterns for specific microbial species, we wished to verify whether the gene synteny network could be replaced. Two types of random networks were designed while keeping the degree distribution constant. In one case, the arrangement of the gene positions in each sample genome was disrupted, and in the other, the positional relationships of all genomes were disrupted. Three strategies were considered for feature selection: not including any adjacency probability matrix, including the adjacency probability matrix of the newly generated network, and including the adjacency probability matrix of the real network. The training results show that (Figure 5C), regardless of which random network was used, the training performance when using a random network was much lower than that achieved using the real network. Interestingly, the true probability adjacency matrix can improve the framework training performance, while including the matrix of a random network actually impairs performance. This further shows that the adjacency probability matrix can capture specific information patterns of bacterial genomes. The gene synteny network and the adjacency probability matrix can provide the framework with different information patterns, and neither can replace the other.